Abstract
SARS-CoV-2 has had an unprecedented impact on human health and highlights the need for genomic epidemiology studies to increase our understanding of the evolution and spread of pathogens and to inform policy decisions. Most efforts have focused on international or country-wide transmission, which are unable to highlight state-wide trends.
We sequenced virus genomes from over 22,000 patients tested at Mayo Clinic Laboratories between 2020-2022 and leveraged detailed patient metadata to describe county-to-county spread in Minnesota. Our findings indicate that spread in the state was mostly dominated by viruses from Hennepin County, which contains the largest metropolis. For many counties, we found that state government restrictions eventually led to a decrease in the diversity of circulating viruses from other counties and that their complete removal in May of 2021 saw a drastic revert to levels at or greater than those observed during the months before.
We also linked over 14,000 genomes with patient risk characteristics and infection-related phenotypes from the Mayo Clinic electronic health record. We found that the genetic relationship of Omicron viruses was structured by clinical outcomes when stratifying by patient risk factor and variant of concern. However, we were unable to identify nucleotide variants that drove this association.
Genomic epidemiology has provided valuable insight into transmission, evolution, and public health surveillance of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), the cause of coronavirus disease 2019 (COVID-19). This has been feasible, in large part, to unprecedented viral genomic sequencing efforts across the globe. As of the 24th of July 2022, there are over 12.1M virus sequences in GISAID 1 and over 5.9M in NCBI Virus 2 and GenBank3. Studies that focus on localized spread such as counties within a state or province can highlight and uncover transmission events that could inform statewide surveillance and prevention efforts. However, there have been limited SARS-CoV-2 genomic epidemiology studies at this geographic level in the United States. Work by Moreno et al. 4 examined evolution and spread of SARS-CoV-2 among two counties (Dane and Milwaukee) in Wisconsin, from the start of the pandemic until the end of April 2020, based on analysis of 247 new full-length SARS-CoV-2 genomes combined with sequences in GISAID. Using this data, they derived county data on synonymous and non-synonymous single nucleotide variants (SNVs), and performed a variety of phylogenetic analyses (using Nextstrain 5 and BEAST2 6), determined R0 for their region, and examined the number and timing of introductions to the two counties and how each introduction subsequently impacted the local transmission 4. The authors were ultimately able to conclude that early transmission within Dane County was not due to its initial introduction followed by local spread, but rather multiple later introductions into the region 4. In other work, Deng et al. 7 sequenced 36 clinical samples from different counties in Northern California, sourced from the California Department of Public Health, Santa Clara County Public Health Department, and the University of California San Francisco. Their phylogenetic analysis revealed multiple California clusters including Santa Clara County, Solano County, San Benito County, as well as lineages from Washington State and Europe 7. They also identified several early notable SNVs including D614G in the Spike protein 7. Work by Müller et al. 8 examined SARS-CoV-2 introduction and spread in the State of Washington including at the county level early in the pandemic from February to July 2020 with a focus on the 614G variant and Google workplace mobility data. Alpert et al. created county-level risk maps for international importation of early Alpha variants via air transport 9. Other work such as Valesano et al. 10, Holland et al. 11, and Currie et al. 12 examined spread on college campuses.
Although these studies have provided within-state snapshots into the genomic epidemiology of SARS-CoV-2, they did not consider how localized evolution and spread within a state changed over multiple years of the pandemic with the introduction and circulation of different variants of concern such as Alpha, Delta, and Omicron. As the virus continues to evolve, such studies are needed and can inform local and state public health response by highlighting the roles of various counties on state-wide transmission. In addition, they can elucidate the impact of out of state introductions on local spread which can inform policy on travel.
SARS-CoV-2 genome distribution by county and patient characteristics
We sequenced SARS-CoV-2 genomes from genomic material isolated from clinical samples of patients tested for SARS-CoV-2 infection at Mayo Clinic Laboratories (Extended Data Figs. 1-2) over a two-year period from March 2020 to March 2022. We combined these sequences with additional genomes generated for surveillance purposes by the Minnesota Department of Health (MDH) and performed Bayesian phylodynamics to understand county-to-county spread as well as the impact and timing of introductions into the State of Minnesota (see Methods). For subgroup analysis, we leveraged the Mayo Clinic’s electronic health record (EHR) to link viral genomes to patient characteristics (demographics) and disease phenotypes (hospitalized ICU, etc.) to describe the genetic structure of variant-specific genomes in relation to clinical severity of disease.
Most of the patients from whom we collected a biological specimen and generated a SARS-CoV-2 genome resided in the State of Minnesota (96%) (Extended Data Table 1). The breakdown by gender was nearly 50/50 between males and females while 50 percent of the patients were between 18-45. Fifteen percent were under 18 while 11 percent were 65 or older.
Hennepin and Ramsey Counties drives in-state transmission
We down-sampled Minnesota genomes with county metadata proportional to their case rates for a final dataset of 4,586 genomes representing 20 counties (Extended Data Table 2). We implemented Bayesian phylodynamic models to examine transmissions in Minnesota from March 2020 to February 2022 (see Methods). We recorded Markov jumps 13 to estimate the timing of introductions and their directionality. Our analysis shows that Hennepin County, the largest population in Minnesota (which includes Minneapolis, the most populated city), drove the transmission of SARS-CoV-2 viruses in the state. This includes the formation of earlier clades for 20A, 20B, 20C, and 20G, as well as variants of concern Alpha, Delta, and Omicron. Ramsey County, the second most populated county (which includes St. Paul, the Capital and second most populated city) also contributed to early evolution and spread (such as 20G) however it had greater impact later in the pandemic including variants of concern Alpha and Delta (Fig. 1A).
A) Maximum clade credibility (MCC) tree of 4,586 SARS-CoV-2 genomes from twenty Minnesota counties. We see that Hennepin County (the root) dominated the early introductions into the other counties. We annotate the different clades by Nextclade-assigned names or VoCs 14. B) Markov jumps between Minnesota Counties as shown via a Chord diagram. The colors for both diagrams represent the counties depicted in the legend.
Markov jump estimates (Fig. 1B) show that transmission of SARS-CoV-2 within the state largely originated from Hennepin and Ramsey counties (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 bordering counties such as Anoka and Dakota. Analysis of the timing of introductions between counties shows the magnitude of Hennepin as a source for transmission throughout the pandemic (Extended Data Fig. 3, red dots).
We measured the ratio of introductions to total viral flow into and out of each county by month from March 2020 to January 2022. The top six sampled counties (Fig. 2A) contain considerably less uncertainty around the posterior mean ratio (depicted by the shaded regions of the 95% Bayesian highest posterior density in Fig2A and Fig2B) than the remaining counties (Supplemental Fig. 1). While Anoka, Dakota, Ramsey, Stearns, and Washington were all fueled by introductions early in the pandemic, they experienced different trends over the remaining time of the pandemic. Anoka, Ramsey, and (to a lesser extent) Stearns experienced a relatively equal divide of both introductions into, and viral flow out of, their counties for the remaining duration of the study period, while Dakota and Washington consistently experienced more introductions than viral flow out of their counties. Meanwhile, Hennepin County showed a drastically different trend than all others as it consistently produced more viral flow into 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 (likely driven by out of state sources [Fig. 3]). This increase in Hennepin likely contributed to a decrease in viral flow into other counties (Fig. 2) as the Alpha variant began to circulate in the state.
A) We show the posterior mean ratio and 95% Bayesian highest posterior density interval. B) Many counties contain considerably more uncertainty around the median ratio given their limited number of sequences as depicted in the ridgeline plot of 95% HPD range. In Supplemental Fig. 1, we show the illustrations for all twenty counties.
The colors correspond the location of origin. The circles represent a county or region in Minnesota. The triangles represent an outside location, USA or International. The x-axis of the scatter plot displays the posterior time of the introduction while the height is jittered for visualization only. We use an arrow to show the first introduction into the state which occurred in Hennepin County at around January 23, 2020 (from a domestic location). We used Baltic to extract introductions (migration events) along the annotated branches of the phylogeographic tree.
While Hennepin County was also contributing to in-state transmission, it was also steadily dominating the harboring of SARS-CoV-2 viruses as illustrated by the Markov reward times (Extended Data Fig. 4). Ramsey county also consistently contributed to the proportion of viruses but to a much less extent.
Restrictions eventually reduced spatial mixing of SARS-CoV-2, but their elimination saw an immediate reverse of its effectiveness
We assessed county-specific virus diversity via a normalized Shannon diversity index (Extended Data Fig. 5) that we computed based on the duration of time associated with continuous partitions of the phylogeographic tree as determined by Markov jumps 15 (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 15. Conversely, a value of 1, suggests significant spatial mixing of SARS-CoV-2 between counties15. Most counties 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 second line on May 28, 2021 indicates the end of all COVID-19 restrictions in the state 16. For many counties, this timepoint is preceded by a drop in diversity suggesting the restrictions were eventually effective in suppressing county-specific diversity. However, removal of state restrictions appears to have resulted in an immediate and sharp rise (to levels at or above before the drop) in diversity. This all precedes a spike in COVID-19 cases and the circulation of the Omicron variant into the state.
Hennepin County received the vast majority of out of state introductions
We analyzed the timing and impact of introductions into the Minnesota by adding 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 additional 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 (Extended Data Table 3). We focused on the timing and source of introductions into the state during the pandemic (Fig. 3) as estimated from our maximum clade credibility tree (Extended data Fig. 6). The earliest introduction into Minnesota was to Hennepin County from a domestic source on around January 23, 2020 (depicted with an arrow in Fig. 3). 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 17. The first international introduction occurred in Hennepin on around February 10 with another shortly thereafter on around February 13. The first two county-to-county introductions were estimated to originate from Hennepin to somewhere in Central Minnesota around February 28 and from Hennepin to Ramsey around March 1. International introductions were most abundant in Hennepin County (home to the Minneapolis/St. Paul International (MSP) airport) totaling 119 separate incidents over the two-year pandemic. The domestic (USA) introductions were most frequent to Hennepin as well (n=89) followed by Southern Minnesota (n=57), potentially from nearby states.
Evolutionary relationship of Omicron viruses structured by clinical outcome
We linked a subset of our SARS-CoV-2 genomes to patient characteristics and COVID-19-related outcomes in the Mayo Clinic electronic health record (see Methods). We first stratified genomes by variant and then the patient’s calculated risk of hospitalization based on the Monoclonal Antibody Screening Score (MASS) score 18 where a patient receiving a score of 1 or 2 was considered “mild” risk, a score of 3 was considered “medium”, and, a score of 4 or more “severe”. We labeled each taxa by outcome category (mild, moderate, severe, and death) based on a classification score from the Adaptive COVID-19 Treatment Trial (ACTT) 19 (Extended Data Table 4). For each VoC-risk group phylogeny, we measured the probability of matching clinical outcome traits to adjoining taxa in the tree20. For all Alpha and Delta scenarios, we found no relationship between the structure of the tree and COVID-19 clinical outcome. However, for our low risk (N=816) and high risk (N=390) Omicron genomes we found significant P-values for the association index (AI) metric (0.03 and 0.04 respectively) suggesting that the phylogeny (genetic relationship) of Omicron viruses, when stratifying for patient risk of COVID-19 hospitalization, is structured by clinical outcome.
While it has been demonstrated that the risk of severe clinical disease is much lower for Omicron viruses overall than previous versions of the virus 21, the possible relationship of genomes between mild and moderate clinical outcomes warrants further investigation. We studied single nucleotide variants (Fig. 4) and maximum clade credibility trees (Fig. 5) of our mild and moderate clinical cases of Omicron to identify any differences that could partially explain disease outcome among patients at the same risk of COVID-19 hospitalization. We note that both the mild (N = 1,259) and moderate (N = 84) clinical outcomes were similar in terms of known Omicron variants in the Spike region including receptor binding substitutions G339D, S373P, S375F, K417N, N440K, G446, S477N, E484A, Q493R, G496S, N501Y, and Y505H. Both groups of outcomes also had S371F and S371P mutations. In addition, there was also evidence of the R346K mutation in both (47% of mild cases and 60% of moderate cases), a known BA.1 mutation 22 which has shown to be associated with a decrease in monoclonal antibody effectiveness 23. In Fig. 5, we show the maximum clade credibility (MCC) trees of sequences from samples taken from patients with low and high risk of hospitalization and color the taxa by clinical outcome. Both show the dominance of mild cases (in red) with more moderate cases (blue) intermingled on the high-risk tree. While the test for trait association via BaTS relies on a posterior collection of trees, the MCC tree represents a single view of the genetic relationship between viruses that contributed to different clinical outcomes.
We show the results via BaTS for association index (AI). The four potential outcome traits were: mild, medium, severe, and death, however only mild and moderate cases were observed among these patients. We did not calculate AI for datasets of Alpha medium (N = 17) and high (N = 46) and Omicron medium (N=91) due to low sample sizes.
We combined sequences from patients of different risk levels for hospitalization (low, medium, etc.) that had either mild or moderate clinical cases.
We color the taxa by clinical outcome (mild, moderate, severe, and death).
Discussion
We analyzed over 22,000 new genomes of patients tested at Mayo Clinic Laboratories during a two-year period in the COVID-19 pandemic. We focused our Minnesota analysis at the county level (2nd administrative boundary) to examine the roles and relationships of virus transmission. Despite numerous efforts in genomic epidemiology, few studies have focused on county-to-county transmission in the US over most of the pandemic (including different VoCs). We expand on earlier efforts such Moreno et al. and 4 and Deng et al.7 but include multiple variants and an extensive timeframe. We found that Hennepin County was the source for most of the county-to-county introductions after its initial introduction with the virus in early 2020 from an out of state sources. We found that state restrictions eventually led to a decrease in county-specific diversity but their complete removal in May of 2021 saw a sharp increase back to levels at or before their elimination.
We leveraged the Mayo Clinic’s EHR to map and patient characteristics and clinical endpoints to over 14,000 SARS-CoV-2 genomes and found that the genetic evolution of sampled Omicron viruses in the state were structured by outcome, the majority of which were mild or moderate clinical cases. We were unable to find a specific, or a combination of, single nucleotide variants that explained the relationship of these viruses. However, the findings demonstrate the importance of linking virus sequences to patient metadata including risk factors (for confounding) and clinical endpoints arising from infection.
We note several limitations in the study including the likelihood of location-specific sampling bias. We attempted to supplement known locations of patients in our study (biased towards southeastern Minnesota) with existing sequences provided via the Minnesota Department of Health. We scaled our number of sequences to the rate of known COVID-19 cases, and, after doing so, omitted counties with a limited number of sequences. Thus, we are unable to account for virus spread from less populated areas of the state. In addition, our use of different versions of the DRAGEN pipeline over the course of our two-year study period, likely led to differences in variant frequencies across virus lineages/VoCs.
Methods
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 24. We used Illumina’s Dynamic Read Analysis for GENomics (DRAGEN) COVID Lineage software and pipeline 25 (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.
For downstream analysis and quality control, we used the following exclusion criteria including consensus sequences that were: 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 14 as bad; 3) potentially contaminated based on presence of unusual allele frequencies (< 0.9); 3) duplicate runs; 4) positive or negative controls.
Single Nucleotide Variant Annotation and Frequency Plotting
We used the Ensembl Variant Effect Predictor (VEP) 26 and a custom gff3 file containing matching gene coordinates used by Nextstrain 5, and Wuhan 1 (NC_045512.2) as the reference, to determine synonymous and non-synonymous variants. We then combined this information into a single file containing county and state information for each sample. We used R and the package dplyr 27 to extract and import the annotated variant data and the package Karyplot 28 to generate plots of synonymous and missense variants across the SARS-CoV-2 genome.
Phylogenetic analysis of SARS-CoV-2
We created three different datasets for phylogenetic analyses including a:
Minnesota counties dataset (N=4,586) that included full genome SARS-CoV-2 sequences from the 20 counties with the greatest number of reported COVID-19 cases as of February 28, 2022. We included sequences from March 2020 through February 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 (Extended Data Table 2 and Supplemental Fig. 2) and used the filter module in augur 29 to distribute (as equally as possible) our heterochronous sequences by month (Extended Data Fig. 7). 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).
International dataset (N=6,479) that contained our Minnesota county samples as well as a global representation of sequences available via GenBank as part of an open access dataset from Nextstrain (Extended Data Table 3) 30. We used the list of accessions to download sequences from NCBI Virus2. 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.
SARS-CoV-2 genomes-clinical phenotype dataset (N=14,810) of sequences with linked Mayo Clinic EHR data including COVID-19 hospitalization risk and patient outcome. We separated these data to variant-specific sets of Alpha, Delta, and Omicron. To partially adjust for confounding, for each VoC, we further divided the data by low, medium, and high risk of COVID-19 hospitalization as determined by the patient’s Monoclonal Antibody Screening Score (MASS) score18.
We aligned all sequences using Mafft 31 and soft masked problematic sites within coding and non-coding regions (UTRs) using a program described in 32,33. We created initial phylogenetic trees via Nextstrain’s augur pipeline 5 with IQTree and rooted based on Wuhan-1 (MN908947) 34. We used TempEST 35 to examine the temporal signal of our heterochronous samples and removed outliers based on visualization of root-to-tip divergence and the residuals. We used least-square dating as done in 36 to refine the root. For our Minnesota-only phylogeny, we additionally pruned the Wuhan-1 outgroup.
Phylodynamics of SARS-CoV-2 in Minnesota
Both of our Minnesota-only (0.94) and international (0.93) datasets showed a positive correlation of sampling date and genetic divergence thus suggesting a strict molecular clock 35. Since our downstream inferencing framework utilizes a rooted and non-bifurcating starting tree, we used the program di2multi in the R package ape 37 to manipulate our refined trees as necessary.
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 38,39. We used our starting trees and specified a GTR + Γ model of DNA substitution and a non-parametric Bayesian SkyGrid coalescent model for our tree prior 40. We ran our Markov-chain Monte Carlo (MCMC) simulation for 5×108 steps and sampling every 5×104 steps for our Minnesota-only dataset and 109 steps with a sampling rate of 105 for our international dataset. We checked for convergence of model parameters via Tracer v1.7.1 41 with an ideal effective sample size (ESS) threshold of 200. We also used Tracer to produce graphs of demographic reconstruction of our estimated virus population size via our Skygrid model 40 (Supplemental Fig. 3). We generated log marginal likelihoods and evaluated population growth priors via a stepping stone and path sampling procedure 42. For both datasets, our results suggested the use of the non-parametric Skygrid tree prior over a constant growth model (Supplemental Table 1).
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. For our Minnesota-only phylogeny, we populated our location trait by the county name of each taxa’s sampling location. For our international set, 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 (Extended Data Table 3). In BEAUti 43, we specified an asymmetric transmission rate matrix of K(K*1) where K is equivalent to the number of discrete locations (20 for Minnesota-only and 10 for international). We recorded Markov jumps13 between locations to estimate the timing and source of introductions and Markov rewards for the time maintained within a location. We combined posterior log data via LogCombiner v.10.4 43 as needed and used TreeAnnotator v.10.4 43 to create a single maximum clade credibility (MCC) tree after 10% burn-in. We used baltic 44 for tree visualization and to extract the timing of discrete location transitions along the branches of the MCC for our estimates of introductions. We used SpreadD3 45 to calculate the Bayes factors to identify the most parsimonious origin-destination scenarios (Supplemental Table 2).
We used two new programs, introduced in 15, as part of our phylodynamic analyses. TreeMarkovJumpHistoryAnalyzer samples from the posterior distribution of trees to collect the timing and location of each Markov jump 15. We used the output from this program to calculate the ratio of introductions to total viral flow into and out of each county as described in Lemey et al. 15 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 state15. We used the output from this program to calculate the normalized Shannon entropy metric as described in Lemey et al.15. 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 46 to calculate normalized monthly diversity metrics for each county and HDinterval 47 for the corresponding 95% highest posterior density region.
SARS-CoV-2 variant-specific genetic characterization and clinical severity
For our 14,810 linked clinical sequences that met our previously stated sequence quality inclusion criteria, we used Bayesian Tip-association Significance (BaTS) testing as introduced in 20 to examine whether our phylogenetic tree of complete SARS-CoV-2 genomes are structured by COVID-19. The null hypothesis of this test is that adjoining tips in the tree are no more likely to have the same trait of interest than by chance 20. We first created a posterior distribution of Bayesian phylogenetic trees via BEAST v1.10.4 43 and Tracer v.1.7.1 41 to check for convergence and ESS values. We wrote a Python script to format our posterior set of trees and assign tips to COVID-19 classification scores (traits) as required by BaTS. We specified a single mode of analysis with 100 randomizations.
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 title 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’s firewall.
Data availability
We have deposited the SARS-CoV-2 genomes and metadata from this study in GISAID with a list available at doi.org/10.55876/gis8.220720me. The Minnesota Department of Health sequences used in this study are available on GISAID with acknowledgments available at doi.org/10.55876/gis8.220709mv. Our GenBank international sequences were identified via the Nextstrain 30 site and obtained from NCBI Virus2.
Author contributions
MS designed the genomic epidemiology study, performed the data analysis, and wrote the manuscript. APN and BSP jointly supervised the study and providing funding, reviewed the manuscript, and approved the final manuscript. CER performed data analysis, assisted in reviewing and writing the manuscript. KL, YC, PTV and EWK implemented the DRAGEN and Nextstrain bioinformatics pipeline and variant analysis pipelines, performed variant analysis, and reviewed the manuscript. JMC, JSL, MM, SJM, VTW, NRS, SKS and EDW supervised and performed the sample extraction, viral genomic sequencing, primary analysis of the sequencing and variant data, and reviewed the manuscript. JJH managed the execution of the study and reviewed the manuscript. JCO extracted the clinical data, performed primary data analysis, and reviewed the manuscript. JDY coordinated the provision of samples for quality control and sequencing, assisted with variant interpretation, and reviewed the manuscript. XW supervised the generation of the Minnesota Department of Health genomic sequences, performed primary analysis of the MDH sequence data, and reviewed the manuscript.
Competing interests
The authors declare no competing interests.
Extended Data
Here, Mayo (N=21,669) represents new sequences generated from this study at Mayo Clinic Laboratories with a known sampling location in a Minnesota county. MDH (N=55,206), refers to a sequence available on GISAID with county metadata provided by the Minnesota Department of Health (MDH).
Branch colors indicate the assigned clade.
We limit this to the period of April 2020 – February 2022 for visualization purposes.
On the top of each panel, we show the normalized Shannon diversity index over time for each of the twenty Minnesota Counties. The shaded areas represent the 95% Bayesian highest posterior density (HPD). Below, we show the seven-day average cases over time of the 20 Minnesota Counties analyzed in this study and provided the CDC 48. The first vertical line indicates the end of lockdown in Minnesota on May 18, 2020 49. The second vertical line indicates the end of all COVID-19 restrictions in the State on May 28, 2021 16.
We show the cumulative number of confirmed COVID-19 cases on February 28, 2022, used to guide the number of sequences (5 per 1,000 cases) for the analysis. We also include population estimates per county based on US Census data 50.
We assigned each numerical score to one of four disease outcome variables.
Acknowledgements
This publication was supported by funding from the Center for Individualized Medicine-Mayo Clinic Research. Research reported in this publication was also supported by the National Institute of Allergy and Infectious Diseases (NIAID) of the National Institutes of Health (NIH) under Award Number R01AI164481 (to MS). The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. The authors would like to thank Gytis Dudas for his assistance with the Baltic Python library. The authors would also like to thank Scott Lunt and Rob Timmer for their assistance with the mforge high-performance computing environment. The authors acknowledge the laboratories that submitted the SARS-CoV-2 genome data to GISAID that were used in this work.
Footnotes
↵# Joint Senior Authors