Baseline Sequencing Surveillance of Public Clinical Testing, Hospitals, and Community Wastewater Reveals Rapid Emergence of SARS-CoV-2 Omicron Variant of Concern in Arizona, USA

ABSTRACT The adaptive evolution of SARS-CoV-2 variants is driven by selection for increased viral fitness in transmissibility and immune evasion. Understanding the dynamics of how an emergent variant sweeps across populations can better inform public health response preparedness for future variants. Here, we investigated the state-level genomic epidemiology of SARS-CoV-2 through baseline genomic sequencing surveillance of 27,071 public testing specimens and 1,125 hospital inpatient specimens diagnosed between November 1, 2021, and January 31, 2022, in Arizona. We found that the Omicron variant rapidly displaced Delta variant in December 2021, leading to an “Omicron surge” of COVID-19 cases in early 2022. Wastewater sequencing surveillance of 370 samples supported the synchronous sweep of Omicron in the community. Hospital inpatient COVID-19 cases of Omicron variant presented to three major hospitals 10.51 days after its detection from public clinical testing. Nonsynonymous mutations in nsp3, nsp12, and nsp13 genes were significantly associated with Omicron hospital cases compared to community cases. To model SARS-CoV-2 transmissions across the state population, we developed a scalable sequence network methodology and showed that the Omicron variant spread through intracounty and intercounty transmissions. Finally, we demonstrated that the temporal emergence of Omicron BA.1 to become the dominant variant (17.02 days) was 2.3 times faster than the prior Delta variant (40.70 days) or subsequent Omicron sublineages BA.2 (39.65 days) and BA.5 (35.38 days). Our results demonstrate the uniquely rapid sweep of Omicron BA.1. These findings highlight how integrated public health surveillance can be used to enhance preparedness and response to future variants.

representative outlook of the virus lineages circulating in a geographic region. Here, we investigated the emergence of the Omicron variant of concern in Arizona by leveraging baseline genomic sequence surveillance of public clinical testing, hospitals, and community wastewater. We tracked the spread and evolution of the Omicron variant as it first emerged in the general public, and its rapid shift in hospital admissions in the state health system. This study demonstrates the timescale of public health preparedness needed to respond to an antigenic shift in SARS-CoV-2.
KEYWORDS emergence dynamics, hospital-associated mutations, Omicron variant, SARS-CoV-2, wastewater surveillance S evere acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is the etiological agent of the Coronavirus Disease 2019 (COVID- 19) pandemic that emerged in December 2019 (1,2). SARS-CoV-2 has evolved variants of concern (VOC) harboring mutations that have increased transmissibility, and potential to mitigate therapeutics, vaccines, and diagnostics (3)(4)(5)(6)(7)(8)(9)(10). In particular, Delta and Omicron VOCs were the first lineages to rise to become dominant circulating variants worldwide. The Delta variant evolved mutations that confer partial escape from neutralizing monoclonal and polyclonal antibodies, as well as exhibit characteristics of increased disease severity (4,6,11). The Omicron variant has increased transmissibility, considerable escape from neutralizing antibodies, and reduced disease severity relative to Delta (8,12,13). Omicron accrued at least 55 mutations relative to the initial Wuhan reference sequence, 30 of which are in the Spike gene. In contrast, Delta had 34 mutations relative to the Wuhan sequence, 8 of which are located in the Spike (14). The combination of Spike mutations in the Omicron variant resulted in a major immune-evasive antigenic shift in SARS-CoV-2 (15,16). Although BA.1 and BA.2 lineages have both been designated the Omicron VOC, phylogenetic studies support that BA.1 and BA.2 branched off independently to form monophyletic (13,17). As the global transmission of SARS-CoV-2 continues, selective pressures exerted by host immunity and vaccines are among the factors that can drive the adaptive evolution of future variants (18,19).
Next-generation sequencing provides important insights into the evolutionary trajectory of circulating SARS-CoV-2 lineages. Genomic sequencing data can be used to improve diagnostic assays (20,21) and vaccines (22,23) and inform outbreak investigations. Although studies comparing COVID-19 clinical outcome typically stratify by SARS-CoV-2 lineage infections, much less is known about mutations associated with severe clinical outcomes (24). Baseline sequencing provides a representative overview of the virus lineages circulating in a geographic region. This is done by random selection of positive diagnostic specimens in a manner representative of the community demographic and clinical disease. In comparison, targeted sampling efforts such as selecting for S-gene target failure specimens or localized outbreak investigations can bias the representation of lineages circulating (25,26). Hence, sequencing surveillance plays an integral role in enhancing public health efforts to the pandemic.
Wastewater-based epidemiology (WBE) is an effective public health approach to monitoring infectious diseases in the community (27). The advantages of WBE in the COVID-19 pandemic include estimating disease burden despite decreased clinical testing rates. Moreover, presymptomatic and asymptomatic infections may not be tested (28)(29)(30). SARS-CoV-2 wastewater sequencing reflects lineages circulating within a community (31). SARS-CoV-2 wastewater-based genomic surveillance can be used in early detection of emerging variants and virus spread that might be missed by clinical surveillance (32)(33)(34)(35). Thus, WBE and clinical surveillance are complementary approaches to public health efforts.
Understanding the dynamics of how new SARS-CoV-2 variants emerge and sweep across the population is of critical importance to formulating informed public health responses to the ongoing pandemic. Here, we investigated the emergence of the SARS-CoV-2 Omicron variant through baseline sequencing surveillance of public clinical cases (n = 27,071) and hospital inpatient cases (n = 1,125) in Arizona at daily resolution between November 1, 2021, and January 31, 2022. We show that Omicron cases presented at hospitals 10.51 days after its emergence in community cases. We identified 4 nonsynonymous mutations in Omicron significantly associated with hospitalization compared to community cases. Wastewater genomic surveillance of Omicron showed a near-synchronous sweep in the community. We demonstrate that in Arizona, Omicron BA.1 rose to become the dominant variant within 17.02 days-more than 2fold faster than the prior Delta variant (40.70 days), or subsequent Omicron sublineages BA.2 (39.65 days) and BA.5 (35.38 days).

RESULTS
Omicron BA.1 displacement of Delta variant leading to surge in COVID-19 cases. As part of public health surveillance, we conducted SARS-CoV-2 baseline sequencing surveillance of SARS-CoV-2 positive diagnostic specimens collected between November 1, 2021, and January 31, 2022, in Arizona. Saliva specimens from voluntary, self-scheduled COVID-19 testing were submitted to the clinical testing laboratory from local and public testing sites across the state (Fig. 1A). In addition, positive specimens from three hospitals in the Arizona health system (Dignity Health, Phoenix Children's Hospital, and Valleywise Health Medical Center) were submitted for sequencing surveillance. Next-generation sequencing was performed on all positive specimens without sampling bias (i.e., baseline surveillance; not targeted by prescreening for variants, travel history, disease severity, etc.). SARS-CoV-2 genome sequences were generated from 27,071 community cases (average 294 specimens per day) and 1,125 from hospital inpatient cases (average 12 specimens per day) ( Table 1). Most of the sequenced specimens were collected from Maricopa County (n = 19,655; 69.7%), with Coconino (n = 4,419; 15.7%) and Pima (n = 1,830; 6.5%) counties having the second and third highest density of sequenced specimens, respectively (Fig. 1B). All other Arizona counties had fewer than 400 sequenced specimens, and 1,618 (5.7%) specimens were collected from patients that registered residence outside Arizona despite sample collection and testing occurring within Arizona. The prevalence of Omicron in baseline surveillance increased from 3.3% of sequenced cases the week of December 6, 2021, to 89.0% of sequenced cases the week of December 27 th , 2021 (Fig. 1C). Distinct viral genomic characteristics of hospital Omicron cases. In Arizona, the "Omicron surge" in COVID-19 cases was reflected in a marked increase in COVID-19 hospitalizations ( Fig. 2A). To compare the temporal emergence of Omicron in the community to hospital admissions, we applied nonlinear regression modeling to the daily Omicron case proportions in each cohort, respectively, and determined the date at which Omicron accounted for 50% of case proportions. We found that Omicron first emerged in community cases, followed by hospital cases 10.51 days later (Fig. 2B). To assess whether there were differences in specimen viral loads, we compared the diagnostic CT values between community and hospital cases. Although there was no statistical difference in CT values between Delta variant specimens, the CT values of Omicron hospital specimens were significantly lower than public specimens (Mann-Whitney U test, P , 0.0001) suggesting that hospitalized cases had a higher viral load (Fig. 2C). While hospital rates are higher in unvaccinated individuals (36), we reasoned that there might also be specific mutations associated with hospital cases compared to community population. We identified seven nonsynonymous mutations in the Delta variant nsp3, nsp6, nsp13, S, and ORF8, that were significantly associated with hospital specimens (Fig. 2D, Fisher's exact test with Benjamini-Hochberg multiple testing correction; Table 2). We identified four mutations in the Omicron variant nsp3, nsp12, and nsp13 that were significantly enriched in hospital specimens (Fig. 2D, Fisher's exact test with Benjamini-Hochberg multiple testing correction; Table 3).
Network analysis models intracounty and intercounty transmission patterns. We sought to understand the transmission dynamics of the Omicron's emergence at the population scale. Existing tools for genomic sequence network analyses do not efficiently scale to large numbers of sequences and are particularly computationally intensive for the large genome size of SARS-CoV-2 (29.9 kb) (37). Therefore, we developed a scalable sequence network analysis methodology (described in Methods). In brief, genome sequences are aligned to the Wuhan1 reference, following which pairwise genetic distances are calculated between all samples and converted to an edge list for network visualization (Fig. 3A). A consideration of this approach is missing data from individuals who did not test and specimens that were not sequenced. For example, single source individuals absent from the network that seeded infections in multiple counties may appear as intercounty transmissions in the network. Hence, to refine transmission relationships, connections (edges) between specimens (nodes) were required to be within 10 days of collection date and only the shortest source-edge is maintained. In other words, the network showed cases in relation to their most-recent potential source of transmission. This resulted in a network from 17,734 Omicron sequence nodes and inference of 21,145 connection edges within the network (average neighbor count of 2.4). The largest Omicron network cluster contained 11,258 sequences ( Fig. 3B). To validate the methodology, we assessed the network by overlaying PANGO lineage classifications and found that clusters grouped by sequence similarity (Fig. S1). We identified large plausible transmission clusters originating from and localized to Maricopa and Coconino counties. Our model suggested that substantial intracounty transmissions occurred within each county, particularly in Maricopa and Coconino counties (Fig. 3C). Further, intercounty transmission in Pima and Pinal counties could feasibly be traced to cases from Maricopa County. Overall, even when biased toward connections with the shortest physical distance, this model illustrates that robust intercounty transmission contributed in part to the spread of the Omicron variant. SARS-CoV-2 wastewater sequencing surveillance corroborates synchronous sweep of Omicron. Wastewater monitoring of SARS-CoV-2 is an effective surrogate of community transmission. We performed wastewater SARS-CoV-2 surveillance on 370 samples collected from 14 sites across the greater Tempe area (population of approxi-  SARS-CoV-2 Omicron Emergence Dynamics mBio mately 700,000) in Arizona (Fig. 4A). SARS-CoV-2 virus RNA was detected by qRT-PCR assays performed on the wastewater specimens during the study window (Fig. 4B). Genomic sequencing indicated that the Delta variant was the dominant variant throughout November and the first week of December 2021. Omicron sequences were first detected in wastewater samples obtained from two sites on December 9th, 2021. We found that the Omicron variant rapidly displaced the Delta variant to become the dominant variant typically within 1 to 2 weeks (Fig. 4B). As a control, we performed sequencing on triplicate independent wastewater extractions demonstrating consistency in wastewater variant calls by Freyja (Fig. 4C). To corroborate Freyja's variant analysis, we compiled a list of variant specific mutations from the five most abundant sublineages of each variant observed in community samples during the investigation period. Fifteen Delta-specific and 17 Omicron-specific SNVs were found across the genome (Fig. 4D). Wastewater sequencing results were queried for the proportion of SNVs found and the mutation frequency at each locus ( Fig. 4E and F). At the start of the investigation period, Delta-specific mutations are found in high prevalence and abundance, but diminished as the Omicron specific mutations increase. We next performed molecular validation of wastewater detection of the Omicron variant using a digital PCR (RT-dPCR) assay by leveraging a 9-nucleotide deletion Taken together, these data supports a date of initial introduction of Omicron no later than December 9th. Further, the molecular RT-dPCR assays and genomic sequencing of wastewater samples applied in tandem expanded the capacity to detect Omicron in wastewater than either method alone. Omicron BA.1 displaces competing variants and rises to dominance faster than Delta and subsequent Omicron lineages. To benchmark the rate of emergence of variants to dominance, we compared Omicron BA.1 to the prior Delta, and subsequent Omicron lineages. Here, we expand the scope of genomes analyzed to include all SARS-CoV-2 baseline sequencing surveillance GISAID submissions originating from our lab between January 1st, 2021, and August 30 th , 2022. Using nonlinear regression models of baseline sequencing surveillance data, we calculated the time it took each variant to increase from 10% to 90% case proportion in nonlinear regression models of baseline sequencing surveillance data. We show that Delta emerged and displaced other competing variants (Alpha, Beta, Gamma, Epsilon) in Arizona in 40.70 days (Fig. 6A). Omicron BA.1 emergence occurred in 17.02 days -2.39-fold faster than Delta. Omicron BA.2 subsequently displaced BA.1 within 39.65 days-an emergence period comparable to Delta. At the time of manuscript submission, Omicron BA.5 lineages only reached a plateau of around 88% prevalence due to cocirculating Omicron BA.4. When adjusted to a window of 10% to 80% prevalence, Omicron BA.5 lineages emergence was 35.38 days, demonstrating similarity to Delta and Omicron BA.2. The hillslope value of the nonlinear regression curves provides an alternative measure of each variant's rate of emergence accounting for differences in the prevalence plateau of each variant. Omicron BA.1's hillslope (n H = 0.12) is more than twice steeper than that of Delta (n H = 0.048), BA.2 (n H = 0.051) or BA.5 (n H = 0.053). The 95% confidence interval for the hillslope of Delta, BA.2, and BA.5 overlap whereas the 95% confidence interval for Omicron BA.1's hillslope excludes that of the other variants (Fig. 6A). These results corroborated the accelerated rate of emergence of Omicron BA.1. Since the introduction of Omicron BA.1 coincided with a period of seasonal holidays in the US, it is possible that social gatherings might have been an extraordinary driver of BA.1's emergence. However, it is also possible that inherent characteristics of Omicron BA.1, such as an immune-evasive antigenic shift in Spike, might be responsible for its rate of emergence. We reasoned that an analysis of the variants from a region that is different in geography and demographics may provide insight. Hence, we compared Arizona genome sequences to those from South Africa. Indeed, we found that the rate of

SARS-CoV-2 Omicron Emergence Dynamics mBio
Omicron BA.1 emergence (13.34 days) was 3.85-fold faster than Delta (51.34 days) suggesting that Omicron BA.1 was uniquely poised for a rapid emergence (Fig. S2). Overall, Omicron BA.1 exhibited an emergence in Arizona and other locations worldwide at a timescale unmatched by competing lineages.

DISCUSSION
In summary, we traced the prominent emergence of SARS-CoV-2 Omicron variant through the community and health system in Arizona. The rate of emergence of Omicron BA.1 (17.02 days) was more than 2.3-fold faster than prior Delta variant and subsequent Omicron lineages BA.2 and BA.5. Our results reveal the uniquely rapid sweep of Omicron BA.1 in Arizona.
Since SARS-CoV-2 genomic sequencing is typically tied to specimens from clinical laboratory-based and point-of-care testing, changes in testing behavior can impact the ability to survey for new variants. Changes to social behavior, funding reimbursements for testing and perceived risks to COVID-19 have led to decreased testing rates (38). COVID-19 self-tests are increasingly used which has made individual health management more accessible (39). However, self-test results are not reported to public health agencies leading to under ascertainment of COVID-19 cases. Further, there is no formal mechanism for genomic sequencing on individuals who self-test. Conversely, individuals with more severe disease are likely to seek care from health care providers. Thus, implementation of public health SARS-CoV-2 sequencing surveillance needs to adapt accordingly.
WBE can effectively estimate COVID-19 community case prevalence (40), infer the variants in communities (32), and identify the presence of cryptic mutations (32,35). Our wastewater sequencing and dPCR molecular assays corroborated the synchronous introduction of Omicron in the community. Interestingly, wastewater surveillance studies in AB, Canada found that the introduction of Omicron was earlier in communities closer to dense population centers and tourist attraction locations which attract international travelers (41). Similarly, the time of initial establishment and subsequent rise of Omicron in MA, USA was found to be earlier and more rapid at institutes of higher education likely owing to dense dormitory style housing and differences in social behavior (42). These findings illustrate the value of deploying wastewater surveillance strategically with on the ground local context and social structure.
Specific mutations in Spike and nonstructural genes were significantly associated with hospital cases in Omicron and Delta variants. Omicron mutations associated with hospital genomes were encoded in the Nsp13, Nsp3, and Nsp12 proteins. Nsp12 and Nsp13 are involved in viral RNA synthesis (43). Nsp3 is a viral transmembrane protein involved in the formation of double-membrane vesicles during SARS-CoV-2 replication (44). While mutations in Nsp12 and Nsp13 encoded proteins may influence viral load via involvement in viral RNA synthesis, we found no statistically significant difference in viral load between hospital cases possessing Nsp13 S301, Nsp3 V1069I, Nsp3 L1244F, or Nsp12 V257 and hospital cases with the reference amino acid at that locus (Fig. S3). No Omicron Spike protein mutations were found to be discordant between hospital and community samples. Conclusions regarding correlations between these mutations and COVID-19 disease severity may be subject to ascertainment bias. Since we do not have access to patient clinical outcomes or disease severity, we are unable to formally demonstrate that the mutations are linked to disease severity. Taken together, these findings merit further investigation of selective pressures other than immune evasion that may mediate SARS-CoV-2 evolution and differential disease severity in the community studied.
Baseline surveillance has immediate applications for public health preparedness. Notably, we found that the shift to Omicron in hospital admissions occurred 10.5 days after the rise in community cases. This means that sequencing surveillance can be used as one of the predictive indicators for emergency surge preparedness. This approach relies on consistent surveillance efforts, such as daily sampling resolution in this study. While Omicron infection is generally associated with less severe disease (45), a future variant with a shorter lag time between community to hospital may indicate more rapid progression to severe disease. Since functional characterization of new variants takes time, the rate of emergence can be benchmarked against Delta and Omicron. For example, an increased rate of emergence comparable to Omicron BA.1 may implicate antigenic shift, while a typical rate of emergence comparable to Delta/ BA.2/BA.5 may suggest antigenic drift in relation to population immunity. Further evolution of Omicron BA.1, BA.2, BA.4, and BA.5 descendant lineages (e.g., BF.7, BA.2.75.2, BQ.1.1, BF.13, BM.1.1.1, etc.) have separately converged on a set of spike gene mutations within the RBD (S: R346, S: L452, S: N460, S: F486), potentially indicating strong positive selection for neutralization resistance and Spike fusogenicity (46). Despite this convergent evolution, no currently circulating lineage has demonstrated the either complete displacement of competing lineages or high rate of emergence as demonstrated by BA.1 in our study. The current landscape of Omicron sublineages currently exhibits similarities to that of Delta sublineages prior to the emergence of Omicronseveral sublineages splintering from the parental lineage with a few demonstrating modest increases in fitness. These findings advance our understanding of the emergence dynamics of SARS-CoV-2 variants and demonstrate the importance of sequencing surveillance in public health preparedness for future variants.

MATERIALS AND METHODS
Study population. This study was approved by Arizona State University, St. Joseph's Hospital and Medical Center, Phoenix Children's Hospital, and Valleywise Health Medical Center Institutional Review Boards. For public health surveillance, all specimens were deidentified. Specimens were not identified in any targeted sampling effort. Community sequences analyzed in the study were from 27,071 saliva specimens submitted for testing at the ASU Biodesign Clinical Testing Laboratory between November 1st, 2021, and January 31 st , 2022. ASU's clinical testing laboratory conducted diagnostic COVID-19 testing that was free to the public. Saliva specimens were submitted from local collection drop-off sites and public testing sites across the state. Hospital inpatient specimens analyzed in the study include 1,125 positive diagnostic nasopharyngeal swab specimens from Phoenix Children's Hospital, St. Joseph's Hospital and Medical Center, and Valleywise Health Medical Center for SARS-CoV-2 genomic sequencing surveillance. RNA was extracted from 250 mL of saliva specimen within 33 h of sample receipt using the KingFisher Flex (Thermo Scientific, Waltham, MA, USA), following the manufacturer's guidelines. Diagnostic testing was performed using TaqPath COVID-19 Combo kit assay (Applied Biosystems, Waltham, MA, USA), following the manufacturer's guidelines.
SARS-CoV-2 sequencing. NGS library preparation for clinical specimens and wastewater samples was performed using the COVIDSeq Test (Illumina, San Diego, CA, USA) with ARTICv4 and ARTICv4.1 primer sets (47). Libraries were sequenced on the Illumina NextSeq2000 instrument using 2 Â 109 paired end reads. Sequencing reads adapter sequences were trimmed using trim-galore, aligned to the Wuhan1 reference genome (MN908947.3) using the Burrows-Wheeler aligner, BWA-MEM version 0.7.17-r1188 (48), and had their primer sequences trimmed using iVAR version 1.3.1 (49). Lineage calling for community and hospital-derived sequencing data were performed with pangolin software (50), with its assignment and designation libraries up to date at the time of analysis. Sequence quality was validated and annotated using VADR version 1.4 (51). Relative SARS-CoV-2 lineage and variant abundance analysis with Freyja (32) was performed only for wastewater samples that achieved $20 depth for $70% of positions across the reference genome.
Variant nonlinear regression analysis. Nonlinear regression of variant emergence was performed in GraphPad Prism Version 9.4.0 for Mac (GraphPad Software, San Diego, CA, USA), using the Sigmoidal, 4PL, X is log(concentration) model with no special handling of outliers, least-squares regression fitting method, and no weighting. For the Arizona rate of emergence experiment, SARS-CoV-2 GISAID submissions were downloaded using text and collection parameters. Text: "Arizona State University." Collection: 2021-01-12 to 2022-07-27. A total of 50,821 Arizona State University GISAID submissions were downloaded on 2022-07-27 matching these criteria, and 2,578 submissions corresponding to hospital specimens were removed from this set using internal specimen IDs. A second GISAID submission download to capture the latest BA.5 submissions was performed on 2022-08-30, using the text and collection parameters. Text: "Arizona State University." Collection: 2021-05-01 to 2022-08-30. 5,204 sequences Arizona State University GISAID submissions were downloaded on 2022-08-30 matching these criteria from which 607 hospital specimens were removed. For the rate of emergence control experiment, SARS-CoV-2 GISAID submissions were downloaded using the location and collection parameters. Location: Africa/South Africa. Collection: 2021-01-12 to 2022-07-27. 36,802 South Africa GISAID submissions were downloaded on 2022-08-07 matching these criteria, and 16 submissions were removed from this set due to incomplete collection date. A second GISAID submission download to capture the latest BA.5 submissions was performed on 2022-08-30, using the location and collection parameters. Location: Africa/South Africa. Collection: 2022-02-18 to 2022-08-30. 8,694 South Africa GISAID submissions were downloaded matching these criteria 5-day sliding averages for each variant (Delta, BA.1, BA.2, BA.5) were calculated for the range of collection dates and were utilized for regression modeling.
Mutation frequency comparisons. Variant mutation frequencies for each genome sequence were tabulated from GISAID annotation metadata. Comparisons between community and hospital abundance were performed using Fisher's exact test for mutations found in 5% or more of samples. In order to confirm sequence ambiguity was not causing mutations to not be annotated, mutations found to be significantly different between cohorts (a = 0.05) had their sequences confirmed for nonambiguous nucleotides at the mutation locus. Sites with a high frequency of ambiguous nucleotide calls were omitted from the analysis due to insufficient sequences (i.e., fewer than 10% of sequences remained after quality filter). Nucleotide frequencies at SNP loci were counted on wastewater bam files using Geneious Prime build 2022-07-07. Mutation abundances were adjusted to remove samples containing ambiguities. Comparisons between adjusted community and hospital abundances were performed using Fisher's exact test with Pvalues corrected using the Benjamini/Hochberg algorithm for multiple testing correction.
Network analysis. To generate viral transmission networks, genome sequences were first aligned using MAFFT software version 7.490 with 26merpair and -addfragments arguments, while specifying the Wuhan1 genome (MN908947.3) as the reference genome. A maximum likelihood distance matrix was created using IQTREE2 version 2.2.0.3, and by specifying a transition substitution model, empirically determined base frequencies, and the FreeRate rate heterogeneity model (-m TIM1F1R argument). The distance matrix was converted to an edge list, and source/target pairs exceeding a genetic distance of 0.00006 were removed. For each node, the most likely source node was determined by first finding the smallest genetic distance between all connected nodes and removing edges with distances exceeding that value. Next, the difference in collection dates between all connected nodes for each target node was calculated. Edges for nodes with collection dates greater than the minimum difference were removed, as well as nodes with collection dates greater than 10 days apart. Next, the Haversine distances between zip codes for each source and target node was calculated and all edges greater than the minimum calculated distance were removed.
The network layout was arranged in Cytoscape (V 3.9.1) using a prefuse force directed layout with genetic distance as edge weights.
Lineage specific mutations. Lineage specific mutations for each variant were compiled from the five most abundant sublineages of each variant found in community samples over the study period (Table S1). Mutation prevalence from each sublineage was obtained from Outbreak.info (52,53). To be considered a lineage specific mutation, mutations were required to have prevalence greater than 95% in one variant and less than 1% in the other. Mutation frequency for each wastewater sample was calculated from variant files created using iVAR (V 1.3.1) with a minimum base quality filter of 20 and the Wuhan1 genome (MN908947.3) as a reference.
Wastewater sample collection, processing, and analysis. Three hundred seventy composite samples of untreated wastewater were collected over 24 h from within the wastewater collection system using high frequency automated samplers at 14 locations in the Greater Tempe area, Arizona (population of approximately 700,000), between November 2, 2021, and January 29, 2022. Samples were transferred to 1L high-density polyethylene (HDPE) bottles and stored at 4°C until nucleic acid extraction. Wastewater samples were processed as previously described (31). Briefly, raw wastewater was filtered through a 0.45 mm polyethersulfone (PES) membrane filter (Fisher Scientific, Lenexa, KS, USA). Viruses present in the resultant solids-depleted filtrate were concentrated on an Amicon ultra 15 centrifugal filter with a 10,000 molecular weight cutoff (MWCO) (Millipore Sigma, Burlington, MA, USA). RNA was extracted from the concentrate using a Qiagen RNeasy minikit (Qiagen, Germantown, MD, USA). Reverse transcription-quantitative PCR (RT-qPCR) was performed on an Applied Biosystems QuantStudio 3 realtime PCR system using SuperScriptIII One-Step RT-PCR System with Platinum Taq DNA polymerase (Invitrogen, Carlsbad, CA, USA) and the Charité/Berlin (World Health Organization) protocol primer and probe E (envelope) gene target.
Omicron RT-dPCR assay. The Omicron RT-dPCR assays (Omn143) was performed using the QuantStudio Absolute Q Digital PCR System (Applied Biosystems, Waltham, MA, USA) in 10 mL reaction mixtures (including 10% overage) containing 2.5 mL of 4Â Combinati one-step RT-dPCR MasterMix, 0.4 mL of each 10 mM primer (Omn143-Fwd and Omn143-Rev; Table S2), 0.2 mL of 10 mM probe (Omn143 Probe; Table S2), and 6.5 mL of sample RNA. Reaction mixtures were loaded into QuantStudio Absolute Q MAP16 plates and overlaid with 15 mL of isolation buffer. The assay was performed with thermal cycling conditions: reverse transcription at 50°C for 10 min, activation/denaturation at 95°C for 5 min, followed by 45 cycles of denaturation at 95°C for 5 s and annealing/extension at 54°C for 15 s. Assay validation was performed in three independent experiments using synthesized gBlock gene fragments (Integrated DNA Technologies, Coralville, IA, USA): Omicron BA.1 control and Delta control (Table S2). Microreaction chambers that passed QC, as determined by ROX signal, were analyzed. FAM fluorescence intensity threshold to call Omicron positive detection per microchamber was set at 10,000 (more than two standard deviations above water and Delta negative controls). Absolute quantification (copies) was reported.
Data availability. Code used for network analysis is available at https://github.com/ASU-Lim-Lab/ Network_analysis. Genome sequences have been deposited to the GISAID repository.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only.