SARS-CoV-2 Genomic Surveillance Reveals Little Spread Between a Large University Campus and the Surrounding Community

University Setting and Epidemic Course

The COVID-19 epidemic in southeastern Michigan evolved over the course of the Fall 2020 semester. Washtenaw County has a population of over 367,000 and is part of Michigan Public Health Preparedness Region 2 South, which includes Monroe and Wayne counties and the city of Detroit. From August 16^th through November 24^th, there were 6707 laboratory-confirmed COVID-19 cases in Washtenaw County. COVID-19 cases increased in Washtenaw County over the fall semester (Figure 1A), rising from 130 cases during the week of August 23^rd (63 daily cases per million) to 1125 cases during the week of November 15^th (405 daily cases per million).

• Download figure
• Open in new tab

Figure 1:

Case curves and sequencing density. (A) New lab-confirmed cases of COVID-19 in Washtenaw County, Michigan from the week of 2020-08-16 through 2020-11-23, displayed by day of symptom onset (as reported by MDHHS). New cases per week are shown on the y-axis and time in weeks on the x-axis. The fraction of new lab-confirmed cases in University of Michigan students is shown in green. (B) Sampling density is displayed as the fraction of new lab-confirmed COVID-19 cases with complete genome sequences (y-axis) per week during the fall term (x-axis). The fraction of student cases sequenced is shown in yellow, all Washtenaw County cases in blue, and all cases in Region 2S in violet (includes Washtenaw, Wayne, and Monroe counties). (C) Genomes sequenced in this study. The number of genomes from cases in students is shown in green and non-students in grey. (D) Of the genomes from cases in students, information on residence hall is displayed (colors).

The University of Michigan-Ann Arbor is a public university with an enrollment of over 47,000 undergraduate, graduate, and professional students. In Fall 2020, the University of Michigan-Ann Arbor opened for on-campus student residence and hybrid in-person and remote learning. In-person instruction was held from August 31^st through November 24^th. Out of 72,798 tests performed from the week of August 16^th through the week of November 22^nd, there were 2374 COVID-19 cases in students, 1064 (44.8%) of which had testing performed at Michigan Medicine laboratories. There was no syndromic case definition for obtaining COVID-19 testing at Michigan Medicine at that time, and testing was available for both symptomatic and asymptomatic cases depending on clinical circumstances. Cases in University of Michigan-Ann Arbor students constituted the majority of cases in Washtenaw County by the end of September 2020 (Figure 1A). In response to this spike, several mitigation efforts were implemented, including stay-in-place orders targeted to specific residence halls, mass testing events, and a broad stay-in-place order for all undergraduates from October 20^th –November 3^rd. The fraction of campus-associated cases in the county declined after mid-October 2020 during a wave of new infections in Washtenaw County in November 2020. The shift in case incidence prompted us to investigate whether the November wave of infections was generated by spillover from outbreaks among students.

Genomic Surveillance in Southeastern Michigan

We processed all SARS-CoV-2 positive samples tested at UHS or the University of Michigan Clinical Microbiology Laboratory collected between August 16^th –November 24^th for genome sequencing. November 24^th corresponded to the end of in-person instruction for the semester and the peak of the November surge of new cases in Washtenaw County. We assembled complete SARS-CoV-2 genomes from 1659 individuals, representing a median of 24% of the confirmed cases per week from Washtenaw County. For comparison, these genomes represent a median of 4.5% of cases per week from the entire Michigan Region 2 South (Figure 1B). We did not have access to SARS-CoV-2 specimens not tested at Michigan Medicine sites, which limited our sampling density (Figure S1A). Despite this, we sequenced 468 complete genomes from University of Michigan-Ann Arbor students, representing 20% of cases from August 16^th through November 24^th (Figure 1C, S1B). Of these 468 students, we were able to determine on-campus residences for 131; the on or off campus residence of the remainder are unknown (Figure 1D). The genomes presented here consisted of several different viral clades, mostly Nextstrain clades 20A, 20C, and 20G (Figure S1C). We did not detect any variants of concern that were circulating in other countries at the time, including B.1.1.7, P.1, or B.1.351.

Transmission Introductions in Students at the University of Michigan-Ann Arbor

We used phylogenetic analysis to characterize the influx of viruses into the student population. We generated a time-calibrated maximum likelihood phylogenetic tree with our sequenced genomes (n = 1659) and additional contextual genomes (See Methods; Figure 2A, S2A). This sequence alignment was a subsampling of data available on GISAID, weighted towards genomes with genetic similarity to those in Michigan. To optimize the inference of transmission lineages in southeastern Michigan, we included all available genomes from Michigan on GISAID that were collected from July – December 2020. We inferred introductions into the student population by discrete trait ancestral reconstruction (see Methods). We inferred traits of ancestral nodes on this time-calibrated phylogenetic tree using a binary model of student vs. non-student, as has been performed in related genomic surveillance studies (Alpert et al., 2021; Lemieux et al., 2020; Müller et al., 2021; Plessis et al., 2021). Genomes from students that shared the same ancestral “student” node were considered part of the same introduced transmission lineage. Genomes from students that were not preceded by a “student” node were considered singleton introductions. We performed the same analysis using a total of ten random subsamples of non-Michigan genomes to verify that our results were not substantially biased by contextual genomes (Lemieux et al., 2020). Due to the large variability in genomic surveillance sampling density from other states in the USA and other countries, we did not attempt to infer the exact geographic sources of transmission introductions into the student population.

• Download figure
• Open in new tab

Figure 2:

Introductions into the student population. (A) Time-calibrated maximum likelihood phylogenetic tree of 7149 SARS-CoV-2 genomes, including 1657 genomes sequenced for this study. Month is shown on the x-axis. Genomes from students are shown in red, other genomes from Michigan in light blue, other genomes from North America in dark blue, and global genomes in grey. (B) Counts of genomes from students (y-axis) per week during the fall term (x-axis) by the inferred transmission lineage group. Singleton introductions are shown in dark blue, genomes from Cluster A in yellow, genomes from Cluster B in violet, and genomes from smaller clusters (2 – 8 individuals) in teal. (C) Inferred transmission lineages in students, with genomes from each lineage or cluster (y-axis) shown as points with time on the x-axis. Only inferred transmission lineages with ≥5 students are shown. Point colors reflect on-campus residences, if known, as in Figure 1 (unknown residence is black).

Using this approach, we inferred 203 distinct transmission introductions into students (Figure S2B). There were two large transmission lineages in students, which we denote as Cluster A (n = 115 students) and Cluster B (n = 73 students). These lineages were the predominant source of cases in students during the middle of the semester, representing >50% of genomes from students from the week of September 20^th through the week of October 18^th. The rest of the introductions were singletons (n = 171) or small clusters of 2 – 8 students (n = 30 introductions). Small transmission lineages (2 – 8 students) were often short in duration, lasting a median of 3.5 days (IQR 1 – 8 days). We do not have information on social contacts of students within these transmission lineages. We detected singleton cases over the entire course of the Fall 2020 semester. The frequency of singletons increased during the latter half of the semester (Figure 2B), indicating new introductions into the student population rather than spread from older transmission clusters, such as Cluster A and Cluster B. These data suggest a shifting epidemiology in the student population across the semester, characterized by two dominant lineages during the early and middle portions of the semester followed by many small introductions later in the semester as community incidence increased.

We further investigated the two largest inferred lineages, which indicate sustained local transmission of SARS-CoV-2 within the student population. Genomes from Cluster A were part of Nextstrain clade 20B and PANGO lineage B.1.1.304 (TMRCA August 26^th, 95% CI August 3^rd –September 26^th). Genomes from Cluster B were part of Nextstrain clade 20C and PANGO lineage B.1.593 (TMRCA July 29^th, 95% CI July 3^rd –August 28^th). We detected cases from Cluster A from October 8^th through November 14^th and Cluster B from September 11^th through October 23^rd (Figure 2C). These clusters were genetically distinct from other lineages in our data and identified in all ten independent subsamples (Figure S2C). Cluster A was defined by mutations G4207T, C11471T, G24816C, G25906T, C28399A, and Cluster B was defined by G3013A, T6769C, G12243A, T15567C, C28435T. PANGO lineages B.1.1.304 and B.1.593 are rare in the GISAID database (0.07% and 0.02%, respectively, out of 532,715 genomes from the USA as of June 7^th, 2021). Among all available genomes from Michigan, these lineages are almost completely unique to this study. Outside of the genomes presented here, B.1.1.304 has been detected in Michigan three times (238 in the USA), and B.1.593 has not been detected elsewhere in Michigan (only twice in the USA). The rarity of these PANGO lineages suggests that these outbreaks among students did not spark significant transmission in southeastern Michigan or in other parts of the country.

We examined the on-campus residence locations for students in these transmission lineages. There was no obvious association with a single residence hall in any transmission lineage in students (Figure 2C). Of the 131 individuals with known on-campus residences, 31% were singletons (n = 40) and 56% (n = 73) were part of Cluster A or B. Individuals from many different residence halls across the University of Michigan-Ann Arbor campus were present in Cluster A and B, including students from nine and seven campus residences, respectively. Of the 188 students in Cluster A and B, 61% (n = 115) did not have an identifiable on-campus residence and likely lived off-campus (Figure 2C, black points). Besides Cluster A and B, there was only one other student lineage that had more than one individual from the same campus residence (2 of 3 students in the lineage were from the same residence). The first individuals detected in Cluster B resided in the same residence hall (Figure 2C), potentially reflecting transmission in a congregate setting. However, without detailed contact tracing to complement the genomic data, it is unclear whether clusters A and B originated in specific residence halls and then spread further among students, or whether the clusters originated in off-campus gatherings. Overall, these data demonstrate local transmission among the student population with intermixing among students from multiple on-campus residences and students residing off-campus.

Limited Spillover from Student Clusters into the Broader Community

Clusters A and B, defined by PANGO lineages B.1.1.304 and B.1.593, waned by mid-November when new COVID-19 cases increased in Washtenaw County (Figure 3), which suggested that these large, student-dominated lineages did not significantly contribute to the rise in community cases. To further examine the extent of COVID-19 spread from students into the broader southeastern Michigan community, we used the same ancestral trait reconstruction method as above, applying a binary student/non-student model. We enumerated genomes from non-students that descended from an ancestral node that was inferred as “student.” There were very few non-students in the community with genomes that descended from student-associated lineages (n = 53), and most were nested within either Cluster A (n = 18) or Cluster B (n = 11, Figure 4A). A total of 24 non-students descended from the rest of the inferred transmission lineages (2 – 8 students each). Out of the 1191 genomes we sequenced from non-students, 96% (n = 1138) were not genetic descendants of detected clusters in students (Figure 4B). There was a positive association between the size of a cluster among students and the number of non-student descendants (Figure 4C). This means that although the total number of non-student descendants was low, larger transmission lineages among students resulted in more spread to the community than smaller lineages. The median age of non-students descending from Cluster A and B was 47 years (IQR 20 – 61 years. We do not have epidemiologic information on the association of these individuals with students, so it is not possible to determine the circumstances of transmission. It is also possible that student status was misclassified for these individuals. This would result in an over-estimation of the amount of spread from outbreaks among students. It is also possible that some of these individuals were campus faculty or staff and therefore had differential exposure compared to the broader community.

• Download figure
• Open in new tab

Figure 3:

All genomes from Michigan available on GISAID from July 2020 – January 2021. Count of genomes is shown on the y-axis with time on the x-axis (binwidth = 1 day). Genomes from lineage B.1.593 are shown in blue, genomes from B.1.1.304 in yellow, and genomes from all other lineages are shown in grey.

• Download figure
• Open in new tab

Figure 4:

Spillover from student-associated clusters. (A) Maximum likelihood phylogenetic trees of student Cluster A (top) and Cluster B (bottom). Time-calibrated trees are displayed on the left and divergence trees on the right. Tip colors reflect genomes from students (light blue) and non-students (red). (B) Bar plot of the number of non-student genomes sequenced (y-axis, n = 1191) per week over the fall term (x-axis). Genomes that are derived from inferred “student” nodes are shown in red and genomes not derived from “student” nodes in dark blue. (C) For each non-singleton transmission lineage in students, the number of non-student descendants is shown (y-axis) by the number of students in the cluster (x-axis). A linear regression line is shown.

To further assess the spread from students to the community, we compared ancestral state models with different transition state matrices. Discrete trait reconstruction models that allowed for different transition rates between students and non-students exhibited better fit to the data compared to equal rates models (AIC 3228 vs. 3951, respectively), but a model that did not allow transitions from students to non-students exhibited a worse fit compared to the all-rates-different model (AIC 3805 vs. 3228, respectively). Overall, these results indicate that spillover from student-associated outbreaks into the broader community was infrequent and did not significantly contribute to the case increase in Washtenaw County in November and December 2020.

Research Ethics and Sample Sources

We obtained residual upper respiratory diagnostic specimens that were positive for SARS-CoV-2 from the University of Michigan Clinical Microbiology Laboratory and the University Health Service. Use of residual SARS-CoV-2 positive specimens and collection of student status and on-campus residence were approved by the Institutional Review Board at the University of Michigan (Protocol HUM185966).

Genome Amplification and Sequencing

We extracted RNA from nasopharyngeal specimens with the Invitrogen PureLink Pro 96 Viral RNA Mini kit, using 140 μL of input sample and eluting into 50 μL. We reverse transcribed RNA with SuperScript IV (ThermoFisher) by adding 1 μL of random hexamers and 1 μL of 10 mM dNTP to 11 μL of RNA extract, heating at 65 °C for 5 min, and placing on ice for 1 min. Then we added a reverse transcription master mix (4 μL of SuperScript IV buffer, 1 μL of 0.1M DTT, 1 μL of RNaseOUT RNase inhibitor, and 1 μL of SSIV enzyme) and incubated at 42 °C for 50 min, 70 °C for 10 min, and held at 4 °C. We amplified SARS-CoV-2 cDNA in two pools using the ARTIC Network version 3 primers. We used Q5 Hot Start High-Fidelity DNA Polymerase (NEB) and the following thermocycler protocol: 98 °C for 30 s, then 35 cycles of 98 °C for 15 s, 63 °C for 5 min, and final hold at 4 °C. We combined PCR products of the two primer pools in equal volumes for a given sample and purified the DNA product with 1X volume of AMPure beads (Beckman-Coulter). We prepared libraries for Illumina sequencing with the NEBNext Ultra II DNA Library Prep Kit (NEB) according to the manufacturer’s instructions. We pooled barcoded libraries for a batch of 96 samples by equal volume and gel extracted the band at 400-500 bp to remove adapter dimers. We quantified pooled libraries with the Qubit 1X dsDNA HS Assay Kit (ThermoFisher) and sequenced on an Illumina MiSeq (v2 chemistry, 2x250 cycles) at the University of Michigan Microbiome Core facility.

Analysis of Raw Sequence Reads

We mapped reads to the Wuhan/Hu-1/2019 reference genome (GenBank MN908947.3) with BWA-MEM (Li, 2013). We trimmed the ARTIC amplification primer sequences with iVar 1.2.1 (Grubaugh et al., 2019). We determined consensus sequences with samtools mpileup and iVar 1.2.1. We assigned the consensus as the base with >50% frequency and placed an N at positions covered by fewer than 10 reads. Genomes with 29000 or more unambiguous bases (> 97%) were considered complete and used in downstream analysis. The consensus genomes are available on GISAID as listed in Supplemental Table 1.

Case Definitions and Metadata

Of the samples from August 16 – November 24^th that yielded complete genomes, we used unique identifiers associated with each sample to obtain a single genome per individual. We used these identifiers to determine which individuals were students at the University of Michigan-Ann Arbor in the fall semester (undergraduate, graduate, or professional). We queried two databases hosted by the University Health Service to identify students and housing status. We identified 468 individuals who were listed as students in either or both databases. There were three individuals with conflicting classification; we investigated these further to determine the most accurate classification. Students with a campus residence hall listed in either database were considered “on-campus,” and the rest were considered “unknown.” We obtained case metrics in students from the University of Michigan COVID-19 Dashboard (https://campusblueprint.umich.edu/dashboard/, accessed on March 4^th, 2021). We obtained data on new confirmed cases in Washtenaw County and Region 2S from the Michigan Department of Health and Human Services website (accessed March 4^th, 2021).

Phylogenetic Analysis and Discrete Trait Reconstruction

To provide phylogenetic context, we obtained a subsample of genomes on GISAID using the augur pipeline (Hadfield et al., 2018). We excluded genomes collected after 2020 and genomes with lengths less than 27000 bases. We aligned genomes with MAFFT as implemented in augur. We performed subsampling with priority on genomes with genetic similarity to Michigan sequences using the following rules: 100 genomes per month in Michigan, 5 genomes per month per division in the USA outside of Michigan, 1 genome per country per month in North America outside of the USA, and 2 genomes per country per month outside of North America. If not captured in the subsampling scheme, we added back the genomes presented here (n = 1659) as well as all other genomes from Michigan collected in August – December 2020. This resulted in an alignment with 7174 sequences, including 3318 genomes from Michigan, 2157 genomes from the USA outside of Michigan, and 1699 global genomes.

We masked the 5’ and 3’ ends of the alignment along with other sites commonly affected by sequencing errors and homoplasies. We inferred a maximum likelihood phylogeny with IQ-TREE with a GTR model and 1000 ultrafast bootstraps (Nguyen et al., 2015). We then used TreeTime to generate a time-scaled phylogeny rooted on Wuhan/Hu-1/2019 with a clock rate of 0.0008 with a standard deviation of 0.0004 substitutions/site/day (Sagulenko et al., 2018). We flagged and removed genomes that exceeded 3 interquartile ranges from the root-to-tip regression (n = 25, including two genomes we generated). We generated a new time-scaled phylogenetic tree with IQ-TREE and TreeTime as above with this filtered alignment (n = 7149). We used TempEst to fit and plot a root-to-tip regression of divergence over time (Rambaut et al., 2016). We used this tree as the basis for discrete-trait ancestral state reconstruction (ASR) using a binary model of student vs. non-student. We performed maximum likelihood ASR with TreeTime with a sampling bias of 2.5. We also performed Bayesian ASR using BEAST 1.10.4 using a fixed tree topology derived from TreeTime (Alpert et al., 2021; Drummond and Rambaut, 2007). We summarized the ASR results in R version 4.0.2. We inferred introductions from non-student to student populations at nodes with a BEAST probability of > 0.9. From these inferred “student” nodes, we counted the number of descendant student and non-student tips. Introduction nodes leading to only one student case were considered singletons, and introductions leading to two or more student cases were termed transmission lineages. To determine whether the contextual genomes were biasing our results, we generated a total of ten random subsamples of the global data using the schema described above and analyzed the data in the same manner (Figure S2C). Contextual genomes from GISAID used in these analyses are listed in Supplemental Table 1. We compared ancestral state models with different transition matrices using the functions implemented in the “castor” R package (Louca and Doebeli, 2018).

Data and Code Availability

The consensus genomes we generated for this study are publicly available on GISAID. Accessions for these sequences as well as other genomes from GISAID used in these analyses are listed in Supplementary Table 1. Analysis code for the generation of consensus sequences and phylogenetic analysis are available on GitHub at the following link: https://github.com/lauringlab/SARS2_Fall_2020.