The phylogenetic and phylodynamic analyses reported here have demonstrated key features of two concurrent North American EV-D68 outbreaks and EV-D68’s epidemiology more widely. Through phylogenetic approaches we have demonstrated the epidemiological connection between the 2022 EV-D68 outbreaks in Maryland and Ontario. Our use of phylodynamic methods have aided in narrowing down the plausible window for EV-D68’s infection period. Furthermore, we have illustrated the practicality of phylodynamic methods in deriving Rt and epidemic origin. Notably, our epidemic origin estimates for the 2022 EV-D68 outbreaks coincide with the removal of restriction aimed at curtailing the spread of COVID-19.
Cross-border disease dynamics between North American countries have been studied for other viral pathogen such as SARS-CoV-2, mumps virus and West Nile viruses42–44. In our analysis we demonstrate close genetic proximity between viruses circulating in Maryland and Ontario in 2022, where isolates captured from each outbreak cluster within the same sub-clade, B3, and share recent common ancestry. These viruses belong to the B3 clade of EV-D68 lineages which have previously been shown to play an important role in outbreaks in the region13,45. The genetic proximity of these isolates indicates that significant epidemiological connections exist between these regions.
Using phylodynamic modelling we were able to estimate important epidemiological parameters: the infection period and the time-varying reproduction number (Rt). Our models showed highest support for a duration of infection of 7.94 days (95% HPD 4.55, 12.8 days) for the ON-2022 dataset (n = 87) and 10.8 days (95% HPD 5.85, 18.6 days) for the MD-2022 dataset. Prior estimates of the infection period are lacking, however, respiratory viral shedding for enteroviruses has been documented to be between 1–3 weeks21. Specifically, the infection period for the EV-D68 ranges from 6–10 days25, and a symptomatic period from 3–10 days26. More recently Nguyen-Tran et al., (2023)8 found that the EV-D68 genome could be detected in the upper respiratory tract for a median of 12 days post symptom onset (7–15 days). Nguyen-Tran et al., (2023)8 point out that these RNA detection period should only be seen as an upper limit for infectious period. Given that infectious period (unlike infection period) does not include the latency period, Nguyen-Tran et al.,’s (2023)8 findings are concurrent with our infection period of 7.94–10.8 days. The concurrence between our findings and Nguyen-Tran et al.,’s (2023)8 brings important specificity, which is clinically relevant in the management of patients, given the previously suggested broad range in infection periods for EV-D68.
Our time-varying reproduction number (Rt) estimates derived through phylodynamic methods produced similar values compared to deriving Rt through case count data and the serial interval (Fig. 4). Of particular note, greater concordance was seen when case count derived Rt values used a higher serial interval (mean = 7 days, SD = 2.6). Our mid-point serial interval for EV-D68 (mean = 3.6 days, SD = 2.6) was based on EV-7141, but given our estimate of the EV-D68 infection period (7.94–10.8 days) and an upper limit for the infectious period of 12 days8, EV-D68’s serial interval is likely to be closer to the higher value used in our sensitivity analysis. Previous epidemiological estimates of Rt from EV-D68 outbreaks (across several US states 2014-17) range between 0.5-1.622. We find that our estimates of the median Rt (Fig. 4) were just over 1 in non-epidemic periods and 2.70 (95% HPD 1.76, 4.08) in Ontario and 2.10 (95% HPD 1.41, 3.17) in Maryland during the respective peak epidemic periods. A build up in the susceptible population due to reduced contacts over 2020–2021 may have led to the increased Rt values observed in Ontario and Maryland 2022, compared to estimates from several US states over 2014-17, a pre-pandemic period22. More generally, our EV-D68 BDSKY derived Rt estimates are consistent with other respiratory pathogens particularly other enteroviruses24,46–49. As with Park et al., (2021) 22, we found delays in the increase of Rt were associated with outbreaks occurring farther north within North America.
The WGS-based substitution rates reported here, 0.0148 substitution per site per year (95% HPD 0.0112, 0.0185) for ON-2022 and 0.0113 substitution per site per year (95% HPD 0.0078, 0.0154) for MD-2022, are substantially higher than reported previously, 0.003 substitution per site per year11. The 38 WGS used in the analysis Eshaghi et al., (2017) 11 come from 14 different countries, span 1960–2014 and therefore come from different EV-D68 clades, whereas, the WGS being used in our analyses are from one regional outbreak of a single sub-clade, B3. Time varying evolutionary metrics have been observed before, with faster rates observed when samples are drawn from shorter time periods50–52. Ghafari et al., (2022)53 demonstrated that during the SARS-CoV-2 and pH1N1 influenza pandemics this time varying evolutionary rate could be attributed to a short-term buildup of mildly deleterious mutations, that were eradicated over a longer term through purifying selection. This process of incomplete purifying selection may be the reason for the discrepancy between the EV-D68 substitution rates reported here and earlier11.
The estimated 2022 EV-D68 epidemic origin and TMRCA statistics from BDSS models for both regions coincide with the periods when measures to reduce social contact, known as Non-Pharmaceutical Interventions (NPIs), were relaxed. Specifically, in Ontario, Canada, the period was from January 31, 2022 to March 14, 2022, coinciding with the decline of the Omicron COVID-19 wave54. Meanwhile, the Maryland data corresponded with the phase of winding down several NPIs aimed at curtailing the spread of COVID-19, from early February 2021 to August 13, 202155–57. NPIs aimed at controlling COVID-19 transmission have also significantly reduced influenza cases, virtually eliminated respiratory syncytial virus (RSV) hospitalization and diminish detectable circulation of several enteroviruses58–61. Therefore, it is possible that the coinciding of our epidemic origin and TMRCA estimates with the lifting of NPIs demonstrates the suppressing effect of NPIs on EV-D68 transmission. However, Fig. 1 depicts 2022 EV-D68 Maryland and Ontario sequences interspersed with each other and sequences from Sweden. This pattern suggests that the 2022 EV-D68 outbreaks in Ontario and Maryland may be the result of several independent introductions into their respective populations, and not a single introduction. This would mean that our Rt estimates are more likely to be for EV-D68 outbreaks in regions greater than Ontario or Maryland the further back in time the estimate is. Likewise, this may mean that our TMRCA and origin estimates are for EV-D68 outbreaks occurring over a much wider region than Ontario or Maryland.
This study has limitations that should be addressed in further research efforts. For instance, the above caveats over Rt, TMRCA and origin estimates have, in part, come about through sampling in acute healthcare settings during an ongoing transmission within the wider community (Fig. 4C). It is important that sampling efforts are broader and capture more localities nationally, as well as broadly in North America, if not globally. As seen in the wider phylogenetic analysis (Fig. 1) there are long branches across the phylogeny which may indicate prolonged periods of within-host evolution, missed infections or un-sampled diversity. Thus, active surveillance is critical in identifying major source and sink populations for the EV-D68 virus, directing intervention efforts effectively. In addition to sampling biases, it is important that clinical observation studies of positive cases are conducted to validate the in-silico estimates of infection period for EV-D68 viruses to robustly model epidemiological dynamics further.
Future study of EV-D68 in a phylodynamic framework will not only be bolstered by wider sampling efforts but will also be aided by the inclusion of secondary metadata to study the importance of different host traits on viral evolution and diffusion. If metadata pertaining to severity of infection, age, and travel history of a patient is available phylodynamic methods can be used to determine the importance of traits in the diffusion process and potentially identify host characteristics that can inform control measures62,63. In summary, this study underscores the importance of pathogen genome surveillance combined with phylodynamics in complementing conventional epidemiological approaches within public health investigations.