Genomic epidemiology reveals the dominance of Hennepin County in transmission of SARS-CoV-2 in Minnesota from 2020-2022

Hennepin County consistently drives in-state transmission

We down-sampled Minnesota genomes nearly proportional to the number of COVID-19 cases per county and then added additional genomes from NCBI GenBank [3] as part of an international dataset used in Nextstrain [5] (see Methods). To address the computational burden of adding sequences to our already large dataset, we aggregated the additional samples into discrete traits International and USA and grouped counties with less sequences into areas in the state such as Southern, Central, and Northern Minnesota (Fig. 1 and Table S2).

• Download figure
• Open in new tab

Figure 1. Map of Minnesota counties included in the phylodynamic analysis.

Counties with the same color are part of a discrete region (northern, central, and southern) used in the analysis. We show an inset of a portion of the map for visualization purposes.

We implemented Bayesian phylodynamic models to examine transmissions in Minnesota from early 2020 to early 2022 (see Methods). We recorded Markov jumps [15] to estimate the timing of introductions and their directionality. After introductions from domestic and international locations, our analysis shows that Hennepin County, the most populous county which includes Minneapolis, the most populated city, drove the transmission of SARS-CoV-2 viruses in the state (Fig. 2). This includes the formation of earlier clades including 20C and 20G, as well as variants of concern Alpha and Delta (Fig. 2A). The counties in Central Minnesota contributed to spread including 20C, Alpha, and Delta, while Southern Minnesota contributed mostly to 20G.

• Download figure
• Open in new tab

Fig. 2. SARS-CoV-2 evolution and spread to and within the State of Minnesota.

A) Maximum clade credibility (MCC) tree of 6,188 SARS-CoV-2 genomes from Minnesota counties and regions as well as international locations and other domestic locations in USA. We manually annotate clades by Nextclade-assigned names or VoCs [16] using their lowercase Greek letter. We do not label less represented VoCs in our tree such as Gamma or Epsilon. Some clades such as 20A are not monophyletic in the tree but we label their most populous clade. B) Markov jumps between locations as shown via a Chord diagram. The colors for both panels represent the locations depicted in the legend. Central MN includes seven Minnesota counties: Benton, Carver, Chisago, Kandiyohi, Sherburne, Stearns, and Wright. Northern MN includes three counties: Clay, Crow Wing, and Saint Louis. Southern MN includes five counties: Blue Earth, Goodhue, Olmsted, Rice, and Scott. USA includes all states except for Minnesota. Abbreviation: MN – Minnesota, α - alpha, δ - delta, ο - omicron. We combined panel A and B and re-created the legend and text labels in Adobe Illustrator for visualization purposes.

Markov jump estimates (Fig. 2B) as shown via a Chord diagram suggest that transmission of SARS-CoV-2 within the state largely originated from Hennepin County (thick arcs and wider fragments at the outer circle). However, we also note existence of transmission back to these areas (white space between arc points and outer fragment) from nearby counties from Central Minnesota.

We measured the ratio of introductions to total viral flow into and out of each county by month from March 2020 to January 2022. A value of 1 suggests a county as solely being a “sink” (accepts SARS-CoV-2 lineages but never exports them to other counties), while a value of 0 indicates a county as solely being a “source”. Anoka, Dakota, Ramsey, Northern Minnesota, Southern Minnesota, and Washington were fueled by introductions mostly throughout the pandemic (Fig. 3). Meanwhile, Central Minnesota (outside of Hennepin and Ramsey) was dominated by introductions early in the pandemic but later in 2020 experienced brief trends of higher virus exportation. Hennepin County showed a drastically different trend than all others as it consistently acted as a source for other Minnesota counties over the nearly two-year period. However, it did experience brief periods of fluctuation such as a spike in the ratio of introductions towards the end of 2020 and early 2021, potentially driven by the dominance of out of state introductions.

• Download figure
• Open in new tab

Fig. 3. Ratio of introductions to total viral flow into and out of each discrete location by month from March 2020 – January 2022.

We show the posterior mean ratio and 95% Bayesian highest posterior density interval. We note that some locations such as Anoka, Dakota, Ramsey, and Washington have wider intervals during certain months, such as January 2022, due to a decrease in local sampling.

Low-to-intermediate spatial mixing within the State of Minnesota

We assessed county-specific virus diversity via a normalized Shannon diversity index (Fig. S3) that we computed based on the duration of time associated with continuous partitions of the phylogeographic tree as determined by Markov jumps [17] (see Methods). The index, in this context, measures the degree of spatial structure (based on counties) during the evolution and spread of SARS-CoV-2 viruses in Minnesota. A value of 0, indicates exclusive spatial structure such as an outbreak contained to only one county [17]. Conversely, a value of 1, suggests significant spatial mixing of SARS-CoV-2 between counties [17]. The counties and regions show low to intermediate (0.25 to 0.5 Shannon) spatial mixing with brief periods of waxing and waning. The two dotted vertical lines indicate changes in state-wide policy. The first vertical line indicates the end of lockdown in Minnesota on May 18, 2020 [18]. The second line on May 28, 2021 indicates the end of all COVID-19 restrictions in the state [19]. Neither of these policy decisions appeared to have significant impact on virus diversity across each individual county. Anecdotally, (looking at the trends of each graph) the changes in case counts over time does not appear to have a relationship with county-specific diversity.

Hennepin County received the vast majority of out of state introductions and was the dominant source for in-state transmission

We focused on the timing and source of introductions into the state during the pandemic (Fig. 4) as estimated from our maximum clade credibility tree (Fig. 2A). The earliest introduction into Minnesota was to Hennepin County from a international source on around January 22, 2020 (depicted with an arrow in Fig. 4). This is about one month before the first patient in the state, a man from Ramsey County (which borders Hennepin), developed symptoms and around six weeks before (March 6, 2020) the Department of Health confirmed the infection [20]. The first county-to-county introductions were estimated to originate from Hennepin to somewhere in Northern Minnesota around February 22 and from Hennepin to Washington County (also in the northern part of the state) around February 24. International introductions were most abundant in Hennepin (home to the Minneapolis/St. Paul International (MSP) airport) totaling 45 (out of 107) over the two-year period. Southern Minnesota counties were most common for domestic introductions with 19 (out of 64), potentially driven by bordering states such as Iowa and Wisconsin as well as Illinois which is nearby. Hennepin also was, by far, the most dominant source of in-state transmissions to other Minnesota locations (n=772) over the two-year period.

• Download figure
• Open in new tab

Fig. 4. Timing and source of international, domestic, and within state introductions for each discrete location.

The colors correspond to the source location. We use an arrow to show the first introduction into Minnesota which occurred in Hennepin County at around January 22, 2020 (from an international location). We used Baltic to extract introductions (migration events) along the annotated branches of the phylogeographic tree for node states with a posterior probability of ≥ 0.90. We combined the location panels into one figure in R and replaced the original month labels (x-axis) in Adobe Illustrator for visualization purposes.

RNA extraction, library preparation and next generation sequencing

From March 2020 to March 2022, we analyzed patient nasopharyngeal or mid-nasal turbinate swabs that tested positive for COVID-19 via RT-qPCR at Mayo Clinic Laboratories and had a Ct value of 28 or lower. We extracted viral RNA on the Hamilton Microlab STAR Automated Liquid Handler system (Hamilton Company, Reno, NV USA) with the use of Promega Maxwell HT Viral TNA Kit (Fitchburg, WI). We generated libraries using the COVIDSeq Test reagent kit from Illumina (San Diego, CA USA) following the manufacturer’s instructions. We sequenced the pooled libraries as 100 × 2 paired end reads using the NovaSeq SP sequencing kit and Xp 2-Lane kit with NovaSeq Control Software v1.6.0. We used the Illumina RTA version 3.4.4 for base-calling.

We de-multiplexed raw sequence data into individual sample fastq files using bcl2fastq2-v2.19.0 [21]. We used Illumina’s Dynamic Read Analysis for GENomics (DRAGEN) COVID Lineage software and pipeline [22] (versions 3.5.1,3.5.3, and 3.5.6) for reference-based alignment to Wuhan-1 (NC_045512.2), quality assessment, variant calling, and generation of consensus sequences.

We excluded sequences from downstream analysis if they met any of the following criteria, including: 1) given an overall score of fail by the DRAGEN pipeline due to having an insufficient amount of detectable viral reads; 2) given an overall quality score by Nextclade [16] as bad; 3) potentially contaminated based on presence of unusual allele frequencies (< 0.9); 3) duplicate runs; 4) positive or negative controls.

Phylogenetic analysis of SARS-CoV-2

We assembled a representative dataset (n = 6,188; Fig. S4) that included SARS-CoV-2 genome sequences from the 20 counties with the greatest number of reported COVID-19 cases as of February 28, 2022 as well as a global representation of sequences available via GenBank as part of an open access dataset from Nextstrain (Table S2) [23]. We used the list of accessions to download sequences from NCBI Virus [2]. We removed sequences less than 29K nucleotides in length as well as duplicates.

We included sequences from December 2019 (including Wuhan-1, GenBank accession MN908947) to February 28, 2022, as well as their sampling location and collection date metadata. To partially address sampling bias, we sampled at a rate of 5 sequences per 1,000 county cases and used the filter module in augur [24] to distribute (as equally as possible) our heterochronous sequences by month (Fig. S4). For each county, we attempted to include SARS-CoV-2 genomes across the two-year timeframe by supplementing our dataset with sequences provided by the Minnesota Department of Health (MDH) (and available in GISAID). The MDH sequences were produced from randomly selected samples from clinics and community testing sites. Sample Ct values were equal or below 30.

We aligned all sequences using Mafft [25] and used trimAl v1.4.rev22 [26] to remove columns which contained more than 70% gaps. We created initial an phylogenetic tree via Nextstrain’s augur tree command [5] with IQTree and rooted the tree based on Wuhan-1 (MN908947) [27]. We used TempEST [28] to examine the temporal signal of our heterochronous samples and removed one sequence as an outlier. We used augur refine and the keep-root option to modify our tree with sequence metadata. We removed additional sequences that had potential misassigned clades or produced inconsistences with the phylogenetic structure as shown in Nextstrain’s global all-time subsampled dataset [29].

Phylodynamics of SARS-CoV-2 in Minnesota

We used R package ape [30] to confirm that our starting tree was rooted and non-bifurcating in order to comply with our downstream inferencing framework. For Bayesian inference, we leveraged a pre-release of BEAST v1.10.5 (ThorneyTreeLikelihood v0.1.1) and BEASTGen v0.3 (pre-thorney) to specify a more efficient likelihood function intended for larger sequence datasets [31, 32]. We used our starting tree and a non-parametric Bayesian SkyGrid coalescent model for our tree prior [33]. We ran two Markov-chain Monte Carlo (MCMC) simulations each for 5 × 10⁸ steps and sampling every 5 × 10⁴ steps. We combined these two runs via LogCombiner v.10.4 [34] after removing 10% burn-in. We checked for convergence of model parameters via Tracer v1.7.1 [35] with an ideal effective sample size (ESS) threshold of 200. We generated log marginal likelihoods and evaluated population growth priors via a stepping stone and path sampling procedure [36]. Our results suggested the use of the non-parametric Skygrid tree prior over a constant growth model (Table S3).

We used LogCombiner to sample 1,000 trees from the posterior distribution and used this as empirical data for ancestral state reconstruction of our location traits. We specified all non-US sequences as “International” and non-Minnesota US states as “USA”. For computational efficiency, we kept the five counties with the greatest number of cases as independent locations and grouped the remaining fifteen counties into three discrete regions including Southern, Central, and Northern Minnesota (Table S2). In BEAUti [34], we specified an asymmetric transmission rate matrix of K(K*1) where K is equivalent to the number of discrete locations (n = 10 for our dataset). We recorded Markov jumps [15] between locations to estimate the timing and source of introductions and specified a MCMC of 5 × 10⁶ sampling every 5 × 10² steps. We used TreeAnnotator v.10.4 [34] to create a single maximum clade credibility (MCC) tree after 10% burn-in. We used baltic [37] for tree visualization and to extract the timing of discrete location transitions along the branches of the MCC for our estimates of introductions. We excluded transmission chains with low support of discrete origin and destination names by including only nodes with a posterior probability of ≥ 0.90 for the location state. We used SpreadD3 [38] to calculate the Bayes factors to identify the most parsimonious origin-destination scenarios (Table S4).

We used two programs of the BEAST library [39], introduced in [17], as part of our Bayesian phylodynamic analyses. TreeMarkovJumpHistoryAnalyzer samples from the posterior distribution of trees to collect the timing and location of each Markov jump [17]. We used the output from this program to calculate the ratio of introductions to total viral flow into and out of each county (number of introductions / (number of introductions + number of exports)) as described in Lemey et al. [17] as well as the visualization of the weights of pairwise transmission between counties via a chord diagram. TreeStateTimeSummarizer, which also samples from the posterior distribution of trees, notes the contiguous partitions for a given discrete state [17]. We used the output from this program to calculate the normalized Shannon diversity metric as described in Lemey et al. [17]. We used this measure to assess the level of location diversity for the viruses within each county during a specified time-period. For our analysis, we used NormShannon method in the R package QSutils [40] to calculate normalized monthly diversity metrics for each county and HDinterval [41] for the corresponding 95% highest posterior density region.

Human subjects and ethics approval

This research was conducted under approval of ethics by the Mayo Clinic Institutional Review Board and assigned a study ID IRB#: 20-005896 and entitled Large Scale Whole Genome Sequencing of SARS-CoV-2. All our published datasets contain randomly generated study IDs and no personal identifiers. All data analysis was performed behind the Mayo Clinic firewall.