Retrospective Surveillance of Wastewater To Examine Seasonal Dynamics of Enterovirus Infections.

Enteroviruses are RNA viruses that are responsible for both mild gastroenteritis and mild respiratory illnesses as well as debilitating diseases such as meningitis and myocarditis. The disease burden of enteroviruses in the United States is difficult to assess because most infections are not recorded. Since infected individuals shed enterovirus in feces and urine, surveillance of municipal wastewater can reveal the diversity of enteroviruses circulating in human populations. Therefore, monthly municipal wastewater samples were collected for 1 year and enteroviruses were quantified by reverse transcriptase quantitative PCR and identified by next-generation, high-throughput sequencing. Enterovirus concentrations ranged from 3.8 to 5.9 log10 equivalent copies/liter in monthly samples. From the mean monthly concentration, it can be estimated that 2.8% of the contributing population was shedding enterovirus daily. Sequence analysis showed that Enterovirus A and Enterovirus B alternate in predominance, with Enterovirus B comprising over 80% of the reads during the summer and fall months and Enterovirus A accounting for >45% of the reads in spring. Enterovirus C was observed throughout the year, while Enterovirus D was present intermittently. Principal-component analysis further supported the date corresponding to enterovirus seasonal trends as CVA6 (Enterovirus A) was predominant in the spring months; CVB3, CVB5, and E9 (Enterovirus B) were predominant in the summer and fall months; and CVA1, CVA19, and CVA22 (Enterovirus C) and EV97 (Enterovirus B) were predominant in winter. Rhinoviruses were also observed. Wastewater monitoring of human enterovirus provided improved insight into the seasonal patterns of enteroviruses circulating in communities and can contribute to understanding of enterovirus disease burden. IMPORTANCE Enterovirus infections are often not tracked or reported to health officials. This makes it hard to know how many people in a community are infected with these viruses at any given time. Here, we explored enterovirus in municipal wastewater to look at this issue. We show that enteroviruses are present year-round in municipal wastewater at levels of up to 800,000 genomic copies per liter. We estimate that, on average, 2.8% of the people contributing to the wastewater shed enterovirus daily. Sequence analysis of the viral capsid protein 4 gene shows that 8 enterovirus types are key drivers of seasonal trends. Populations of Enterovirus A members peak in the spring, while Enterovirus B types are most prevalent during the summer and fall months and Enterovirus C members influence the winter months. Enterovirus D was observed sporadically and did not influence seasonal trends.

INTRODUCTION

Enteroviruses (EV) are enteric viruses that predominantly replicate in the gut of infected animals, including humans. During human infection, virus particles are shed into the sewerage system for a period of time lasting as long as 16 weeks (1) in concentrations as high as 106 particles per gram of feces (2). The EV consist of about 116 distinct members, including serotypes of echoviruses (e.g., echovirus 30), polioviruses, coxsackievirus A and B, and the numbered enteroviruses (e.g., EV71) (3). They are categorized into four species, currently named Enterovirus A (EV-A), Enterovirus B (EV-B), Enterovirus C (EV-C), and Enterovirus D (EV-D) in the genus Enterovirus, family Picornavirdae (4). Nine additional species are found in the genus Enterovirus, including Enterovirus E, Enterovirus F, Enterovirus G, Enterovirus H, Enterovirus I, and Enterovirus J, which primarily infect animals, as well as Rhinovirus A, Rhinovirus B, and Rhinovirus C, which includes viruses that cause respiratory disease in humans.

EV are transmitted through the fecal-oral route, and clinical symptoms range from mild infections that include acute gastrointestinal and respiratory symptoms to severe, debilitating conditions, such as encephalitis and paralysis, as reported for the 2014 outbreak of EV-D68 in the United States (5, 6). EV have also been associated with chronic diseases, such as type 1 diabetes (7). It has been estimated that EV cause 10 to 15 million infections and tens of thousands of hospitalizations annually in the United States (8). In 2012, two EV (poliovirus and enterovirus 71) were listed in the five top global infectious disease threats being monitored by the Global Disease Detection Operations Center of the Centers for Disease Control and Prevention (CDC) (9). Among EV infections, only paralytic polio is a nationally reportable disease (10), which makes it difficult to assess the actual disease burden attributable to EV. However, some nationwide data on EV infections has been collected through the National Enterovirus Surveillance System (NESS) of the CDC (11). Through this passive surveillance system, diagnostic laboratories voluntarily report diagnoses of EV infections. NESS data collected from 1970 to 2013 reveal a consistent peak in EV infections from June through October. Additionally, serotypes of EV-B were reported to NESS most frequently, followed by EV-A and then EV-C and, finally, EV-D (10, 12, 13). Furthermore, EV-A, EV-B, and EV-C have been associated with a variety of clinical sample types, including cerebrospinal fluid and throat swabs (14–17).

Municipal wastewater (sewage) surveillance of EV can uncover the enteroviruses serotypes circulating in a community and serve as a surrogate for clinical data of EV infections. Additionally, identifying the sources of genetic diversity (i.e., mutation and recombination) of enteroviruses can provide relevant information regarding the emergence of pathogenic forms of EV. Detection of EV serotypes in sewage has been reported (18–33). However, the results of those studies show few observations and a narrow range of EV serotypes present, possibly due to the culture-dependent approaches used to enrich for EV (18–30) and/or to the use of traditional sequencing techniques that have limited capacity and result in detection of only abundant taxa. Furthermore, the recovery efficiency of the methods used in those studies was not reported, and the methods used may have favored recovery of specific EV species from sewage. Nevertheless, these studies have provided valuable data, including the fact that EV concentrations are typically highest in the summer-fall time frame, which aligns with the NESS reports of peak clinical infection. Additionally, most studies have reported patterns of prominence in municipal wastewater similar to those reported in the clinical NESS data, with members of EV-B identified most frequently, followed by EV-A and then EV-C; members of EV-D are rarely reported in municipal wastewater. Due to the limited range of EV species observations in each wastewater surveillance study and the dominance of EV-B detections, the overall EV species distribution has not been well described on a seasonal basis.

The ability to accurately identify enterovirus in municipal wastewater is dependent on the effectiveness and limitations of the methods used. Historically, typing of EV from clinical and sewage samples was performed by first propagating virus on various cell lines and then identifying EV with virus-specific antisera (34–37). These methods can be cost prohibitive and labor-intensive since there is not a single cell line that supports propagation of all enterovirus serotypes. Furthermore, the antisera used for EV identification are no longer commercially available and are not likely existent for the newer identified types. More-recent typing has relied on gene sequencing of the VP1 (36) or VP4 (37) viral capsid protein. Culture-independent strategies such as next-generation, high-throughput sequencing facilitate the interrogation of a particular sample with the depth necessary to observe the presence of both abundant and minority members (17, 38–40). Next-generation sequencing (NGS) technologies have allowed the identification of viral pathogen groups present in environmental samples, such as biosolids (38) and municipal wastewater (41), despite the high background of nonpathogenic viruses. NGS has also increased the discovery rate of human-pathogenic viruses in feces (39) and has revealed that culture-based enrichment can indeed alter the taxonomic profiles of environmental samples (40). Recently, NGS approaches have been described that examine the diversity of adenovirus (42), norovirus (43), and astrovirus (43) in sewage. Those studies highlighted the utility NGS approaches for community-wide enteric virus surveillance and offered advantages over Sanger sequencing, which provides information only on the dominant circulating strain (42), and can facilitate identification of virulent recombinant virus strains that may negatively impact public health (43).

The objective of this study was to quantify the load of EV in municipal wastewater and to use NGS to determine the taxonomic distribution and seasonal dynamics of EV using two gene targets, VP1 and VP4.

RESULTS

qPCR of EV in monthly municipal wastewater samples.

Reverse transcriptase quantitative PCR (RT-qPCR) and linear regression analysis were used to determine the monthly quantities of total EV collected from two wastewater treatment plants (WWTPs) over the course of 1 year. The mean EV concentrations for March (2010) to February (2011) are shown in Fig. 1. The lowest concentration of EV, 7.05 × 103 equivalent copies/liter, was observed in February, and the highest concentration, 8.3 × 105 equivalent copies/liter, was measured in July. Friedman repeated-measure analysis of variance (ANOVA) on ranks showed borderline significance by month (P value = 0.046), but the use of Tukey's multiple-comparison procedure could not resolve any significant differences among the monthly ranked sums. When the data were categorized by season (summer, June, July, and August; fall, September, October, and November; winter, December, January, and February; spring, March, April, and May), a one-way repeated-measure ANOVA showed that the seasons were significantly different (P value < 0.001), with summer months exhibiting significantly higher EV than any other season (Holm-Sidak multiple-comparison procedure; P value < 0.05).

FIG 1 

Mean quantities of enterovirus detected in monthly primary effluent from two WWTPs by RT-qPCR. Quantities are reported in log10 equivalent copy numbers per liter, and error bars represent standard deviations. Mar, March; Apr, April; Jun, June; Jul, July; Aug, August; Sep, September; Oct, October; Nov, November; Dec, December; Jan, January; Feb, February.

Taxonomic analysis of VP1 and VP4 sequence reads.

The results from processing and taxonomic classification of VP1 and VP4 amplicon reads are presented in Tables 1 and 2, respectively. The number of raw reads obtained from VP1 PCR products ranged from 5,239 to 13,761, and after processing, 16% to 81% of the raw reads were retained and determined to be unique. Various proportions (1% to 97% of unique reads by month) of these reads were classified by BLAST analysis using The NCBI Virus Reference Database, and 29% to 100% of these reads were assigned to EV. The proportionate distributions of VP1 read assignments to EV-A, EV-B, EV-C, and EV-D by month are shown in Fig. S1 in the supplemental material. Other assignments were made to Rhinovirus A, Rhinovirus B, and simian picornavirus (see Table S1 in the supplemental material). The number of raw reads obtained from VP4 PCR products was slightly higher than the number obtained from VP1, ranging from 7,687 to 24,025. Processing of the reads resulted in 26% to 60% retention, and almost all (≥96%) reads were classified at the species level; 19% to 95% were EV-A, EV-B, EV-C, and EV-D. Other taxonomic assignments were Rhinovirus A, Rhinovirus B, Rhinovirus C, canine picornavirus, and tomato mosaic virus (Table S1).

TABLE 1 

Analysis of VP1 reads obtained from sequencing

Moc	No. of raw reads	Filter processa		Dereplication processb			% of unique reads classifiede	% of unique classified as EV A– EV D
		No. of reads	% retained	No. of unique reads	No. of reads in largest clusterd	% of unique reads
April	9,719	2,977	31	2,515	13	26	1	100
May	6,538	1,077	17	1,011	4	16	33	99
June	10,340	4,654	45	3,991	15	39	53	100
July	10,387	4,188	40	3,895	6	38	78	69
August	5,239	4,725	90	4,266	9	81	68	98
September	7,052	4,565	65	4,015	9	57	97	87
October	13,761	3,321	24	2,912	10	21	14	29
November	12,909	8,080	63	6,329	24	49	76	100

Reads were filtered as follows: minimum read length, 300 nt; average Q score, ≥30; <8 runs of homopolymers.

Dereplication was performed to find unique reads by clustering at 100% sequence identity.

Samples collected in March, December, January, and February were not analyzed.

Data represent the number of reads in the largest cluster after dereplication.

Classification was performed by BLAST analysis of the NCBI Virus Reference Database.

TABLE 2 

Analysis of VP4 reads obtained from sequencing

Mo	No. of raw reads	Filter processa		Dereplication processb			% of unique reads classifiedd	% of unique classified as EV A–EV D
		No. of reads	% retained	No. of unique reads	No. of reads in largest clusterc	% of unique reads
March	11,117	5,952	54	4,578	374	41	96	19
April	22,099	15,110	68	11,730	461	53	99	95
May	16,281	7,303	45	6,395	104	39	99	76
June	16,452	9,576	58	8,640	62	53	99	80
July	18,380	11,340	62	10,028	116	55	100	88
August	7,687	3,190	41	3,132	5	41	99	90
September	24,025	17,516	73	14,314	160	60	100	66
October	19,468	12,306	63	10,585	208	54	100	30
November	13,968	7,400	53	6,611	115	47	99	30
December	11,761	5,832	50	5,209	99	44	99	46
January	10,766	5,094	47	4,563	38	42	99	82
February	11,090	2,284	26	2,389	47	22	98	53

Reads were filtered as follows: minimum read length, 300 nt; reads must contain the forward primer sequence; average Q score, ≥30; <8 runs of homopolymers.

Dereplication was performed to find unique reads by clustering at 100% sequence identity.

Data represent the number of reads in the largest cluster after dereplication.

Classification was performed by BLAST analysis of the NCBI Virus Reference Database.

Monthly (April to November) distribution of EV-A, EV-B, EV-C, and EV-D based on Roche 454 GS FLX sequencing of VP1. Download FIG S1, PDF file, 0.1 MB.

Read assignments to other viruses. Download TABLE S1, PDF file, 0.1 MB.

Correlation of EV sequence abundances to RT-qPCR quantities.

The mean quantities of EV measured in monthly WWTP samples by RT-qPCR were compared to the number of unique reads assigned to the EV species from sequence analysis of VP1 and VP4. The RT-qPCR quantities had a positive correlation with the number of VP1 reads (Pearson product-moment correlation coefficient = 0.574) and the number of VP4 reads (Pearson product-moment correlation coefficient = 0.666). However, the correlation between numbers of VP1 reads and RT-qPCR quantities was not significant (P value = 0.137), while significance was observed for VP4 reads (P value = 0.018). Additionally, there was not a significant relationship among the numbers of VP1 and VP4 reads in monthly samples, as the Pearson product-moment correlation coefficient was determined to be −0.264 (P value = 0.527). Thus, VP4 was selected for subsequent analysis because of its taxonomic precision and significant correlation with EV quantities obtained by RT-qPCR.

Monthly distribution of EV species based on VP4.

The proportionate distributions of VP4 read assignments to EV-A, EV-B, EV-C, and EV-D by month are shown in Fig. 2. EV-A, EV-B, EV-C, and EV-D were detected in 92%, 100%, 100%, and 42% of the months examined, respectively. EV-D was rarely detected from May through November, and EV-C was present at relatively low levels in all of the months. EV-A and EV-B tended to alternate in predominance. EV-A was predominant from March through May, with a significant decreasing trend (Mann-Kendall, P value = 0.0075). EV-B abundance was greatest from June through February and exhibited an increasing trend during May and a decreasing trend during early winter but was not unidirectional (Mann-Kendall, P value = 0.30). Overall, these two distributions were negatively correlated (r = −0.74, P value = 0.0058), with maximum correlation during March (r = 0.97, P value = 0.00093) and minimum correlation during September (r = −0.845, P value = 0.00054). A borderline significant upward trend was observed for EV-C abundance (Mann-Kendall, P value = 0.064) from November to February, and the dearth of data for EV-D (7 months with 0% relative abundance) resulted in a nonsignificant trend (Mann-Kendall, P value = 1.0).

FIG 2 

Monthly distribution of EV-A, EV-B, EV-C, and EV-D based on classification after Roche 454 GS FLX sequencing of the VP4 gene.

Monthly distribution of EV serotypes.

Sequences matching 85 EV serotypes were detected in the monthly sewage samples, and the numbers of reads assigned to each serotype are shown according to species in Tables S2, S3, S4, and S5. Relative abundances of the 85 sewage serotypes (plus 3 untyped bins of EV-B, EV-C, and EV-D serotypes) are illustrated as a matrix plot in Fig. 3A. The overall predominant trends were EV-A serotypes in the spring; EV-B serotypes in the summer, fall, and early winter; and EV-C serotypes in the winter. Principal-component analysis (PCA; see section below) delineated 8 serotypes (i.e., CVA6, CVB3, CVB5, E9, EV97, CVA19, CVA22, and CVA1) that played a significant role in the prevalence of EV. The relative abundances of these 8 serotypes are shown in Fig. 3B. The spring months (March, April, and May) were noticeably predominated by EV-A serotype CVA6, representing 86.3%, 57.6%, and 35.6% of the reads during these months and 18.1%, 4.5%, and 0.6% in the subsequent months of June, July, and August. There was a significant decreasing trend for CVA6 (Mann-Kendall, P value = 0.002). In contrast, EV-B serotypes CVB3, CVB5, and E9 levels predominated in the summer and fall months of June through December, with peak levels in June (E9), July (CVB3), and November (CVB3 and CVB5). The abundance of serotype CVB3 increased from May to July and decreased from November to February, and it had no significant unidirectional trend (Mann-Kendall, P value = 0.63). A steady increase followed by a steep decline was observed for the abundance of CVB5 from June to November, and it exhibited a significant trend (Mann-Kendall, P value = 0.03). E9 was observed in April, and its abundance peaked in June and followed a steady decline through November, with no significant unidirectional trend (Mann-Kendall, P value = 0.22). EV-C serotypes CVA19, CVA22, and CVA1 predominated in December, January, and February. In addition, EV-B serotype EV97 was also observed during the winter. The Mann-Kendall trends for the abundances of these 4 serotypes were not significant (CVA19, P value = 0.077; CVA22, P value = 0.072; CVA1, P value = 0.27; EV97, P value = 0.29).

FIG 3 

(A) Matrix plot of all EV serotypes detected (85 serotypes and 3 untyped bins of EV-B, EV-C, and EV-D). (B) Matrix plot of 8 impactful serotypes classified in the 12 monthly municipal wastewater samples. The EV serotype abbreviations are as follows: CVA, coxsackievirus A; CVB, coxsackievirus B; EV, enterovirus; E, echovirus; SP, simian picornavirus; SVDV, swine vesicular disease virus. The monthly samples are indicated by the 3-letter abbreviation. The scale represents percentage of reads per wastewater sample.

Number of reads assigned to EV-A serotypes in monthly municipal wastewater samples. Download TABLE S2, PDF file, 0.1 MB.

Number of reads assigned to EV-B serotypes in monthly municipal wastewater samples. Download TABLE S3, PDF file, 0.1 MB.

Number of reads assigned to EV-C serotypes in monthly municipal wastewater samples. Download TABLE S4, PDF file, 0.1 MB.

Number of reads assigned to EV-D serotypes in monthly municipal wastewater samples. Download TABLE S5, PDF file, 0.1 MB.

Detection of rhinovirus species in municipal wastewater by VP4 sequencing.

Although not a focus of this study, rhinoviruses were also detected in all monthly samples. The relative abundances of rhinovirus species RV-A, RV-B, and RV-C are shown in Fig. S2. RV-A and RV-C exhibited high levels in spring and fall, respectively, but were nearly absent in the summer. There was a significant decreasing trend for RV-A (Mann-Kendall, P value = 0.019) and a significant increasing trend for RV-C (Mann-Kendall, P value = 0.042). In contrast, RV-B was present throughout the year and exhibited no significant unidirectional trend (Mann-Kendall, P value = 0.837).

Monthly distribution of rhinoviruses based on Roche 454 GS FLX sequencing of VP4. Download FIG S2, PDF file, 0.1 MB.

Principal-component analysis of EV species and serotypes.

Principal-component analysis was performed on the VP4 reads to examine the associations of the abundance of each EV species and the month of detection (Fig. 4A). The first and second principal components (i.e., PC1 and PC2) contained 80.0% and 19.9% of the variance and together comprised 99.9%. The spring, summer, and fall months were separated along PC1, whereas the winter months were separated from those seasons on PC2. The influential species loadings on PC1 were EV-A (correlation to PC1 = −0.95) and EV-B (correlation = 0.97), while EV-C (correlation = −0.14) and EV-D (correlation = −0.08) were less influential. EV-C was significantly correlated with PC2 (correlation to PC2 = 0.99) followed by EV-A (correlation = −0.32), EV-B (correlation = −0.24), and EV-D (correlation = −0.04). Vectors of these loading coefficients (purple and green arrows, Fig. 4A) indicated that EV-A and EV-B together had the greatest influence on the dynamics of enterovirus infections for the nonwinter months but that they were negatively correlated to each other (Spearman −0.74, P value = 0.006). EV-C had the greatest influence on EV infections during the winter months.

FIG 4 

Principal-component analysis of the EV species (A) and serotype distributions (B) in monthly samples. The monthly samples are shown as abbreviations and grouped according to season (blue, winter; pink, spring; red, summer; orange, fall). The biplot arrows were proportionally rescaled according to the right triangle ratio tan(α) = O/A, where O and A are opposite and adjacent sides of angle α of the arrow relative to the center of the coordinate system. (A) The influence of each of the 4 EV species is indicated by colors (purple, EV-A; green, EV-B; royal blue, EV-C; black, EV-D). The biplot arrow for EV-D did not project in the coordinate system and is shown as a black box. (B) Serotypes with loading coefficients greater than 0.1 or less than −0.1 are indicated by colored arrows for each EV species (purple, EV-A; green, EV-B; royal blue, EV-C). The sign of the axis 2 component values was reversed to rotate the plot 180° for comparison with the data in panel A.

Principal-component analysis was also performed on the percentage of reads associated with 85 EV serotypes (and 3 untyped EV bins) for each of the 12 monthly municipal wastewater samples. PC1 and PC2 explained 64.2% and 16.3% of the variance, respectively, with a cumulative variance of 80.5% (Fig. 4B). The monthly samples clustered by season, with the variance of spring, summer, and fall maximized along PC1 and the variance of winter on PC2. Of the 85 EV serotypes observed, 8 serotypes (i.e., CVA6, CVB3, CVB5, E9, EV97, CVA19, CVA22, and CVA1) had loading coefficients higher than 0.1 or lower than −0.1 and are displayed in Fig. 4B (see Fig. S3 for a plot of the coefficients). The loading coefficients of the remaining 78 serotypes (and 3 untyped EV bins) were not clearly distinguishable on the PCA plot and are not shown. Biplot PCA revealed that these 8 serotypes had significant relationships with the seasons (Fig. 4B). Serotype CVA6 (purple arrow) of EV-A was negatively correlated with PC1 (correlation = −0.99) and was associated with the spring season. In contrast, EV-B serotypes CVB3 and CVB5 (green arrows; PC1 correlations of 0.70 and 0.57, respectively) were positively correlated with this axis and were associated with the late summer, fall, and early winter months. EV-B serotype E9 (green arrow pointing downward) was not significantly correlated with PC1 (correlation = 0.02) but was negatively correlated with PC2 (correlation of −0.48) and was associated with the early summer months. Finally, 4 of the 8 serotypes were significantly correlated with PC2 and associated with the middle to late winter months. These included EV-C serotypes CVA1 (PC2 correlation = 0.63), CVA19 (PC2 correlation = 0.93), and CVA22 (PC2 correlation = 0.94) and EV-B serotype EV97 (PC2 correlation = 0.88). Taken together, these analyses showed that CVA6 had the greatest influence on EV infections occurring during spring whereas 3 of the 4 EV-B serotypes were driving EV infections during summer, fall, and early winter. Finally, 3 EV-C serotypes as well as EV-B EV97 contributed to the EV infections occurring during January and February.

Bar chart of the EV serotypes with PC1 and PC2 loading coefficients greater than 0.1 and less than −0.1. Black bars are PC1 coefficients, and gray bars are PC2 coefficients. Download FIG S3, PDF file, 0.1 MB.

Visual inspection of Fig. 4 suggested a high level of concordance between species and serotype PCA plots. To confirm this observation, the similarity between the PCA plots was determined by Procrustes analysis (PA) and the Mantel r test. The objective of PA is to match the structure of one PCA plot to another by rotating, translating, and dilating their configurations, and the Mantel r test determines the correlations between the positions of 2 Euclidean distance matrices. Initially, these statistical tests were used to compare a data set containing all 85 EV serotypes (and 3 untyped EV bins) to one containing 8 EV serotypes with loading coefficients higher than 0.1 or lower than −0.1. The Procrustes rotation correlation was 0.9997 (P value = 0.001), where a correlation of 1.0 would indicate identical PCA structures. The minimum, median, and maximum residual errors were 0.0016 (April) and 0.0045 and 0.016 (January), and the Mantel r statistic was 0.98 (Spearman’s rank correlation, P value = 0.001; number of permutations, 999). Thus, these plots had nearly identical Euclidean matrices and PCA configurations as shown in Fig. S4, top. Next, PA and the Mantel r test were performed on the species PCA plot (i.e., Fig. 4A) and a PCA plot was generated from the 8 serotypes with influential loading coefficients. The minimum, median, and maximum residual errors were 0.022 (August) and 0.073 and 0.227 (January). The data from the summer (June, July, and August) and fall (September, October, and November) months showed greater concordance, whereas the data from some of the winter (January) and spring (March and May) months were slightly discordant (Fig. S4, bottom). Overall, the species and serotype PCA plots were very similar to each other, with a Procrustes rotation correlation of 0.91 (P = 0.001) and a Mantel r value of 0.81 (Spearman’s rank correlation, P value = 0.001; number of permutations, 999).

Procrustes analysis of superimposed PCA plots. (Top) PCA plots of the 85 serotypes and the 8 most significant serotypes. (Bottom) PCA plots of the species and the 8 most significant serotypes. The month labels indicate the locations of the serotype PCA scores, and the blue arrows indicate the vector distance between these scores (base of arrow) and the species PCA scores (arrowhead) due to rotation and scaling. No arrows are shown in panel A because the vector errors were too small. The solid and hashed cross-hair lines centered at 0/0 represent PCA plots for the 8 most significant serotypes and 85 serotypes (top) and the species and serotypes (bottom). Download FIG S4, PDF file, 0.1 MB.

Comparison of municipal wastewater EV serotypes to those in the NESS database.

The rank abundances of the wastewater EV serotypes were compared to serotype data reported to NESS in 2010 and 2011 (12) (Table 3). The top 5 wastewater serotypes, which comprised 60% of the reads in all the monthly samples, were CVB3 (19.1%), CVA6 (17.6%), E9 (10.1%), CVB5 (7.2%), and E7 (5.3%). The top 5 serotypes reported to NESS in 2010 were E6 (17.5%), E18 (15%), CVB5 (11.5%), E9 (11.5%), and E7 (10%), while E6 (15.9%), CVB3 (15.2%), E30 (10.8%), CVB1 (8.2%), and E18 (8.2%) were the top 5 reported in 2011. The total abundances for the top 5 serotypes reported to NESS in 2010 and 2011 were 65% and 58%, respectively.

TABLE 3 

Rank abundance of EV serotypes identified in municipal wastewater samples and reported to NESS in 2010 to 2011

Serotype	Species (EV)	Wastewater (12 mo)		NESS 2010 (n = 199 cases)		NESS 2011 (n = 232 cases)
		% reads	Rank	% cases	Rank	% cases	Rank
CVB3	B	19.1	1	4	7	15.5	2
CVA6	A	17.6	2	2	13	3	11
E9	B	10.1	3	11.5	4	3.8	10
CVB5	B	7.2	4	11.5	3	6.5	6
E7	B	5.3	5	10	5	1.3	17
CVB1	B	4.1	6	—a	—	8.2	4
E25	B	3.4	7	—	—	—	—
Untyped EV-B	B	3	8	—	—	—	—
CVB4	B	2.5	9	3	11	1.3	12
E11	B	2.5	10	2	15	—	—
E1	B	1.9	11	—	—	—	—
CVA1	C	1.72	12	—	—	—	—
E14	B	1.67	13	—	—	—	—
CVB2	B	1.53	14	2	14	1.3	16
CVA22	C	1.45	15	—	—	—	—
CVA10	A	1.29	16	3	10	—	—
CVA19	C	1.26	17	—	—	—	—
CVA21	C	1.01	18	—	—	—	—
EV107	B	1.00	19	—	—	—	—
E30	B	0.94	20	2.5	12	10.8	3
CVA16	A	0.93	21	3.5	9	3.4	9
E6	B	0.819	22	17.5	1	15.9	1
CVA13	A	0.817	23	—	—	—	—
EV97	B	0.79	24	—	—	—	—
Untyped EV-C	C	0.77	25	—	—	—	—
CVA7	A	0.76	26	—	—	—	—
CVA9	B	0.73	27	6.5	6	6	8
CVA5	A	0.70	28	—	—	—	—
CVA20	C	0.59	29	—	—	—	—
E16	B	0.52	30	—	—	—	—
E3	B	0.48	31	—	—	—	—
CVA2	A	0.32	32	—	—	—	—
E21	B	0.31	33	—	—	—	—
EV78	B	0.28	34	—	—	—	—
SVDV-2	B	0.26	35	—	—	—	—
EV98	B	0.26	36	—	—	—	—
E15	B	0.21	37	—	—	—	—
CVA17	C	0.2	38	—	—	—	—
E20	B	0.19	39	—	—	—	—
EV71	A	0.17	40	—	—	1.3	15
E12	B	0.17	41	—	—	—	—
CVA24	C	0.14	42	—	—	—	—
CVB6	B	0.13	43	—	—	—	—
SP	NA	0.13	44	—	—	—	—
CVA8	A	0.12	45	—	—	—	—
E19	B	0.08	46	—	—	—	—
EV77	B	0.07	47	—	—	—	—
EV88	B	0.07	48	—	—	—	—
CVA3	A	0.05	49	—	—	—	—
EV109	C	0.04	50	—	—	—	—
PV2	C	0.03	51	—	—	—	—
EV99	C	0.023	52	—	—	—	—
PV1	C	0.021	53	—	—	—	—
EV68	D	0.019	54	2	16	6.5	7
EV75	B	0.016	55	—	—	—	—
EV73	B	0.015	56	—	—	—	—
EV79	B	0.015	57	—	—	—	—
E18	B	0.013	58	15	2	8.2	5
EV81	B	0.013	59	—	—	—	—
EV82	B	0.013	60	—	—	—	—
EV74	B	0.012	61	—	—	—	—
EV80	B	0.012	62	—	—	—	—
E5	B	0.01	63	—	—	—	—
E4	B	0.009	64	4	8	—	—
E13	B	0.009	65	—	—	—	—
EV104	C	0.008	66	—	—	—	—
EV70	D	0.008	67	—	—	—	—
Untyped EV-D	D	0.008	68	—	—	—	—
CVA14	A	0.007	69	—	—	—	—
EV87	B	0.007	70	—	—	—	—
EV100	B	0.006	71	—	—	—	—
SV43	A	0.005	72	—	—	—	—
E2	B	0.005	73	—	—	—	—
PV3	C	0.005	74	—	—	—	—
CVA12	A	0.004	75	—	—	—	—
CVA4	A	0.004	76	—	—	1.3	14
E26	B	0.004	77	—	—	—	—
SV46	A	0.002	78	—	—	—	—
E24	B	0.002	79	—	—	—	—
EV76	A	0.002	80	—	—	—	—
E29	B	0.002	81	—	—	—	—
E32	B	0.002	82	—	—	—	—
EV85	B	0.002	83	—	—	—	—
EV101	B	0.002	84	—	—	—	—
EV90	A	0.001	85	—	—	—	—
EV91	A	0.001	86	—	—	—	—
EV84	B	0.001	87	—	—	—	—
EV69	B	0.0006	88	—	—	—	—
E17	B	—	—	—	—	2.2	13

—, lack of detection in wastewater or reports by NESS.

Approximately 19% of the 85 wastewater serotypes were observed in the NESS database for 2010 and 2011. Of the 5 most abundant serotypes in wastewater, 4 EV-B serotypes (CVB3, E9, CVB5, and E7) were also ranked in the top 5 most frequently reported cases to NESS in either 2010 or 2011. Conversely, EV-A serotype CVA6 was ranked second in overall abundance in wastewater, and although it was reported to NESS in 2010 and 2011, it was not ranked in the top 5 reported EV serotypes (rank 13, 2010; rank 11, 2011). Serotypes E6, E18, and E30 were among the top 5 EV serotypes reported to NESS in 2010 and 2011, but each of these serotypes accounted for less than 1% of the total EV VP4 reads observed in the monthly wastewater samples.

DISCUSSION

In this study, EV was concentrated and extracted from primary effluent of municipal wastewater using a Celite (diatomaceous earth) method shown previously to recover EV serotypes representing the individual species with statistically similar rates (44). In addition, the wastewater samples were processed without cell culture enrichment, which also has been suggested to introduce a bias for the selection of a subset of EV (45). Furthermore, the use of a next-generation sequencing platform facilitated species identification of less-abundant members of the community since the number of sequence reads available is orders of magnitude greater than the practical limits of traditional cloning and sequencing efforts.

Two gene targets, VP1 and VP4, were evaluated for their ability to represent the seasonal distribution of EV species in municipal wastewater. To do this, published RT-PCR assays were implemented for VP1 (46) and VP4 (37) that have been shown to work with both prototype and clinical strains of EV. These two gene targets represent a hypervariable region (VP1) and a more conserved region (VP4) of the viral capsid. Although the use of VP1 RT-PCR for molecular typing of EV isolates from clinical samples has been established, it was found to be inefficient for targeting EV in wastewater samples, resulting in fewer classified reads than are obtained by targeting the VP4 region. This was likely due to the cross-reactivity of the degenerate VP1 primers (432 and 648 possible outcomes with primer sequences 292 and 222, respectively) with unknown bacteriophages present in the complex background of wastewater and the finite number of amplicons that can be sequenced in a single run. In fact, sequencing VP1 amplicons without prior DNase treatment of the extracts resulted in few reads classified as EV and many classified as phage (data not shown). Therefore, it was not surprising that the VP1 reads assigned to EV did not significantly correlate with the EV quantities determined for each sample by RT-qPCR. In contrast, VP4 reads significantly correlated with EV quantities, with >96% of processed monthly reads classified with the Virus Reference Database. These data demonstrated that the VP4 primers were more accurate than the VP1 primers for determining EV species distributions in complex environmental samples such as wastewater.

Unexpectedly, more than half of the VP4 reads for March, October, November, and December were classified as rhinovirus (Table S1) and, as expected, rhinovirus exhibited a characteristic spring and fall seasonal trend (see Fig. S2 in the supplemental material) (47). Although the detection of rhinoviruses in municipal wastewater is intriguing, it was not surprising that the VP4 primers used in this study cross-reacted with rhinoviruses because the MD91 forward primer matched perfectly to many NCBI rhinovirus genotypes (data not shown) and the OL-68 reverse primer has been shown to anneal to rhinovirus genomes (48). Determination of rhinovirus occurrence using VP4 RT-PCR and sequencing has been described for human feces (49) and sewage samples from Finland, Latvia, and Slovakia (50). Nevertheless, it is currently challenging to understand the significance of the rhinovirus species seasonal dynamics in wastewater since the concentration, extraction, and RT-PCR methods implemented in this study have not been evaluated for rhinoviruses. Further studies into the implications of the presence of rhinovirus in wastewater and into efficient methods to isolate them may be warranted. Such studies would require the use of rhinovirus-specific primers, such as those described by Guan et al. (51).

The occurrence of total EV in municipal wastewater as measured by RT-qPCR showed that densities were lowest in the winter and highest from April through October, peaking in July. This was a trend that has been reported in wastewater elsewhere (19, 24, 27, 29). Additionally, this trend is observed clinically as 78% of cases of EV infections reported to NESS also occurred from June through October (10, 12, 13).

Principal-component analysis of EV species and serotypes revealed that the monthly samples clustered according to season, with an overlap among the data from the summer and fall months. Procrustes analysis showed that species and serotype PCA plots were highly similar to each other and that 8 serotypes (i.e., CVA6, CVB3, CVB5, E9, EV97, CVA19, CVA22, and CVA1) were significant drivers of the EV seasonal trends. The spring months (March, April, and May) were predominated by EV-A, particularly CVA6, with abundances that then declined in the subsequent months of June, July, and August, suggesting a limit to the seasonal trend. CVA6 was reported to the NESS in 2010 and 2013, although it accounted for only 1.4% and 3.2% of the respective cases in those years. Most notably, infection by CVA6 was the EV infection most frequently reported to NESS in 2012, accounting for 44.7% of the 421 cases. In contrast, EV-B serotypes CVB3, CVB5, and E9 predominated in the summer, fall, and early winter months of June through December. All 3 of these EV-B serotypes were reported to NESS in 2010 and 2011. EV-C serotypes CVA19, CVA22, and CVA1 were predominant in December, January, and February. In addition, EV-B serotype EV97 was also observed during the winter. However, none of these serotypes were reported to the NESS in 2010 or 2011. This suggests that while EV-C is circulating, resulting infections do not necessarily require medical treatment or may be asymptomatic. EV-D does not appear to influence the seasonal dynamics of EV infections. EV-D serotypes EV68 and EV70 were detected in municipal wastewater but at very low abundances and sporadically throughout the sampling time frame. EV68 was reported to NESS in both 2010 and 2011, accounting for 1.4% and 5.9% of cases, respectively.

The distributions of EV observed could be a direct proportional result of loads deposited by infected individuals, including those with symptomatic and asymptomatic infections. Using either of the two wastewater processing plants studied as an example, it can be estimated on the basis of the mean RT-qPCR concentration of 2 × 105 equivalent copies/liter and the 5 million gallons of waste treated per day that 3.7 × 1012 equivalent genomic copies of EV pass through the plant daily. From this value, it can be calculated that between 12.3 and 1,237 individuals are infected in the local community, using a shedding rate from EV-infected individuals of 300 million infectious virus particles shed daily (2) and an estimation of between 10 and 1,000 genome-containing physical particles per infectious particle (52). The upper estimate of the number of infected individuals in the communities served by these WWTPs equates to about 2.8% of the local community shedding EV daily. This rate of community shedders is lower than the EV-A, EV-B, and EV-C shedding rate of 11.8% in healthy Norwegian children reported by Witsø et al. (45), but children would be expected to have a higher rate than the total population (53).

The EV distributions described here could also be reflective of differential persistence rates among EV species in municipal wastewater between households and primary treatment. Previous studies have reported variable rates of persistence among EV in wastewater sludge (54–56), but controlled studies comparing species-specific differences are lacking. Such information is necessary to help determine the relative levels of importance of seasonal shedding rates versus persistence factors (e.g., temperature, resistance to grazing, etc.) for the temporal trends in abundance observed in the present study.

It should be noted that the seasonal dynamics of EV studied here are limited to 2 WWTPs serving communities in southwest Ohio and span only a single year of monthly measurements. Although the observation of EV-B dominating wastewater in the summer-fall time frame mirrors the peak in cases reported to NESS, the seasonal trends observed here may not extrapolate to other geographical locations or time frames. To address whether EV-A, EV-B, and EV-C conform to regular cyclical patterns of peak occurrence in spring, summer-fall, and winter, respectively, additional work is necessary and could be achieved by expanding the sample locations over multiple years to discern the seasonality of the various EV species.

Here, we report the seasonal dynamics of EV infections from surveillance of municipal wastewater, using a culture-independent, deep-sequencing approach. EV-A, most notably CVA6, dominates during the spring months, while EV-B is the most abundant species in June through December, due in part to a relatively high abundance of CVB3, CVB5, and E9. EV-C changes the EV distribution profile in January and February, driven by an increased abundance of EV97, CVA1, CVA19, and CVA22. These data reveal that the EV present in municipal wastewater comprise all four species and reflect a seasonal shift in human infections from the contributing human population. Culture-independent deep sequencing is a powerful tool that could be used for studies of additional enteric viral pathogens in sewage, such as human parechoviruses, noroviruses, and adenoviruses.

MATERIALS AND METHODS

Sample collection and concentration and extraction of DNA.

Monthly wastewater samples (1 liter) were collected from the primary clarifier of two wastewater treatment plants (WWTP) in Ohio over the period of 1 year, beginning in March 2010 and ending in February of 2011. Both WWTPs treat approximately 5 to 15 million gallons of waste daily, with the majority contribution coming from residential facilities. Enteroviruses in primary effluent samples were concentrated and their genomic material extracted, as previously described (44). Briefly, nonflocculating beef extract (BD Biosciences, Sparks, MD) (1.5% [wt/vol]) and Celite 577 (Sigma-Aldrich, St. Louis, MO) (0.1% [wt/vol]) were added to each sample. After mixing, the pH of the sample was lowered to 4.0, and further mixing was then performed for 10 min at room temperature. The samples were filtered on an AP20 filter (Millipore Corporation, Bedford, MA) using a vacuum manifold to collect virus-bound Celite. AP20 filters were eluted with 30 ml of 1× Dulbecco’s phosphate-buffered saline (PBS) (United States BioLogicals, Swampscott, MA) (without MgCl2 and CaCl2; pH 9.0). The eluate pH was adjusted to 7.0 to 7.5 and then filter sterilized through a 0.22-µm-pore-size Acrodisc filter (Pall, Ann Arbor, MI). The 30-ml samples were stored at −70°.

Nucleic acids from viruses in 10 ml of concentrated samples were extracted with a QIAamp DNA Blood Maxi Extraction kit (Qiagen, Valencia, CA) according to the manufacturer’s instructions, with modifications. Buffer AVL (Qiagen) was used for chemical lysis, and protease digestion was omitted. Nucleic acids were eluted from Maxi columns with 1 ml of AE buffer (10 mM Tris-Cl; 0.5 mM EDTA; pH 9.0) augmented with 400 units of recombinant RNasin RNase inhibitor (Promega, Madison, WI). To increase yield, a second elution was performed with the volume of eluate collected from the first elution. Turbo DNA-free DNase I (Life Technologies, Inc., Grand Island, NY) was used to purify RNA genomes, as instructed by the manufacturer. DNase-treated extracts were stored at −70°C.

Quantitative reverse transcriptase PCR for enterovirus.

EV were enumerated in each sample using RT-qPCR as described above and primers targeting a conserved region of the 5′ untranslated region (forward primer, 5′-CCCTGAATGCGGCTAAT-3′; reverse primer, 5′-TGTCACCATAAGCAGCCA-3′; probe, 6-carboxyfluorescein [FAM]-ACGGACACCCAAAGTAGTCGGTTC-6-carboxytetramethylrhodamine [TAMRA]) (44). Reverse transcription was carried out in 30-µl reaction mixtures containing 5 µl of sample, 1× PCR buffer II (Life Technologies, Inc.), 1.5 mM MgCl2, 0.67 mM deoxynucleoside triphosphates (dNTPs) (Promega), 0.83 µM reverse primer, 50 U MuLV RT (Life Technologies, Inc.), and 30 U recombinant RNasin RNase inhibitor (Promega). RT reaction mixtures were incubated at 43°C for 1 h, and then RT was inactivated by heating to 94°C for 5 min. The 30-µl RT reaction mixture was included in 50-µl qPCR reaction mixtures. The final qPCR reaction mixtures contained 1× PCR buffer II, 5 mM MgCl2, 1 µl ROX dye (Life Technologies, Inc.), 0.5 µM forward primer, 0.1 µM probe, and 2.5 U AmpliTaq Gold (Life Technologies, Inc.). Reactions were processed in a 7900HT Fast real-time PCR system (Life Technologies, Inc.) at 95°C for 10 min followed by 40 cycles of 95°C for 15 s and 60°C for 1 min. To enumerate EV in wastewater samples, custom-made Armored RNA (Asuragen, Austin, TX), described in EPA method 1615 (57), was used to generate a standard curve. The custom Armored RNA was made from a construct described previously (44) containing a 195-nucleotide (nt) sequence of the 5′ untranslated region that is amplified by the primer and probe assay, along with 96-nt and 89-nt sequences of the norovirus ORF1-ORF2 region. Tenfold serial dilutions of the Armored RNA were extracted using a DNA Blood Mini Extraction kit (Qiagen). The quantities of target sequences extracted at each dilution were determined by RT-droplet digital PCR (RT-ddPCR; described below). Linear regression analysis of the standard curve was used to determine concentrations of EV in wastewater samples; standard curves exhibited amplification efficiency of at least 80% and R2 values of ≥0.97 (57).

RT-ddPCR.

Dilutions of the custom Armored RNA used for the RT-qPCR standard curve were quantified by RT-ddPCR with a Bio-Rad QX200 droplet digital PCR system. Reactions were prepared by adding 5 µl of extract to a 25-µl reaction mixture containing 600 nM primers, 250 nM probe, 0.75 mM manganese acetate (Bio-Rad Laboratories, Inc., Hercules, CA), and 2× One-Step RT-ddPCR SuperMix (Bio-Rad Laboratories, Inc.). The RT-PCR reaction mixture was sequestered into droplets using a QX200 droplet generator (Bio-Rad Laboratories, Inc.) with 20 µl of the reaction mixture and 70 µl of droplet generation oil (Bio-Rad Laboratories, Inc.). The droplets were transferred to a 96-well plate (Fisher), and RT-PCR was performed in a C1000 thermal cycler (Bio-Rad Laboratories, Inc.) using the following conditions: 60°C for 30 min and 95°C for 5 min followed by 40 cycles of 94°C for 30 s and 55°C for 1 min. At the completion of the cycling protocol, the enzymes were inactivated at 98°C for 10 min, and then the reactions were held at 4°C. The amplification products were then loaded into the droplet reader (Bio-Rad Laboratories, Inc.) to count positive and negative droplets. Quantities were determined using QuantaSoft analysis software (Bio-Rad Laboratories, Inc.). The 10-fold dilution series of the custom Armored RNA standards consisted of 4.3 × 104, 4.3 × 103, 4.3 × 102, 4.3 × 101, and 4.3 × 10° equivalent copies (per reaction).

RT-PCR of VP1 and VP4 and amplicon isolation for sequence analysis.

DNase-treated samples were used in 2-step RT-PCR reactions to amplify part of the VP1 gene (approximately 340 nt) and the entire VP4 gene (approximately 650 nt). VP1 primers previously described by Oberste et al. (46) were modified by adding sequence adapters for use with Roche’s Lib-A GS FLX Titanium emulsion-based PCR (emPCR) system (forward primer [primer A+292], CGTATCGCCTCCCTCGCGCCATCAGMIGCIGYIGARACMGG; reverse primer [primer B+222], CTATGCGCCTTGCCAGCCCGCTCAGCICCIGGIGGIAYRWACAT). VP4 primers described by Ishiko et al. (37) were modified with the adapter sequences for use in a Lib-L GS FLX Titanium emPCR system (forward primer [primer A+MD91F], CCATCTCATCCCTGCGTGTCTCCGACTCAGCCTCCGGCCCCTGAATGCGGCTAAT; reverse primer [primer B+OL68-1R], CCTATCCCCTGTGTGCCTTGGCAGTCTCAGGGTAAYTTCCACCACCANCC). RT reactions were prepared with 10 µl of DNase-treated nucleic acid extracts, 1× reaction buffer (Roche), 1.5 mM MgCl2, 0.5 mM dNTPs, 1.6 µM reverse primer, 32 U of murine leukemia virus (MuLV) RT (Life Technologies, Inc.), and 20 U of recombinant RNasin RNase inhibitor (Promega). The 20-µl reaction mixtures were incubated at 43°C for 1 h and then at 94°C for 5 min. The RT reaction mixture was added to a 50-µl PCR reaction mixture comprised of 1× reaction buffer, 2.6 mM MgCl2, 0.2 mM dNTPs, 2 µl dimethyl sulfoxide (DMSO), a 0.5 µM concentration of each primer, and 2.5 U high-fidelity enzyme blend (Roche). For amplification of the 5′ end of VP1, PCR cycling conditions were 94°C for 5 min and 40 cycles of 94°C for 30 s, 42°C for 30 s, and 72°C for 45 s with a ramp time of 0.6°C/s. Extension continued at 72°C for 5 min. VP4 was amplified by heating to 95°C for 5 min followed by 40 cycles of 95°C for 30 s, 55°C for 30 s, and 72°C for 1 min and extension at 72°C for 7 min. Amplicons were collected using an E-gel Size Select system (Life Technologies, Inc.), and replicate sample amplicons were pooled. Positive controls for PCR were prepared from stocks of a serotype to represent the 4 targeted EV species (for EV-A, coxsackievirus A7; for EV-B, coxsackievirus B1; for EV-C, poliovirus 1; for EV-D, enterovirus 70). Stocks were prepared from infection of Buffalo green monkey kidney cells as previously described (44).

Quality assurance/quality control for RT-PCR.

The guidance proposed by Sen et al. (58) was strictly enforced to minimize potential PCR contamination. Separate and dedicated laboratories were used for preparation of master mixes, nucleic acid extraction and reaction preparation, PCR amplification, and post-PCR analysis.

Sequencing of purified amplicons.

VP1 amplicons were pooled by month, and the April-November samples were prepared with a Lib-A GS FLX Titanium emPCR kit (Roche) according to the manufacturer’s instructions. VP4 products were pooled by month, and libraries were prepared with a Lib-L GS FLX Titanium emPCR kit (Roche) as instructed by the manufacturer. VP1 and VP4 libraries were sequenced with a 454 GS FLX kit (Roche) according the manufacturer’s instructions.

Bioinformatic and statistical analysis.

The raw reads were processed with Mothur (59) and USEARCH (60). Reads were filtered for quality (Phred quality score of 30 with window length of 50 nucleotides), homopolymers (maximum of 8 nucleotides), length (minimum of 300 nucleotides), and ambiguous base calls (0 allowed). Filtered reads were clustered at 100% sequence identity to find unique sequence reads that were then used in discontinuous megablast analysis (61) of the Virus Reference Database in GenBank for species assignment. The sequence reads were uploaded to MG-RAST (62) and are available as project number 12933 using the search icon at http://metagenomics.anl.gov. Serotype assignments were made using best-hit classification and the M5NR annotation source in MG-RAST (minimum percentage identity cutoff of 60%, maximum E value cutoff of 1e−5, and minimum alignment length cutoff of 15 bp). The sequence read archives are available at NCBI.

Statistical comparisons among EV species measured each month by RT-qPCR were assessed by repeated-measure ANOVA using SigmaPlot (version 12.5, Systat Software, Inc., Chicago, IL). Correlation analysis of quantities determined by RT-qPCR and calculations of the number of VP1 and VP4 reads obtained for each sample were also performed with SigmaPlot. The nonparametric Mann-Kendall test was performed on the time series to determine if the read abundances were unidirectional. This test and cross-correlation analyses were performed using Past3 version 3.13 software (63). Principal-component analysis (PCA) using a variance-covariance format and a matrix plot was performed using Past3. Procrustes analysis and Mantel correlation analysis were performed with the vegan package for R (64) using Euclidean distances.

Accession number(s).

The sequence read archives are available at NCBI under BioProject identifier (ID) PRJNA384418.