A complex mosaic of enteroviruses shapes community-acquired hand, foot and mouth disease transmission and evolution within a single hospital.

Human enteroviruses (EV) pose a major risk to public health. This is especially so in the Asia-Pacific region where increasing numbers of hand, foot and mouth disease (HFMD) cases and large outbreaks of severe neurological disease associated with EV-A71 have occurred. Despite their importance, key aspects of the emergence, epidemiology and evolution of EVs remain unclear, and most studies of EV evolution have focused on a limited number of genes. Here, we describe the genomic-scale evolution of EV-A viruses sampled from pediatric patients with mild disease attending a single hospital in western Sydney, Australia, over an 18-month period. This analysis revealed the presence of eight viral serotypes-Coxsackievirus (CV) A2, A4, A5, A6, A8, A10, A16 and EV-A71-with up to four different serotypes circulating in any 1 month. Despite an absence of large-scale outbreaks, high levels of geographical and temporal mixing of serotypes were identified. Phylogenetic analysis revealed that multiple strains of the same serotype were present in the community, and that this diversity was shaped by multiple introductions into the Sydney population, with only a single lineage of CV-A6 exhibiting in situ transmission over the entire study period. Genomic-scale analyses also revealed the presence of novel and historical EV recombinants. Notably, our analysis revealed no association between viral phylogeny, including serotype, and patient age, sex, nor disease severity (for uncomplicated disease). This study emphasizes the contribution of EV-A viruses other than EV-A71 to mild EV disease including HFMD in Australia and highlights the need for greater surveillance of these viruses to improve strategies for outbreak preparedness and vaccine design.

1. Introduction

Human enteroviruses (EV) are important viral pathogens that pose major public health risks and substantial disease burden, particularly in the Asia-Pacific region. EV infections are highly contagious and manifest as mild self-limiting disease with symptoms including fever, malaise and hand, foot and mouth disease (HFMD). HFMD is characterized as a non-itchy skin rash on the palms of the hands and soles of the feet, which presents as erythematous papules, some of which may form blisters, and small ulcers in the mouth and throat that can progress to herpangina—the coalescence of multiple ulcers on the soft palate and posterior pharynx. The annual incidence of HFMD in China is reported to be 1.2 per 1,000, with an estimated 7.2 million cases between 2008 and 2012 (Xing et al. 2014). A small proportion of EV infections are associated with severe neurological disease, including meningo-encephalitis, which can result in long-term neurological morbidity and in rare cases death. Between 2012 and 2016, an average of 5.63 per cent (ranging from 2.81 to 9.55) of HFMD cases in Shanghai, China, resulted in severe illness (Wang et al. 2018a).

RNA viruses of the genus Enterovirus belong to the family Picornaviridae that is divided into fifteen species including poliovirus, rhinovirus and enterovirus (EV) A71 (EV-A71). Members of the EV A (EV-A) species are considered the main etiological agent of HFMD and are responsible for more than 90 per cent of reported cases (Mirand et al. 2012). Currently, there are twenty-five identified EV-A serotypes, nineteen of which solely infect humans, including Coxsackievirus (CV) A2-8, 10, 12, 14, 16 and EV-A71. EV-A viruses are positive-sense single-stranded RNA viruses approximately 7.4 kb in length. EV-A genomes consist of a single open reading frame, flanked by 5′ and 3′ untranslated regions (UTR). The genome is initially translated into a single polyprotein of ∼2,185 amino acids, encoding eleven proteins. The polyprotein is then cleaved into three precursor proteins: P1 that is further cleaved into VP1–VP4 encoding the capsid proteins and the P2 and P3 segments that encode the non-structural proteins, including the viral polymerase (3Dpol gene). As with all RNA viruses, the Picornaviridae experience high rates of evolutionary change (Jenkins et al. 2002), reflecting frequent mutation per replication (Sedivy et al. 1987; de la Torre et al. 1992). In addition, recombination at both the intra- and inter-serotypic levels is commonplace and an important source of genetic variation (Kyriakopoulou et al. 2015). Although recombination is considered in Enterovirus classification at the species level, classification at the serotypic level is defined by the genetic diversity of the capsid genes alone, such that the impact of recombination is ignored (Lukashev 2005). The loss of such large amounts of genetic information in serotype classification may in part explain our current inability to connect viral serotype with disease profile.

Historically, the coxsackievirus A16 (CV-A16) and EV-A71 serotypes have been considered the primary agents of HFMD (Xing et al. 2014). More recently, however, they have been displaced by the CV-A6 serotype that has become the prominent EV-A serotype associated with major HFMD outbreaks worldwide (Bian et al. 2015). Additional serotypes, including CV-A2, 4, 5, 8, 10, 12 and 14, are also responsible for sporadic cases of HFMD globally (Park et al. 2011; Lukashev et al. 2014; Chen et al. 2016). Although CV-A6 is the most common causal agent of HFDM outbreaks, EV-A71 is still considered the greatest threat due to its association with severe disease including large outbreaks of encephalitis and high levels of mortality. Large EV-A71 outbreaks include those in Taiwan in 1998 (Ho et al. 1999), China in 2008 (Zhang et al. 2010) and Vietnam in 2011 (Nguyen et al. 2014). In Australia, four EV-A71 outbreaks associated with severe illness have been reported in Victoria (1986) (Gilbert et al. 1988), Perth (1999) (McMinn et al. 1999) and two in Sydney (2000–2001 and 2013) (Prager et al. 2003; Zander et al. 2014). In the most recent 2013 Sydney outbreak, 119 cases of suspected EV-A71 infection were identified in the Sydney Children’s Hospital Network, a pediatric healthcare network incorporating Sydney’s two tertiary/quaternary children’s hospitals and associated services (Zander et al. 2014). Viruses of the CV-A2 (Chen et al. 2010; Yip et al. 2013; Yen et al. 2017) and CV-A4 (Chen et al. 2018) serotypes have also been described as causing a small number of severe disease cases.

The significant population growth of Sydney, 18.2 per cent in the 10 years between 2006 and 2016 (http://www.abs.gov.au/census), and the dramatic increase in the number of Australian border crossings (73% increase in the same time frame; http://www.abs.gov.au/census), contribute to a rise in EV-A introductions into Australia and potentially an increase in epidemic frequencies. Despite the disease burden of EV-A71 globally, little is known about the threat that these viruses pose to Australia, with those studies performed so far focusing primarily on routine serotyping (McMinn et al. 2001; Sanders et al. 2006). While HFMD is considered a common viral infection in Australia, no detailed incidence data is available and as such it remains unclear what EV-A serotypes routinely circulate in Australia and cause sporadic HFMD outbreaks. Additionally, it is not known whether the EV-A71 vaccines currently licensed in China could provide protection against severe disease caused by EVs in Australia. To address these questions, we determined the genomic composition of EV-A infections in children who had laboratory confirmation of EV infection after presenting to The Children’s Hospital at Westmead (CHW), in western Sydney, Australia, over an 18-month period (February 2016–July 2017). Westmead is located 26 km west of Sydney’s central business district and has a high level of cultural diversity with nearly 75 per cent of residents born outside of Australia and only 20 per cent of households having English as the only language spoken at home. From these data, we described a wide diversity of EV-A species collected through a single hospital and investigated their potential source populations as well as the association of specific viruses with different disease outcomes.

2. Materials and methods

2.1 Sample collection

This study was performed under the ethical approval of the Children’s Hospital at Westmead Ethics Committee (Approval number SCHN LNR/14/SCHN/528). Clinical samples used in this study were collected from patients who presented to CHW and tested positive for picornaviruses (EV or human paraechovirus) by the diagnostic laboratory at CHW. All samples were collected between February 2016 and July 2017. Following collection, ‘retrospective opt-out’ consent was obtained from participants or their guardians.

2.2 Disease severity classification

Samples were grouped according to the severity of enteroviral disease using a schema based on those for EV-A71 by Huang et al. (2017) and the WHO HFMD guidelines (WHO 2011). Accordingly, patients were classified into five stages: (i) Stage 0A: febrile illness or respiratory tract infection, and EV detected from a non-sterile site; (ii) Stage 0B: febrile convulsions and EV detected from a non-sterile site; (iii) Stage 1: HFMD or febrile illness with rash or herpangina, and EV detected from a non-sterile site; (iv) Stage 2: aseptic meningitis (detected from cerebrospinal fluid (CSF) with or without CSF pleocytosis) or myoclonic jerk and EV detected from any site; (v) Stage 3: complicated CNS infection (EV detected from any site, and seizures (excluding febrile convulsion) or focal neurological signs or altered LOC/behavior), with specific syndrome designated as encephalitis (Stage 3A) or cardio-pulmonary failure (Stage 3B); and (vi) Stage 4: other severe organ disease including, but not limited, to hepatitis or myocarditis. The non-sterile sites sampled included feces, nasopharyngeal aspirate, throat swab or skin, vesicles or ulcer swabs, with CSF the only sterile site sampled.

2.3 Viral genome sequencing and assembly

Total nucleic acid was isolated from 140 μl of clinical material using the QIAamp viral RNA kit (Qiagen) and eluted in 50 μl of water and stored at −80°C. Two overlapping fragments were amplified to recover complete viral genomes. First, Super-Script III One-Step reverse transcription polymerase chain reaction (RT-PCR) System with PlatTaq HiFi was used to amplify the last ∼3 kb of the genome (Invitrogen), using primers designed to recognize regions conserved between all EV-A viruses, EV-2F1 and EV-2R, the sequences of which are given in Supplementary Table S1. For amplification of the first ∼4.5 kb of the genome, a two-step RT-PCR process was used with primers designed from the sequences obtained. cDNA was generating using the SuperScript IV VILO system, as per manufacturer’s instructions (Invitrogen). The first 4.5 kb were subsequently amplified with Platinum SuperFi DNA Polymerase. The forward primer (EV-1F1) was designed to the highly conserved 5′ UTR, but due to the high level of diversity between samples, reverse primers were designed specifically to the samples, EV-1R1 to EV-1R39, shown in Table 1.

Table 1.

Epidemiological data of pediatric patients presenting to the Children’s Hospital at Westmead (Sydney, Australia) between February 2016 and July 2017 from which viral genomes were obtained.

Sample #	Age (months)	Sex	Collection date (d/m/y)	Serotype	Disease severity
C002	13	F	29/02/16	CV-A6	Stage 0A
C007	15	M	16/03/16	CV-A6	Stage 1
C019	22	M	05/04/16	CV-A4	Stage 1
C026	10	M	18/04/16	CV-A5	Stage 0A
C027	123	F	27/04/16	CV-A5	Stage 1
C028	16	M	22/04/16	CV-A16	Stage 0A
C031	12	M	06/05/16	CV-A6	Stage 1
C037	14	M	28/05/16	CV-A4	Stage 0A
C038	10	F	30/05/16	CV-A6	Stage 1
C040	14	M	10/06/16	CV-A4	Stage 0B
C045	53	F	23/06/16	CV-A6	Stage 1
C047	11	F	23/06/16	CV-A6	Stage 0B
C048	63	M	28/06/16	CV-A6	Stage 1
C051	16	M	14/07/16	CV-A6	Stage 1
C056	66	M	02/08/16	CV-A6	Stage 1
C060	17	M	18/08/16	CV-A6	Stage 1
C062	20	F	11/08/16	CV-A4	Stage 0B
C067	86	M	23/08/16	CV-A6	Stage 0A
C068	19	M	14/08/16	EV-A71_C2	Stage 2
C069	25	M	24/08/16	CV-A6	Stage 0B
C071	10	M	28/08/16	CV-A6	Stage 0A
C074	11	M	24/08/16	CV-A4	Stage 0A
C079	5	F	09/06/16	CV-A4	Stage 1
C080	33	M	05/09/16	CV-A4	Stage 0B
C082	20	M	05/09/16	CV-A4	Stage 0B
C090	61	F	14/09/16	CV-A2	Stage 0A
C096	25	F	22/09/16	CV-A6	Stage 1
C098	181	M	22/12/16	CV-A6	Stage 0A
C099	10	F	29/10/16	CV-A6	Stage 1
C104	11	M	12/01/17	CV-A6	Stage 1
C105	20	F	18/01/17	CV-A16	Stage 1
C107	11	F	18/09/16	CV-A16	Stage 1
C108	33	F	09/10/16	CV-A2	Stage 0B
C109	2	F	02/12/16	CV-A6	Stage 1
C115	24	M	23/09/16	CV-A6	Stage 1
C116	5	M	18/08/17	CV-A10	Stage 0A
C123	7	F	17/02/17	CV-A16	Stage 1
C125	18	F	22/02/17	CV-A10	Stage 0A
C131	17	F	06/02/17	CV-A2	Stage 0A
C136	7	M	13/12/16	CV-A10	Stage 1
C138	2	M	14/12/16	CV-A16	Stage 1
C140	30	M	04/11/16	CV-A10	Stage 1
C142	1	F	22/10/16	CV-A10	Stage 2
C146	52	F	04/04/17	CV-A4	Stage 0B
C149	15	M	02/04/17	CV-A10	Stage 0A
C150	4	M	14/04/17	CV-A10	Stage 1
C151	11	M	16/04/17	CV-A16	Stage 0A
C152	30	M	27/03/17	CV-A2	Stage 0B
C157	46	M	12/03/17	CV-A10	Stage 0A
C158	23	M	17/04/17	CV-A6	Stage 1
C160	22	M	09/04/17	CV-A6	Stage 1
C162	12	M	29/03/17	CV-A10	Stage 0A
C177	15	F	15/05/17	CV-A8	Stage 0A
C179	8	F	25/04/17	CV-A4	Stage 0A
C183	48	M	09/07/17	CV-A6	Stage 0A
C184	22	F	29/06/17	CV-A8	Stage 0B
C185	6	F	02/07/17	CV-A6	Stage 1
C191	145	F	14/07/17	CV-A6	Stage 0A

DNA libraries of each amplicon separately were produced using the NexteraXT DNA sample preparation kit (Illumina) and libraries were confirmed by quantitative PCR using the universal KAPA library quantification kit for Illumina platforms (Kapa Biosystems). Quality and fragment size were then estimated using the Agilent 2100 Bioanalyser. Equimolar amounts of each sample were pooled for 2 × 150 bp paired-end sequencing performed on the MiSeq platform using a 300-cycle reagent kit (v2), at the Australia Genome Research Facility.

Geneious 7.1.9 (http://www.geneious.com; Kearse et al. 2012) was used to analyze and process sequencing reads. First, forward and reverse barcode reads were paired. Sequences were then trimmed to remove poor-quality data according to their quality scores using the modified Mott algorithm, with an error probability cut-off of 0.001 and at least 15 and 2 bp at the 5′ and 3′ ends were removed, respectively. Trimmed and filtered sequences <100 bp were discarded.

A total of 10,000 random reads from each library were sampled and assembled de novo. Consensus sequences were extracted and an alignment of the two amplicons was undertaken to obtain a full-length consensus sequence. Libraries of both amplicons were then combined and reads remapped to the consensus sequence.

2.4 Phylogenetic analysis

Overall, fifty-eight complete genomes (i.e. consensus sequences) of EV-A viruses were sequenced as part of this study. Preliminary genotyping was undertaken using the Enterovirus Automated Genotyping Tool (Kroneman et al. 2011). Full-length genomes were constructed containing the entire genome with the exception of the first and last 23 nt of the 5′UTR and 3′UTR. The missing nucleotides represented the conserved sequences used for primer design. For more in-depth evolutionary analysis (see Results), the full-length viral genome, the VP1 gene sequence and 3Dpol gene sequence were aligned individually using MAFFT (Kuraku et al. 2013) and inspected manually, producing data sets of 7,612 bp, 894 bp and 1,386 bp in length, respectively. For all three data sets the optimal nucleotide substitution model identified using jModelTest (Posada 2008) was GTR+I + Γ4, which was used in the subsequent phylogenetic analysis. Phylogenetic trees were estimated using the maximum likelihood (ML) procedure in PhyML 3.0 (Guindon and Gascuel 2003). A combination of SPR and NNI branch-swapping was used to estimate the ML tree, and bootstrap support values were generated using 1,000 replications.

To place our specimens within the context of viruses sampled globally, phylogenetic analysis using the full-length VP1 gene was undertaken for each of the coxsackievirus serotypes, with background sequences compiled from GenBank (http://www.ncbi.nlm.nih.gov/genbank/). After the removal of duplicate sequences this resulted in data sets of: CV-A2 (n = 89 sequences), CV-A4 (n = 132), CV-A5 (n = 15), CV-A6 (n = 1589), CV-A8 (n = 20), CV-A10 (n = 322) and CV-A16 (n = 1913). Initially, the sub-genogroup of the single EV-A71 specimen was assigned using the Enterovirus Automated Genotyping Tool (Kroneman et al. 2011) before the same phylogenetic procedure as described above was performed on only EV-A71 C2 sub-genogroup sequences, resulting in a data set of n = 196 sequences. The very large size of some of the VP1 data sets (i.e. up to 1,913 sequences) meant that ML phylogenetic trees were necessarily estimated using the RAxML program (version 8.2.4) (Stamatakis 2006), employing the GTR + Γ4 nucleotide substitution model and 100 bootstrap replicates.

2.5 Recombination analysis

An analysis of putative recombination events within our EV data set was performed using the RDP, GENECOV and Bootscan methods available within the Recombination Detection Program v4 (RDP4) program (Martin et al. 2015). Default parameters were used in all cases with a window size of 30 pb.

2.6 Phylogenetic trait analysis

To determine whether sequences with particular phenotypic traits, specifically disease severity, patient age and patient sex, clustered together more often on phylogenetic trees than might be expected by chance alone, we used the Bayesian approach in MrBayes (Ronquist et al. 2012) to infer phylogenies of the VP1 gene, full-length genomes and 3Dpol gene (using the same substitution models as in our PhyML analyses). Chain length was set to 1,000,000 for VP1 and 3Dpol, or 2,000,000 for full-length genomes, sampling every 500 with the first 10 per cent of generations discarded as burn-in. Convergence and effective sample size were examined using Tracer to confirm that estimated sample size exceeded 200. To determine the association significance between a particular disease characteristic and its distribution on a phylogeny, and hence including serotype, we performed a series of phylogeny-trait tests using the Association Index, Parsimony Score and maximum clade size statistic using the Bayesian Tip-Significance (BaTS) testing program (Parker et al. 2008) that corrects for phylogenetic uncertainty by using a posterior distribution of trees. Statistical significance was defined as P < 0.05 and the program was run for 1,000 replications.

2.6 GenBank accession numbers

All consensus sequences generated here have been submitted to GenBank and assigned accession numbers MH111016 to MH111073.

3. Results

3.1 Demographic characterization of EV-A viruses in patients presenting to the Children’s Hospital at Westmead

Between February 2016 and July 2017 a total of 8,246 samples were tested, of which 328 tested positive for EV at the diagnostic microbiology laboratories at CHW. Of these, 174 samples, collected in an unbiased manner, were salvaged and analyzed for this study (Supplementary Fig. S1). Demographic and disease data were collected with consent for 170 patients who tested positive for EV. The male-to-female ratio of infections was 1.15:1. The average age of patients was 25 months (range, 8 days to 15 years). EV was detected from a non-sterile site in 60 per cent of patients, 52 per cent with non-specific febrile illness (i.e. Stage 0A) and 8 per cent with febrile convulsions (Stage 0B). Of the remaining 40 per cent, 25 per cent had EV only detected from a non-sterile site but with HFMD, herpangina or febrile illness with a rash (Stage 1), and 12 per cent with aseptic meningitis or myoclonic jerk and EV detected in a sterile or non-sterile site (Stage 2). Two patients presented with complicated CNS infections which resulted in cerebellitis (Stage 3) and the remaining three patients progressed to Stage 4 with severe organ disease including hepatitis, hepatic failure and pancreatitis.

Of the 174 EV samples collected, full-length human EV-A genomes were produced from fifty-eight different uncultured patient samples. These were all from children with mild (Stage 0A–2) disease (Table 1). Using the Enterovirus Automated Genotyping Tool (Kroneman et al. 2011) that utilizes the consensus sequence of the VP1 gene, eight unique serotypes were identified: CV-A2 (n = 4), CV-A4 (n = 10), CV-A5 (n = 2), CV-A6 (n = 24), CV-A8 (n = 2), CV-A10 (n = 9), CV-A16 (n = 6) and EV-A71 sub-genogroup C2 (n = 1) (Table 1).

3.2 High levels of geographical and temporal mixing of EV serotypes

To determine the seasonal patterns of EV infections at CHW and identify any unusual spikes in cases, the number of EV-positive samples was examined each month. Over the period of the study between 2 and 20 EV-positive samples were collected each month (Fig. 1). A spike of thirty-eight cases was observed at the end of winter 2016 (August and September) initially indicating a potential outbreak. However, serotype characterization revealed the presence of five different serotypes (EV-A71, CV-A6, CV-A4, CV-A16 and CV-A2) circulating during this period, suggesting that this elevated number of infections was a result of multiple viruses and not a single viral outbreak. Viruses of the CV-A6 serotype were present throughout the entire study suggesting potentially sustained transmission within Sydney, while CV-A10 appeared to enter the population, circulate for a limited time, and then disappear. To identify localized outbreaks, the residential postcode of patients was acquired and mapped. From this it was visually apparent that there were no centralized outbreaks with single serotype infections in any locality (Fig. 2). Together, these data suggest a surprisingly high level of serotype mixing both geographically and temporally. Strikingly, while a high number of serotypes were identified throughout this study, we found no evidence of mixed infection.

Figure 1.

Epidemiological details of laboratory confirmed enteroviruses at the Children’s Hospital at Westmead (CHW), Sydney, Australia for which complete viral genomes could be obtained. The number of laboratory confirmed EV cases collected from patients who presented to CHW and were tested for EV between February 2016 and July 2017 is shown. The number of each serotype is shown and colored according to the key provided, un-serotyped infections are presented in white.

Figure 2.

Maps showing the residential postal codes within Sydney, Australia of patients with laboratory confirmed EV at the Children’s Hospital at Westmead (CHW) between February 2016 and July 2017. The black circle depicts the location of CHW.

3.3 Diverse EV-A viruses circulate in Western Sydney

To characterize the genetic diversity and evolutionary relationship of the EV-A viruses sequenced in this study, we inferred phylogenetic trees of the VP1 gene that encodes the main antigenic sites and from which serotypes are defined. Accordingly, the EV-A viruses sampled in this study formed eight well-supported distinct monophyletic groups containing the prototype strain of the corresponding serotype (Fig. 3A). Notably, within each serotype (i.e. clade), the viruses from CHW had only 75.4 per cent–84.7 per cent nucleotide sequence similarly to the corresponding prototype viruses in the VP1 gene, confirming that they are evolutionarily divergent from the prototype strains.

Figure 3.

Phylogenetic trees based on the (A) VP1 gene, (B) full-length genome and (C) 3Dpol gene of EV-A. Prototype strains of each serotype downloaded from GenBank were included as reference strains. For those viruses sampled from children who presented to CHW, the first number indicates the sample reference and the last number indicates the date of collection (d/m/y). Branches were colored by sample serotype and according to the key provided. Asterisks (*) represent statistical support for individual nodes >70 per cent estimated from 100 bootstrap replicates. Scale bars are proportional to the number of nucleotide substitutions per site. Trees were mid-point rooted for clarity.

To further investigate the evolutionary relationship of the viruses isolated at CHW, we performed additional phylogenetic analyses using the full-length genome sequences (Fig. 3B). Viruses isolated in this study again formed well-supported monophyletic groups at a full-length genome scale, but the relationship between the serotypes broke down, with the CV-A2 and CV-A6 prototype strains not clustering with those CHW samples of the same serotype. The phylogenetic differences between the full-length genome and the VP1 genes suggest a high level of genetic diversity between the prototype strains and the CHW isolates exists outside of the VP1 gene and that there has been frequent recombination. To address this possibility, we inferred phylogenetic trees of the highly conserved 3Dpol gene locate in the P3 region at the opposite end of the viral genome to the VP1 gene (Fig. 3C). The phylogenetic relationship among the 3Dpol genes was markedly different to that of both VP1 and the full-length genome. In particular, two well-supported groups were identified each with a mixed serotype composition and clustering. The clear clustering of viruses isolated from CHW of the same serotype disappeared and viruses often occupied different phylogenetic positions to those observed in VP1. While there were examples where the CHW isolates clustered based on serotype, such as CV-A2 and CV-A10, viruses of the CV-A4 and CV-A6 serotypes appeared to fall in to two and three distinct lineages, respectively. The inconsistencies between the VP1 and 3Dpol phylogenies highlight the mosaic nature of EV-A viruses and strongly supports a history of inter-serotypic recombination.

To analyze these recombination events in more detail we compared whole EV genome sequences using the RDP4 program (Martin et al. 2015). Accordingly, the RDP, GENECONV and Bootscan methods consistently identified eight putative recombination events with strong statistical support, of which two were newly identified here (Supplementary Table S2, Supplementary Fig. S2). These novel recombination events were inferred to occur between viruses of the CV-A4 and CV-A10 serotypes, and between the CV-A6 and CV-A4 serotypes.

Phylogenetic analysis also allowed the identification of potential localized transmission events in which very closely related groups of viruses were sampled within narrow time periods. For example, a small monophyletic group of closely related CV-A4 viruses (C074, C079 and C080) were sampled within a month of each other, indicative of localized transmission (Fig. 3B). However, a lack of additional sampling indicates that this cluster did not lead to a substantial outbreak. In addition, two pairs of CV-A6 viruses (C045 and C051, and C069 and C067) were sampled within a short time period and similarly did not result in large-scale disease outbreaks. A monophyletic group of viruses was also observed in the CV-A10 serotype (Fig. 3), with at least one virus sampled each month between October 2016 and May 2017, and hence compatible with a limited outbreak of at least eight months duration. Together, these data suggest that while it is possible to identify potential localized transmission events, few developed into major EV outbreaks.

3.4 Characterization of EV-A virus source populations

To place viruses isolated in this study in the context of the global EV genetic diversity and investigate the potential source populations of those viruses, we conducted expansive phylogenetic analyses using the (full-length) VP1 gene (Fig. 4). These analyses highlight the overall genetic diversity of the viruses circulating in western Sydney and sampled within a single hospital. In the case of CV-A2, CV-A8 and EV-A71 (C2), a single virus introduction was identified, manifest as tight monophyletic clusters of sequences. There were also only a small number of samples for each of these viruses, and only one virus in the case of EV-A71. In contrast, multiple introductions of highly diverse CV-A4, A5, A6, CV-A10 and A16 viruses—at least six, two, nine, two and three events, respectively—were identified despite the relatively small data set. All these introductions seemingly resulted in only limited outbreaks. In contrast, one CV-A6 introduction resulted in a monophyletic group of eleven closely related viruses. As this cluster was sampled over the entire 18-month study period, it is indicative of persistent transmission within Sydney and hence may comprise a potentially endemic strain. Overall, the genetic diversity both within and between serotypes, and the lack of evidence of major disease outbreaks, suggest that the sustained transmission of EV-A within western Sydney is relatively uncommon.

Figure 4.

Phylogenetic trees of the VP1 coding gene of EV-A71 sub-genogroup C2, CV-A2, CV-A4, CV-A6, CV-A5, CV-A8, CV-A10 and CV-A16 sampled globally. All full-length VP1 gene sequences were downloaded from GenBank with duplicate sequences removed. Black squares represent prototype strains, prototype strains were not included in the EV-A71, CV-A10 and CV-A16 trees to improve phylogenetic resolution. Trees were mid-point rooted for clarity only. Branch color represents the continent from which samples were collected and circles depict sequenced obtained in this study. Asterisks (*) represent statistical support >70 per cent for major clades and those directly related to clades containing samples from this study estimated from 100 bootstrap replicates. Scale bars are proportional to the number of nucleotide substitutions per site.

It is difficult to determine source populations of many of the viruses isolated in CHW, in large part due to the substantial genetic diversity between the viruses sampled in this study and those previously sampled globally. This is especially true for serotypes that remain poorly sampled, particularly CV-A2, CV-A5 and CV-A8. One of the CV-A5 viruses isolated in this study had less than 85 per cent sequence similarity in the VP1 gene to their closest relatives, making it impossible to identify a source population and highlighting the diversity of globally circulating viruses that are still uncharacterized. With the more routinely sampled serotypes like CV-A6 and CV-A10 we were able to determine that each of the CHW viruses grouped with viruses sampled from Asian countries, including Japan, Thailand, India, Viet Nam and China (Fig. 4), suggesting that they are of Asian origin. Similarly, each of the three CV-A16 introductions grouped with viruses sampled in France, China and India, while the single EV-A71 strain clustered closely with viruses sampled in Poland at a similar time; however, we were unable to determine whether this strain was introduced into Australia from Poland, or vice versa, or if they were both derived from a common source population.

3.5 Phylogenetic correlates of disease severity

We performed a series of Bayesian phylogeny-trait association tests to determine if patient age or sex was correlated with serotype infection. As expected, this analysis revealed that there was no significant association between phylogenetic position (including virus serotype) and either the age or sex of the patient (P > 0.05 in all analyses). As the genetic diversity of EV-A is not well described by the VP1 gene alone, cluster associations were also explored for full-length genomes and the 3Dpol gene. Again, no cluster association was present (P > 0.05 in all analysis). Hence, on these data, there was no association between viral genetics and the age or the sex of the patient.

To determine whether there was any phylogenetic correlate of disease severity (i.e. the existence of distinct viral lineages of a particular severity stage), samples were grouped according to the severity of enteroviral disease using the schema described in Section 2. Of the viruses for which full-length genomes were obtained, disease phenotypes ranged from Stage 0A to Stage 2, that is for febrile illnesses with or without rash and meningitis without neurological symptoms. Again, Bayesian phylogeny-trait association tests of this sample of EV genomes revealed no significant association between disease severity among mild stages and phylogenetic position in the VP1, full-length genome or 3Dpol gene trees (P > 0.05), suggesting that specific disease outcomes can be caused by a range of EVs.

4. Discussion

A better understanding of the emergence and evolution of EV is required to determine the factors that impact the spread of these viruses through populations and from this improve strategies for outbreak mitigation and preparedness, including vaccine design. Although the majority of EV infections result in mild self-limiting febrile illnesses or HFMD, severe clinical manifestations have been associated with some EV serotypes, although the mechanisms of pathogenesis remain unclear.

Through an analysis of fifty-eight new full-length genomes we identified a remarkably diverse array of EV-A viruses collected within a single Sydney hospital over an 18-month sampling period. The sampled diversity comprised eight serotypes, two of which, CV-A4 and CV-A8, have not been previously characterized in Australia. Currently, only fifteen full-length EV-A genomes, ten non-EV-A71 and ninety-three full-length VP1 gene sequences, of which only fourteen are non-EV-A71, from Australian clinical samples are in the public domain (Clancy et al. 2005; Perera et al. 2007; Zhou et al. 2011). Between two and twenty EV-A-positive samples were collected each month throughout the study, with up to four different subtypes circulating in any 1 month.

The overall serotype diversity observed generally mirrored those characterized recently in a study of EV causing HFMD in Wenzhou, China (Guo et al. 2015). However, in contrast to that study, we found a lower proportion of EV-A cases caused by CV-A6 and EV-A71 and a higher proportion caused by CV-A2, A4, A5, A8, A10 and A16. While at a lower proportion than observed in China, CV-A6 was also the most frequently sampled serotype in this study representing 48 per cent of EV-A infections. This supports the large number of studies that describe viruses of the CV-A6 serotype as the current etiological agents of HFMD (Han et al. 2014; Hongyan et al. 2014; Li et al. 2014; Anh et al. 2018). Interestingly, viruses of the CV-A4 serotype were the second most characterized virus, constituting nearly 17 per cent of cases, although this serotype has not been detected previously in Australia. The diversity within the CV-A4 serotype and within other serotypes was evident in the 3Dpol phylogeny which revealed the presence of multiple sub-genogroups.

Phylogenetic analysis also revealed that the high level of genetic diversity present at the Children’s Hospital at Westmead was due to the introduction of multiple viruses into Australia each year, although it is notable that only rarely did these introductions lead to large disease clusters. Indeed, the majority of introductions resulted in no or little sustained transmission and these viruses quickly become extinct, as in the case of the two CV-A5 strains which were only detected once each, although this may in part reflect limited sampling. In contrast, a single CV-A6 strain was present in the population over the entire 18-month sampling period of the study, providing evidence of sustained transmission of what appears to be an endemic strain. The molecular epidemiology of EV-A in Australia can therefore be depicted as a pattern of recurrent introduction and continual extinction, combined with the persistent transmission of some endemic strains. More generally, these data suggest that novel EV-A strains are continually entering the Sydney community and have the potential to become endemic, although we are unable to predict which will achieve sustained transmission.

The identification of EV-A source locations may help in developing strategies for outbreak preparedness by improving early detection of viruses with greater epidemic potential. While a number of studies have attempted to identify source locations, it remains extremely difficult because of vagaries in sampling which is often both small and biased. Although the sample size of our study was relatively small, the genetic diversity of EV-A viruses at an epidemiological level sampled at a single hospital in Western Sydney mirrored the genetic diversity of viruses circulating globally. It was also clear that many of the EV-A viruses introduced into Australia group with viruses that originated from Asia and Europe. In other cases, however, such as CV-A2 and CV-A5, we were unable to identify source locations due to the low numbers of viruses isolated previously and the substantial genetic diversity between these viruses and those sampled here. Clearly, wider and more structured sampling is necessary to fully resolve the molecular epidemiology of human EV. Interestingly, the CV-A16 strains isolated in this study were most closely related to viruses isolate in Asia and Europe, while all CV-A6 strains were most closely related to Asian viruses. We therefore propose that viruses of the CV-A6 serotype are endemic in Asia, which seemingly acts as a transmission hub for this virus, yet with only relatively small numbers of infections occurring outside of this region. Hence, with respect to CV-A6, Asia constitutes a ‘source’ population that continues to re-seed the global ‘sink’ population. In contrast, CV-A16 appears to circulate on a more global level in a manner that shows some similarity to the migratory dynamics of influenza (Rambaut et al. 2008; Bahl et al. 2011), and as such viruses appear to have originated from a wider range of global localities including Europe.

High levels of historic intra- and inter-serotypic recombination were evident in the EV data collected here, such that EV-A viruses effectively constitute a single gene pool (Lukashev 2005; Simmonds and Welch 2006; Mirand et al. 2007; Zhang et al. 2010; Yip et al. 2013; Lukashev et al. 2014; Kyriakopoulou et al. 2015; Yen et al. 2017; Wang et al. 2018b ). The pattern of recombination in this study matched that observed in a recent study of viruses isolated in Wenzhou, China (Wang et al. 2018a). Specifically, we identified a novel CV-A4 recombinant, of which around 80 per cent of its genome was most closely related to viruses of the CV-A4 serotype but that carried a 3Dpol gene closely related to that present in CV-A10 viruses circulating in China in 2015 (Fig. 3, Supplementary Table S2 and Fig. S2). The second novel recombinant identified was of the CV-A6 serotype that appears to be endemic in Sydney. The recombination break point of this virus was identified within the 2C gene (Fig. 3, Supplementary Table S2 and Fig. S2), with the last 3 kb of the genome (which encodes non-structural proteins) most closely related to viruses classified as CV-A4, CV-A2 and EV-A71. The high level of recombination also highlights the complexity of using serotyping for virus classification. Indeed, because serotyping is based on VP1 alone it only represents around one-third of the genetic diversity of the genome, such that a significant amount of genetic information is lost in this classification process (Lukashev 2005). In particular, viruses of the same serotype can exhibit substantial genetic differences, while viruses of different serotypes can show high levels of sequence similarity in the non-structural proteins.

EV infections result in a very wide range of clinical outcomes, from mild self-limiting HFMD to severe disease outbreaks and meningo-encephalitis. Previous studies have associated severe disease phenotypes with viruses of the EV-A71 serotype, although attempts to identify virulence determinants either in vitro or through genetic characterization of clinical samples have been largely unconvincing (Li et al. 2011, 2016; Chang et al. 2012; Chen et al. 2012; Wen et al. 2013). We were also unable to associate disease phenotype with any particular phylogenetic clustering of sequences including serotype. This was true of whether the VP1, 3Dpol gene or full-length sequences were analyzed, and looking for any such association is clearly complicated by frequent recombination (so that, eg the VP1 mutations that define serotype can occur in a variety of genetic background). Unfortunately, the range of disease phenotypes from which we were able to obtain full-length genomes sequences was limited to mild disease so that the possible genetic basis of severe disease could not be investigated here. Larger studies with full-length sequences from patient with a wider range of disease phenotypes are therefore required to improve our ability to correlate disease severity with virus or host genetics.

EV-A71 is generally considered as the EV-A virus responsible for the largest number of outbreaks of severe neurological disease in young children globally. In this study, only a single EV-A71 infection resulted in Stage 2 disease, characterized by febrile illness and myolonic jerk. South-East Asia has often been suggested as the major source population for EVs to enter Australia due to its proximity and the number and scale of outbreaks identified (McMinn et al. 2001; Zander et al. 2014). The EV-A71 virus isolated in this study was of the C2 sub-genogroup, one of the rarer EV-A71 sub-genogroups that has previously been isolated in Australia (Brown et al. 1999; McMinn et al. 2001; Sanders et al. 2006; Zhou et al. 2010), but which is most commonly associated with disease in Europe (Hassel et al. 2015). In particular, four cases of aseptic meningitis caused by viruses of the C2 sub-genogroup were identified in Poland in mid-late 2016 (Wieczorek et al. 2018). These viruses exhibited high sequence similarity in the VP1 and 3Dpol genes (ranging between 98.2–99.0% and 97.4–98.1%, respectively) to the virus identified in this study suggesting that they share a recent common ancestor, although their precise origin was uncertain.

In sum, we have revealed a remarkably high level of EV-A genetic and serotypic diversity within a single Australian hospital, some of which has the potential to seed new epidemics. Further research on a larger sample is needed to better understand the source location of these viruses and to determine their mechanisms of pathogenesis to improve strategies for outbreak preparedness and vaccination.

Click here for additional data file.