ABSTRACT
Human enteroviruses of species A (EV-A) are the leading cause of hand-foot-and-mouth disease (HFMD). EV-A71 is frequently implicated in HFMD outbreaks and can also cause severe neurological manifestations. We investigated the molecular epidemiological processes at work and the contribution of genetic recombination to the evolutionary history of EV-A in Madagascar, focusing on the recently described EV-A71 genogroup F in particular. Twenty-three EV-A isolates, collected mostly in 2011 from healthy children living in various districts of Madagascar, were characterized by whole-genome sequencing. Eight different types were identified, highlighting the local circulation and diversity of EV-A. Comparative genome analysis revealed evidence of frequent recent intra- and intertypic genetic exchanges between the noncapsid sequences of Madagascan EV-A isolates. The three EV-A71 isolates had different evolutionary histories in terms of recombination, with one isolate displaying a mosaic genome resulting from recent genetic exchanges with Madagascan coxsackieviruses A7 and possibly A5 and A10 or common ancestors. The engineering and characterization of recombinants generated from progenitors belonging to different EV-A types or EV-A71 genogroups with distantly related nonstructural sequences indicated a high level of permissiveness for intertypic genetic exchange in EV-A. This permissiveness suggests that the primary viral functions associated with the nonstructural sequences have been highly conserved through the diversification and evolution of the EV-A species. No outbreak of disease due to EV-A has yet been reported in Madagascar, but the diversity, circulation, and evolution of these viruses justify surveillance of EV-A circulation and HFMD cases to prevent possible outbreaks due to emerging strains.
IMPORTANCE Human enteroviruses of species A (EV-A), including EV-A71, are the leading cause of hand-foot-and-mouth disease (HFMD) and may also cause severe neurological manifestations. We investigated the circulation and molecular evolution of EV-A in Madagascar, focusing particularly on the recently described EV-A71 genogroup F. Eight different types, collected mostly in 2011, were identified, highlighting the local circulation and diversity of EV-A. Comparative genome analysis revealed evidence of frequent genetic exchanges between the different types of isolates. The three EV-A71 isolates had different evolutionary histories in terms of recombination. The engineering and characterization of recombinants involving progenitors belonging to different EV-A types indicated a high degree of permissiveness for genetic exchange in EV-A. No outbreak of disease due to EV-A has yet been reported in Madagascar, but the diversity, circulation, and evolution of these viruses justify the surveillance of EV-A circulation to prevent possible HFMD outbreaks due to emerging strains.
INTRODUCTION
Human enteroviruses of species A (EV-A), which belong to the Picornaviridae family, are a leading cause of hand-foot-and-mouth disease (HFMD) (1, 2). HFMD is a self-limited febrile illness characterized by mouth ulcers, a spotty rash, and blisters on the palms and soles. Nineteen types of EV-A viruses are known to infect humans (e.g., coxsackievirus [CV] A4 to A7, A10, and A14 and EV-A120). The viruses most frequently implicated in HFMD outbreaks are CV-A6, CV-A16, and EV-A71 (3, 4). EV-A71 can also cause severe neurological manifestations, such as meningitis, encephalitis, meningoencephalitis, Guillain-Barré syndrome, and polio-like acute flaccid paralysis, particularly in children (5). The most severe cases may result in death due to cardiopulmonary failure. Major HFMD outbreaks associated with EV-A71 have regularly been reported in the Asia-Pacific region, in China, Malaysia, Taiwan, Vietnam, and Cambodia (6–8), making EV-A71 a major health concern in these countries. EV-A71 infections have been reported regularly in European countries, but outbreaks have been infrequent over the last 20 years (9–11).
EV-A viruses are nonenveloped viruses with a single-stranded (+) RNA genome. This genome consists of a single open reading frame (ORF) flanked by two untranslated regions (UTRs) at its 5′ and 3′ ends (here referred to as the 5' UTR and 3' UTR) (12). The ORF contains three coding regions (P1 to P3). P1 encodes the four viral capsid proteins (VP1 to VP4). P2 and P3 encode nonstructural proteins involved in intracellular viral replication. The error proneness of the viral RNA-dependent RNA polymerase drives the evolution of EV genomes through point mutations. In addition, recombination between EVs is also a major genetic mechanism in the emergence of viral lineages (13–16). Recombination frequently occurs between EV-B or EV-C viruses (16–21). Recombination is also involved in the evolution of EV-A viruses, including EV-A71 (22–27), and has been shown to be associated with the emergence of different (sub)lineages of EV-A71 (28).
Molecular epidemiology studies based on genomic sequence 1D encoding the VP1 protein (here referred to as 1DVP1) led to the classification of EV-A71 strains into phylogenetic lineages or genogroups A to G. Genogroups B and C have been implicated in most of the reported outbreaks worldwide over the last few decades. These two canonical genogroups have been further subdivided into subgenogroups B0 to B5 and C1 to C5. Phylogeographic studies of genogroups B and C revealed transmissions between Asia and Europe (29). Genogroup A includes the prototype strain (BrCr) isolated in 1969 and a few recent isolates from China (30). Genogroups D and G were recently reported in India (31). Genogroup E has been isolated exclusively in Africa (Central African Republic, Cameroon, and Senegal) (32–34), and genogroup F has been reported only in Madagascar (30). The EV-A71 genogroups D, E, F, and G were identified in patients with acute flaccid paralysis and healthy individuals but were not reported in any specific outbreaks. The rare identification of the EV-A71 genogroups D, E, F, and G raises important questions about their genetic origins, their links to the canonical genogroups B and C, and their pathogenicity.
We have published several reports on the diversity and evolution of cocirculating human enteroviruses of species C (EV-C) in Madagascar, including pathogenic circulating recombinant vaccine-derived polioviruses (cVDPVs) (35–37). The diversity and complexity of the recombinant EV-C forms, particularly for the cVDPVs that have been circulating for only a short period (2 to 3 years), suggested that recombination is frequent and that a given recombinant form has a short life span (38, 39). This may reflect the considerable diversity of EV-C types circulating in the population and the high frequency of intra- and intertypic recombination events (35, 36, 39). In addition, our previous studies suggested that new recombinant pathogenic cVDPV lineages emerge through the acquisition of sequences from EV-C, such as CV-A viruses (40, 41). We therefore investigated the occurrence of similar intertypic recombination events between EV-A strains circulating in Madagascar. In particular, we investigated the underlying molecular processes and the contribution of genetic recombination to the evolutionary history of EV-A in Madagascar by analyzing EV-A isolates circulating on the island over short time periods and by engineering recombinants, focusing specifically on EV-A71 from genogroup F.
RESULTS
Phylogenetic relationships between EV-A71 genogroup F and heterotypic EV-A strains circulating locally and worldwide.We explored the diversity of EV-A in Madagascar through the whole-genome sequencing of 23 EV-A isolates from different regions of the island. Most were isolated from the stools of healthy children in 2011 (30, 37). These viruses were typed on the basis of full 1DVP1 sequences (42), and eight different EV-A types were identified (Table 1).
Molecular typing of enterovirus A isolates collected in Madagascar
The complete genomes of Madagascan EV-A isolates were compared with a set of 810 publicly available genomes from 21 different EV-A types. We also included the full genome sequences of two representative EV-A71 isolates from genogroup E (Central African Republic and Cameroon) and one EV-A71 isolate from subgenogroup C1 from Cameroon. The initial data set was downsized to 153 representative sequences (see Table S1 in the supplemental material) for comprehensive phylogenetic analyses with the 5′ UTR, P1, P2, and P3 regions. The phylogenies reconstructed by the maximum likelihood (ML) method are shown in Fig. 1.
Maximum likelihood analysis of the nucleotide sequences of different genomic regions of the enterovirus A species. A sequence alignment including the complete genomes of 153 virus strains from 21 types of enterovirus A was constructed and split into four partition data sets corresponding to the capsid-encoding region P1 (A), the 5′ untranslated region (5′ UTR) (B), and the nonstructural protein-encoding regions P2 (C) and P3 (D). Phylograms were reconstructed by the maximum likelihood (ML) method. Bootstrap values of ≥70% are indicated at key nodes. For clarity, consistent clusters of strains are boxed with the name of the corresponding EV-A type or EV-A71 genogroup. Isolate naming is according to the following scheme: EV type, laboratory name, country code (e.g., MDG), and the last two numbers of the year of isolation. Madagascan sequences are highlighted in red. Types of isolates are indicated on the right (colored vertical lines and type names). The two main clusters of P2 and P3 sequences are indicated as clades I and II. Scale bars indicate the ML of substitution per nucleotide position. Gen, genogroup; Sgen, subgenogroup.
In the P1 capsid-encoding region, all sequences clustered to form a bootstrap (BST)-supported monophyletic cluster assigned to their respective type affiliations (BST value of 100%) with ≥76.2% nucleotide identity and ≥91.4% amino acid identity per type cluster (Fig. 1A). Within the EV-A71 type cluster, all EV-A71 sequences segregated into consistent subclusters (BST value of ≥99%) corresponding to the 1DVP1-based genogroups A, B, C, E, and F. Most EV-A isolates from Madagascar formed a distinct subcluster within their respective type cluster (BST value of ≥95%, ≥82.2% nucleotide identity, and ≥97.7% amino acid identity for each Madagascar subcluster).
The phylograms inferred from the 5′ UTR, P2, and P3 sequences displayed several bootstrap-supported topological incongruences (BST value of ≥90%) (Fig. 1B to D) relative to the phylogram inferred from P1 sequences (Fig. 1A). The clustering of EV-A71 genogroups A and E and of individual subgenogroups B and C (B0 to B5 and C1 to C5) was conserved in all trees, but incongruent clustering was observed for EV-A71 genogroup F and other Madagascan EV-A types. The 5′ UTR sequences of most Madagascan isolates segregated into clusters with sequences of other isolates of the same EV type (Fig. 1B), some of which were related to sequences from isolates collected in China. On the trees inferred from the P2 (Fig. 1C) and P3 (Fig. 1D) sequences, the EV-A71 isolates of subgenogroups B and C were scattered over different, distant parts of the tree and were clustered with other EV-A viruses. These patterns are consistent with the previously described recombinant forms of genogroups B and C (22, 24, 26, 28, 43). All but two of the EV-A isolates from Madagascar clustered in a large clade (Fig. 1C and D, clade I), including EV-A71 genogroups E and F and CV-A5 and CV-A10 sequences from China recovered in 2009 and 2013 (22). Only two CV-A10 isolates from Madagascar (MAD3995 and MAD9506) had P2 and P3 sequences distantly related to those of other Madagascan isolates. They clustered in a second large clade (clade II) including EV-A71 sequences from subgenogroups B3 and C4 (BST values of 96% and 95% for the P2 and P3 phylograms, respectively).
Overall, our data showing that the P1 structural genomic regions of most EV-A isolates from Madagascar constitute distinct subclusters within their type clusters suggest that these isolates have been circulating and evolving locally in the island population for a relatively long period of time. In addition, the clustering of the sequences of most of the P2 and P3 nonstructural genomic regions of EV-A isolates from Madagascar, including those of EV-A71 genogroup F, strongly suggests that a local long-term evolutionary process involving intertypic recombination events has been at work.
Recombination in Madagascan virus isolates and the evolution of EV-A71 genogroup F.Recombination was investigated by comparing the full-length genomic sequences of all Madagascan EV-A isolates, EV-A71 genogroup E, the EV-A71 prototype strain BrCr, and other isolates displaying evidence of genetic relationships in the phylogenetic trees described above. Similarity analyses showed that the Madagascan EV-A isolates had signatures of genetic exchange in the P2 and P3 regions, as indicated by genomic fragments with nucleotide identify of ≥95% (Fig. 2A and Fig. S1).
Genomic structure of mosaic recombinant genomes of enterovirus A isolates. (A) Similarity plot analyses of different EV-A genomic sequences, with the genomic sequence of EV-A71 MAD7842-MDG10 as a reference. A scheme indicating the genomic organization of EV and the structural and nonstructural regions is shown. Similarity plot profiles are color coded according to EV type or the genogroup of EV-A71. Numbers of isolates for each type and genogroup are indicated in brackets. The nucleotide sequence identity threshold indicating evolution through genetic exchanges in the noncoding and nonstructural coding regions was set at 95%. Additional similarity plots with this and other EV-A genomes as references are shown Fig. S1 in the supplemental material. (B) Putative genomic crossovers with recombinant partners and positions in the nonstructural protein-coding regions of a set of EV-A genomes. The consistency of recombination events was assessed with seven different statistical methods implemented in a recombination detection program (RDP4) (44). Genomes for which significant evidence of recombination (P value of <0.001) was obtained with 7 methods (****), 6 methods (***), 5 methods (**), and 4 methods (*) are indicated. The nucleotide positions of the recombinant parts are indicated by dots if determined by the program (99% confidence interval) or by arrows if undetermined. Detailed values obtained with the program are presented in Table S2.
Genetic exchanges and their boundaries were considered validated if supported by at least four of the seven statistical methods used in the recombination detection program RDP4 (44). Intratypic exchanges were observed for CV-A7 isolates, and multiple intertypic exchanges between CV-A4, -A5, -A6, -A7, and -A10, EV-A71, and EV-A120 isolates were detected (Fig. 2B and Table 2 for statistical data).
RDP analysis of recombination events in enterovirus A genomes
The genome of the EV-A71 genogroup F isolate MAD7842 displayed significant evidence of evolution through recombination with genomic fragments of CV-A5 isolates (MAD2940 or MAD3144) in the 2C protein-encoding region and the 3′ part of the 3D region (Fig. 2 and Table 2). A large part of the genome of this isolate, encoding proteins 2B to 3D, presented significant evidence of recombination with the genomes of CV-A10 (MAD9856 and MAD9937) and CV-A7 (MAD9750) isolates. We cannot exclude the possibility that these recombination events involved unknown recent common ancestors.
The high degree of similarity between the nucleotide sequences of defined genomic fragments encoding nonstructural proteins of certain cocirculating Madagascan EV-A isolates provides strong evidence for the occurrence of recent and frequent intra- and intertypic genetic exchanges leading to mosaic recombinant genomes.
Comparative analysis of the evolutionary history and temporal origin of EV-A71 genogroup F for two different genomic regions.The evolutionary history of the capsid-encoding region of EV-A71 was investigated by performing Bayesian analysis with an EV-A71 1DVP1 sequence data set (Fig. 3A and B and Table S1). The mean evolution rate was 3.97 × 10−3 (95% highest posterior density [HPD], 3.47 × 10−3 to 4.45 × 10−3) substitutions per site year. The most recent common ancestor of all EV-A71 lineages was dated to 1901 (95% HPD, 1870 to 1931). The three EV-A71 sequences sampled in Madagascar between 2004 and 2011 (genogroup F) shared a common ancestor that arose in 1995 (node N10; 95% HPD, 1990 to 2000; posterior probability [PP] = 1) (Fig. 3B). Interestingly, the isolates of genogroup F and those of the cosmopolitan genogroup B diverged from a common ancestor in 1946 (node N2; 95% HPD, 1937 to 1955; PP = 0.98). The isolates of the other newly described EV-A71 genogroups D and G (reported in India) and genogroup E (reported in sub-Saharan Africa) diverged at a similar time from another common ancestor (node N3; PP = 0.95) in 1950 (95% HPD, 1935 to 1964).
Evolutionary history of EV-A71 genogroup F, as assessed by maximum clade credibility phylogenetic analysis. (A) Patterns of evolution were reconstructed for EV-A71 from the nucleotide sequences encoding the capsid protein 1DVP1 and for a segment of the P3 (nonstructural protein-encoding) region, as indicated on the schematic diagram of the EV genome. (B) Maximum clade credibility tree based on the capsid protein 1DVP1 sequence, including EV-A71 of various geographic origins from all described EV-A71 (sub)genogroups. (C) Maximum clade credibility tree based on the nonstructural P3 region specifically including EV-A71 strains within genogroup F and EV-A sequences belonging to clade I (Fig. 1D). The maximum clade credibility trees were inferred by Bayesian Markov chain Monte Carlo analysis, with a discrete model, to determine the virus corresponding to tree branches. Tips are scaled according to the date of virus sampling (timescale on the x axis). The main lineages are color coded according to EV-A71 (sub)genogroup and EV-A type. Time to the most recent common ancestor (tMRCA) and highest posterior probability density (HPD) intervals (blue lines) are indicated for the main nodes. The size of the circles at the nodes is proportional to posterior probability (PP).
We investigated the evolutionary history of lineages genetically related to recombinant EV-A71 genogroup F (clade I in the above ML trees) (Fig. 1C and D) by Bayesian analysis of a segment of the P3 region (nucleotide positions 5200 to 6200). We found no evidence of a recombination breakpoint common to different recombinant lineages within this segment (Fig. 2 and data not shown). The mean evolution rate was 6.54 × 10−3 (95% HPD, 4.30 × 10−3 to 8.90 × 10−3) substitutions per site year. The three EV-A71 genogroup F sequences had different evolutionary histories (Fig. 3C). The P3 segment of the EV-A71 isolate MAD7842 had a common direct ancestry (PP = 1) with segments of six other heterotypic virus strains sampled in Madagascar and was assigned to EV-A types CV-A7 and CV-A10. The common ancestor of this cluster was dated to 2007 (95% HPD, 2005 to 2009) and was assigned to CV-A7 with a probability of 0.84. The other putative EV types were EV-A71 (P = 0.15) and CV-A10 (P = 0.01). This phylogenetic pattern suggested that the EV-A71 genogroup F isolate MAD7842 acquired all or part of its P3 genomic region through recombination with a CV-A7 strain. The other two EV-A71 genogroup F sequences had different unresolved evolutionary patterns and were temporally distant from the EV-A71 isolate MAD7842 cluster. According to P3 sequences, the origin of the genogroup F lineage was dated to 1994 (95% HPD, 1988 to 2001), consistent with the 1DVP1-based phylogeny (1995; 95% HPD, 1990 to 2000) (compare Fig. 3B and C).
Overall, Bayesian phylogenetic analysis indicated a recent divergence of the Madagascan EV-A71 genogroup F isolates (common ancestor dating back to 1995) and a genetic link between this genogroup and the cosmopolitan genogroup B (common ancestor dating back to 1946). The nonstructural P3 sequences of the Madagascan EV-A71 isolates followed a different evolutionary history. This pattern confirmed the existence of active evolution through frequent intertypic recombination events in Madagascar, with a recent event dating back to about 3 years before isolation.
Conservation of functional compatibility among genetic domains of EV-A71 and other EV-A viruses.We investigated the propensity for recombination of EV-A71 and other EV-A viruses by constructing recombinant genomes with 5′ and 3′ halves originating from different EV types. We combined the 5' UTR-P1-2A genomic region of a given EV-A (strain X) with the 2BC-P3-3′ UTR region of another isolate (strain Y) from the same or another phylogenetic clade. We selected the progenitors from clade I (EV-A71 genogroup F, isolate MAD72341; CV-A7, isolate MAD9750) and clade II (EV-A71 subgenogroup C4, isolate SEP06; CV-A10, isolate MAD3995) (Fig. 4). A recombinant of the X and Y strains (X/Y) and its converse counterpart Y/X were constructed by a fusion-PCR-based procedure developed for the engineering of chimeric genomes (45). The transcription of the viral chimeric cDNA was controlled with the T7 RNA polymerase promoter, and genomic RNAs and viruses were produced in human 293T/T7 cells expressing this polymerase following the transfection of viral cDNAs. We engineered recombination breakpoints at the junction of viral proteins 2A and 2B because this genomic stretch is subject to frequent recombination in natural isolates (Fig. 2A and Fig. S1). We produced the following chimeric X/Y genomes: EV-A71(C4)/CV-A10, EV-A71(F)/CV-A10, EV-A71(C4)/CV-A7, and EV-A71(F)/CV-A7; their converse YX genomes were also produced. We also used the full-length genomes of the four primary progenitors to transfect cells to obtain infectious virus controls.
Comparative plaque phenotypes of the progenitor and chimeric viruses. Plaque phenotype (A to D) was assessed for the progenitor and chimeric viruses on Vero E6 cells or RD cells for viruses with the 5′ half of the CV-A10 genome. Negative and positive controls (coxsackievirus B3) are presented (E). Schematic diagrams of the progenitor and chimeric genomes are shown. Virus names and genomes are color coded according to genetic background and affiliation to the nonstructural CDS clade I (red) or clade II (blue) (Fig. 1C and D). Titers of viral suspensions are indicated.
All of the generated chimeric genomes yielded viable recombinant viruses after the transfection of 293T/T7 cells with viral cDNA and one passage in sensitive cell lines, as indicated by cytopathic effects (CPEs) in infected cells. The nucleotide sequences of the recombinant genomes were checked by Illumina high-throughput sequencing. The recombinant genomes displayed 0- to 3-nucleotide (nt) differences relative to the progenitor sequences (Table 3). These differences were located mostly in the VP1 coding sequence. We hypothesized that these mutations were introduced during gene amplification or selected from the genome pool present in the original stock of progenitor viruses. We cannot rule out the possibility that some of these substitutions were selected to adapt viral functions to a given recombinant.
Nucleotide and amino acid substitutions in the reconstituted progenitors and engineered chimeric viruses
The phenotypic characteristics of progenitors and chimeric viruses were assessed through plaque assays and replication kinetics in Vero E6 cells. Viruses with the 5′ half of CV-A10 genomes were analyzed in rhabdomyosarcoma (RD) cells, the only cells in which they replicate efficiently. Most of the progenitors and recombinant viruses studied formed turbid plaques of various sizes (Fig. 4A to D), contrasting with the round plaques generally produced by coxsackievirus B3 (EV-B) (Fig. 4E). Titers of viral stocks were assessed by determining the 50% tissue culture infective dose (TCID50) per milliliter and the number of PFU per milliliter. The differences between these two values for a given stock varied between viruses, suggesting differences in the efficiency with which clear-cut CPEs were induced in infected cells. The plaques of some of the recombinant viruses were similar in size and appearance to those of the donor of the 5′ half of the genome (Fig. 4D, CV-A7/X). In other recombinants, the presence of a certain 3′ half of the genome modified the plaque phenotype, e.g., EV-A71(C4)/CV-A7 versus A71(C4)/SEP06, and EV-A71(F)/CV-A7 versus EV-A71(F)/MAD7234 in Fig. 4A and B, respectively.
We then compared the growth kinetics of the progenitor and recombinant viruses (Fig. 5A to D). Progenitor viruses and most of the chimeric viruses induced complete cell lysis within 24 to 48 h of infection. All but one [EV-A71(C4)/CV-A7] of the chimeric viruses used had growth kinetics similar to the kinetics of the progenitors, with only minor differences in virus yields and genome copy numbers. The EV-A71(C4)/CV-A7 chimeric virus induced only partial CPE, consistent with its tiny plaques (Fig. 5A) and its titer and genome copy number, which differed from those of its EV-A71 progenitor.
Replication kinetics of the progenitor and chimeric viruses. Replication kinetics of the progenitor (black) and chimeric (gray and white) viruses were analyzed in triplicate following infection of Vero E6 cells or RD cells. RD cells were used for viruses carrying the 5′ half of the CV-A10 genome. Infected cells and supernatants were collected at 0, 6, 16, 24, and 48 h postinfection (hpi). Virus production is expressed as the number of genome copies and TCID50, as indicated. The data are presented are mean values ± standard deviations. The occurrence of a cytopathic effect and complete cell lysis are indicated. Virus names and genomes are color coded according to genetic background and affiliation to the nonstructural CDS clade I (red) or clade II (blue), as indicated in the legend of Fig. 4.
Overall, these results indicate a high degree of conservation and compatibility of viral functions through genetic exchanges between EV-A71 genogroups C4 and F and CV-A7 and CV-A10, allowing intertypic recombination. In certain cases, recombination can lead to phenotypic modifications.
DISCUSSION
Characterization of a set of EV-A isolates from Madagascar demonstrated the diversity of viral types cocirculating throughout this country and provided evidence of local evolutionary processes characterized by frequent recent intertypic genetic exchanges leading to mosaic recombinant genomes. EV-A71 isolates from the recently discovered genogroup F were subject to a similar active evolution process, strongly suggesting that EV-A strains have been circulating and interacting genetically with each other for many years on this island. The construction and characterization of viable recombinants suggested that EV-A evolution has kept primary viral functions compatible for intertypic recombination.
At least eight of the known EV-A types were cocirculating all over the island (Fig. 6) during the short period considered here (Table 1), highlighting the diversity of the EV-A strains circulating in the population. Similar levels of EV-A diversity have been reported in Asia and Europe for clinical specimens, stools from refugees, and wastewater (1, 22, 46–49). Such high levels of diversity have already been reported in Madagascar for EV species B and C (35, 36, 50). The phylogenetic analysis of the capsid-encoding region P1, based on 810 available published EV-A genomes, revealed that all reported EV-A isolates from Madagascar belonged to distinct lineages different from the other EV-A isolates circulating worldwide. Furthermore, the phylogenetic links between their nonstructural genomic regions suggested that most of the Madagascan isolates shared a recent common ancestor, resulting in a robust clade of related sequences (clade I) (Fig. 1C and D). The phylogenetic analysis of whole virus genomes from different Madagascan districts provided further evidence of the extensive circulation of EV-A isolates throughout the country, constituting a dynamic local EV-A ecosystem.
Circulation of enterovirus A isolates in different geographic areas of Madagascar. The Madagascan isolates used in this study are color coded by type and shown on the map of the island (map from Wikimedia Commons Atlas of the World [https://commons.wikimedia.org/wiki/Atlas_of_the_world]).
However, despite the insularity of the EV-A ecosystem in Madagascar, it displays relationships with other EV-A strains from other countries. Our Bayesian phylogenies indicated a common ancestry for the genes encoding the VP1 capsid proteins of EV-A71 genogroups F and B in the 1940s even though EV-A71 genogroup F has only ever been reported in Madagascar. The nonstructural regions of a number of EV-A isolates recovered from Madagascar are related to those of Chinese CV-A5 and -A10 strains (22) and the recently reported African EV-A71 genogroup E (32–34). In addition, the two CV-A10 isolates in clade II (MAD3995 and MAD9506) were collected in the Mahajanga district and were distantly related to the other virus isolates collected in Madagascar. The Mahajanga district is a port of entry and destination for people originating from the Comoros Islands, Europe, India, and, more recently, China, so these CV-A10 isolates may belong to a recently imported lineage. Virus importations may, therefore, contribute to the diversity of the local EV-A ecosystem. Only a few EV-A sequences from Africa are available in databases. Given the proximity of this continent to Madagascar, we cannot rule out the possibility that the discovery of new isolates in Africa in the future will modify our view of the insular nature and molecular evolution of isolates from Madagascar.
The phylogenetic incongruences between structural and nonstructural genes suggest that multiple intertypic recombination events have occurred within the Madagascan EV-A ecosystem. Recombination has already been documented in the evolutionary history of EV-A, and the recombinant forms described here are consistent with those described in previous studies (22, 24, 28, 49). Most EV-A strains from Madagascar appear to have mosaic genomes combining fragments with nucleotide sequences displaying a high degree of identity (>95%) to two to four different viral isolates. EV-A71 isolates from genogroup F have been subject to a similar active evolution process because the three EV-A71 strains studied had different patterns of recombination. In particular, the EV-A71 MAD7842 isolate had a striking mosaic genome that emerged following several recombination events with partners related to the cocirculating CV-A7 and, possibly, CV-A5 and CV-A10 viruses. Previous studies in Madagascar reported similar genome plasticity and mosaicism within EV-C species (38) and provided evidence of local evolutionary processes, with frequent recent intertypic genetic exchanges leading to mosaic recombinant genomes. Present data are consistent with a model of modular evolution in which the diversification of the capsid genes determining EV type and genogroup is dissociated from that of the nonstructural genes (21, 41, 51–55). Previous studies have highlighted the occurrence of much lower rates of recombination (by a factor of about five) in EV-A71 than in other species or types, such as CV-A2, -A4, or -A10, and suggested a high degree of conservation of EV-A71 among enterovirus types and species (16, 24, 26, 28, 43). However, this study and other recent studies analyzing the genetic features of EV-A71 and other EV-A strains circulating in Africa, Taiwan, and elsewhere suggested that EV-A71 and cocirculating EV-A strains may evolve more frequently through intertypic recombination in certain contexts (49, 56).
We evaluated these possible intertypic genetic exchanges and the genetic features limiting recombination further by engineering eight different intertypic chimeric viruses from four EV progenitors. These progenitors belonged to three different types and two different EV-A71 genogroups (F and C4) and had closely or distantly related nonstructural sequences (Fig. 1, clades I and II). All chimeric combinations and counterparts gave rise to viable infectious viruses, highlighting the tremendous genome plasticity of EV-A. The primary viral functions of most of the chimeric viruses were similar to those of their progenitors. These findings merit confirmation through the engineering of other kinds of EV-A recombinants, but the implications of our findings are already far-reaching. Indeed, they suggest that despite long-term diversification, the primary viral functions associated with the structural and nonstructural genomic regions of EV-A species have been conserved and remain compatible for intra- and intertypic recombination. Both characteristics, conservation of viral functions and compatibility for recombination, may have remained closely linked during evolution, as proposed for EV species (55).
Only two of the engineered chimeras, EV-A71(F)/CV-A7 and EV-A71(C4)/CV-A7, gave rise to viruses causing weaker cytopathic effects and smaller plaques than their progenitors. These data indicate that, in certain cases, intertypic recombination can lead to phenotypic modifications. By engineering intratypic recombinants of EV-A71 subgenogroups B3 and C2, Phuektes et al. showed that the capsid-encoding region (P1) has the strongest effect on phenotype and that the P2/P3 regions have a more moderate effect (57). These phenotypic differences may be associated with changes in the pathogenicity of recombinants in infected hosts. However, we cannot exclude the possibility that recombinants displaying no differences relative to their progenitors in cultured cells have different virulence characteristics in infected hosts. Recombination has been shown to modify EV-C pathogenicity in animals without affecting the capacity of the virus to replicate in cultured cells (40, 41, 58).
The EV-A types detected in Madagascar have been reported in diverse countries in sporadic cases and outbreaks of disease, including febrile illness, herpangina, HFMD, and more severe diseases, such as encephalitis and polio-like syndrome (59, 60). No outbreak linked to EV-A has yet been reported on the island. However, phylogenetic analysis revealed extensive circulation of EV-A in the island population. In particular, our findings provide evidence for the evolutionary diversification of EV-A71 genogroup F in Madagascar (Fig. 1, 2, and 6) due to mutation and frequent genetic exchanges with other EV types. Virus migration also plays an important role, as shown by the recent importation of CV-A10 and the phylogenetic links between EV-A isolates from Madagascar and virus isolates from other countries. The phenotypes of the engineered recombinants indicate that local isolates can evolve through intertypic genetic exchanges with distantly related EV-A. All these data suggest that surveillance of human EV-A and HFMD cases in Madagascar should be implemented to detect possible HFMD and encephalitis outbreaks due to emerging strains.
MATERIALS AND METHODS
Ethics statement.Field missions were performed in 2011, with the collection of stool specimens from healthy children under the age of 5 years. The study protocol was reviewed and approved by the National Research Ethics Committee (Madagascar Ethics Committee agreements 011-MSANP/CE and 046-MSANP/CE). Written informed consent was obtained from the parents or legal guardians of the children enrolled in the study.
Isolation of enterovirus A in Madagascar.EV-A was isolated on human rhabdomyosarcoma (RD) cells from stool specimens treated with chloroform, as described in the standard procedures of the laboratory manual for the WHO Global Polio Laboratory Network (61). Stool specimens were sampled in 2011 from healthy children or, in a few cases, from children with acute flaccid paralysis. Isolates with cytopathic effects on RD cells were primarily typed by molecular methods through partial 1DVP1 sequencing, as previously described (50, 62).
Whole-genome sequencing.Whole-genome sequences were determined for EV-A from Madagascar and for previously reported EV-A71 isolates from genogroup F (MAD72341 and MAD3126), genogroup E (CAF-NMA-03-008 and CAE-C08-146), and subgenogroup C1 (CAE-C08-041) (30, 32, 34). RNA genomes were extracted with a High Pure Virus RNA kit (Roche), according to the manufacturer’s instructions. Complete genome sequences were determined by genome walking. The terminal sequences of the viral genomes were confirmed by rapid amplification of the cDNA ends (RACE) with a 5′/3′ RACE kit (Roche). Sanger sequencing was performed with a BigDye Terminator kit on an ABI 3730xl 96-capillary DNA Analyser (Applied Biosystems). Alternatively, viral genome sequences were obtained by the Illumina high-throughput technique with generic primers specific for EV-A (63).
Data collection and nucleotide sequence data sets.Complete genome sequences for all enteroviruses from species A were retrieved from the GenBank database (to March 2016). Wrongly annotated sequences, partial genomes, coding sequences (CDS) including internal termination codons, and sequences without a fully specified date (year) and country of origin were discarded. A sequence alignment including the complete genomes of 810 virus strains from 21 enterovirus types was generated with ClustalW implemented in BioEdit, version 7.2.3, software (http://www.mbio.ncsu.edu/bioedit/bioedit.html) and split into four partitions corresponding to the 5′ UTR, P1, P2, and P3 genomic regions. For more accurate investigation of the phylogenetic relationships by maximum likelihood (ML) analysis, the initial global data sets were downsized to 153 sequences representing genetic clusters and various outbreaks at different time periods and locations (see Table S1 in the supplemental material).
Maximum likelihood and neighbor-joining phylogenetic analysis.Phylograms based on the sequence partitions corresponding to the 5′ UTR and P1, P2, and P3 regions were reconstructed with different methods. Maximum likelihood (ML) trees were generated with PhyML software (64), with the general time-reversible plus invariant site and gamma distribution model (GTR plus I and γ) and 500 bootstrap replicates for the estimation of node consistency. The ML best-fit model for nucleotide substitution was validated statistically with JModel Test software (65, 66). The resulting ML trees were annotated with FigTree, version 1.4.2, software (http://tree.bio.ed.ac.uk/). Neighbor-joining (NJ) trees were generated with MEGA6 software (67), using the Kimura two-parameter evolutionary model and 1,000 bootstrap replicates for the estimation of node consistency.
Recombination analysis.The degree of nucleotide sequence identity between complete genome sequences was analyzed with SimPlot, version 3.5.1, software (68). The pairwise similarity between all sequences in a multiple sequence alignment was calculated with a 200-nt window moved along the sequence in 20-nt steps. The consistency of recombination events was evaluated with seven different statistical methods implemented in RDP4 software (44). In some cases, we checked that the recombinant genome sequences did not result from fortuitous cross-sequencing of viral mixtures present in isolates by carrying out control Illumina high-throughput sequencing following reverse transcription-PCR (RT-PCR) with generic primers specific for EV-A (63).
Sequence data for Bayesian phylogenetic analyses and tree reconstruction.For the inference of coalescent times for the EV-A71 isolates from genogroups A to G, a sample of 190 1DVP1 sequences, including the three EV-A71 sequences (MAD7842, MAD3126, and MAD72341), was randomly selected from the sets of sequences used in our earlier analyses (29, 30). We subsampled sequences of strains recovered between 1963 and 2015 from the seven EV-A71 genogroups and all subgenogroups of genogroups B and C to obtain this 1DVP1 sequence data set. We demonstrated that including more sequences in the data set had no effect on the coalescent times calculated from the genealogy.
The sequence samples are described in Table S1. The EV-A71 genealogy was inferred with the Bayesian Markov chain Monte Carlo (MCMC) methods available in BEAST, version 1.8.2, software (69). Phylogenetic inference was performed with the general time reversible (GTR) substitution model with gamma rate heterogeneity across sites. Coalescent times were inferred under a relaxed molecular clock model (uncorrelated log-normal distribution across lineages). A Bayesian Skyline prior was chosen as the tree prior, with 10 different population size intervals and a piecewise linear model (70). The tree was reconstructed with a Markov chain Monte Carlo process of 120 million states, sampled every 12,000 states. A second sequence data set (n = 30 taxa) was selected from the ML phylogeny inferred from the P3 sequences. This sample included 27 taxa genetically related to the three EV-A71 genogroup F isolates in a monophyletic cluster. The phylogenetic analysis was performed as described above, except that a discrete trait analysis with a Bayesian stochastic search variable selection (BSSVS) procedure was included to investigate the changes in EV types across lineages. The MCMC process was performed with 50 million states and a sampling frequency of 1/5,000.
Cell culture.Monkey kidney epithelial Vero E6 cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM) with GlutaMAX (Gibco 31966-021) and supplemented with 5% fetal calf serum. Human rhabdomyosarcoma RD cells (WHO) were cultured in DMEM (T15-106; PAA) supplemented with 2 mM l-glutamine and 5% fetal calf serum. Human embryonic kidney 293T/T7 cells (Pierre Charneau, Institut Pasteur, Paris, France), which constitutively express the T7 RNA polymerase, were cultured in minimal essential medium (MEM) (M5650; Sigma) supplemented with 2 mM l-glutamine and 10% fetal calf serum. All cell lines were cultured at 37°C under a humid atmosphere containing 5% CO2.
Engineering of chimeric genomes.Chimeric genomes were engineered with a previously described fusion-PCR-based procedure (45). The RNA genomes of the progenitor viruses were denatured (5 min at 65°C). Then they were subjected to full-length reverse transcription (30 min at 50°C) with a Maxima H Minus First Strand cDNA synthesis kit (Thermo Scientific) with specific reverse primers (20 pmol) containing a 5′ poly(T) tail (n = 24) followed by the last 24 complementary nucleotides of the considered genome (a full list of primers is given in Table S2). Then full-length cDNA genomes were amplified using a Phusion High-Fidelity DNA polymerase (Thermo scientific) with the previous reverse primers and a forward primer containing a T7 promoter and the first 25 nt of the virus genome (25 pmol each). The amplification conditions were as follows: 30 s at 98°C, a preamplification step (7 s at 98°C, 20 s at 63°C, and 3.5 min at 72°C), and 20 cycles of 7 s at 98°C, 20 s at 65°C, and 3.5 min at 72°C. After a final elongation for 5 min at 72°C, the progenitor T7 promoter/cDNAs were purified with a QIAquick PCR purification kit (Qiagen). The full-length T7 promoter/cDNAs of the progenitor viruses were used both to regenerate progenitor viruses through cell transfection and as a template to produce overlapping segments for chimeric genome construction. Overlapping segments were produced by nested PCR with a Phire Hot Start Pol II kit (Thermo Scientific) and overlapping primers (10 pmol each) containing 50 nt of nonpriming 5′ tails designed to match the sequence of one progenitor genome upstream or downstream from the 2A/2B breakpoint, followed by 25 nt matching the second progenitor genome. The amplification conditions were as follows: 30 s at 98°C, followed by 30 cycles of 5 s at 98°C, 5 s at 60°C, and 40 s at 72°C. After a final elongation phase for 3 min at 72°C, the amplified overlapping segments were subjected to electrophoresis in an agarose gel, purified with a QIAquick gel purification kit (Qiagen), and quantified with a NanoDrop spectrophotometer (Thermo Scientific).
Equimolar amounts of purified overlapping amplicons were mixed (10-μl final volume) and ligated with Phire Hot Start Pol II (Thermo Scientific) in a final volume of 50 μl without added primers. The fusion conditions were 30 s at 98°C followed by 10 cycles of 5 s at 98°C, 1 min at 56°C, and 3 min at 72°C. The forward and reverse primers previously used for full-length genome amplification (20 pmol each) were then added to the fusion mixture, and PCR resumed, with heating for 30 s at 98°C followed by 20 cycles of 5 s at 98°C, 5 s at 65°C, and 3 min at 72°C. After a final elongation for 7 min at 72°C, the T7 promoter/chimeric genome amplicons were purified with a QIAquick PCR Purification kit (Qiagen) and quantified with a NanoDrop spectrophotometer (Thermo Scientific).
Recovery of chimeric viruses.Human embryonic kidney 293T/T7 cells (Pierre Charneau, Institut Pasteur, Paris, France), which constitutively express the T7 RNA polymerase, were transfected with 1 μg of purified T7 promoter/cDNAs from progenitor or chimeric genomes with FuGENE HD transfection reagent (Promega). Transfected cells were incubated at 37°C until complete cell lysis. Cell supernatants were collected and clarified by centrifugation, and 200 μl of the resulting suspension was passaged once on 80% confluent Vero E6 or RD cells grown in T-25 flasks incubated at 37°C. When Vero E6 or RD cell lysis was complete, flasks were frozen and thawed three times. Cell supernatants were collected, clarified by centrifugation, and stored at −20°C until use. Virus titers were determined in microplates by calculating the 50% tissue culture infective dose (TCID50) per milliliter according to the Kärber method. All recovered progeny virus genomes were sequenced with Illumina high-throughput sequencing technology.
Viral plaque assay.Vero E6 or RD cell monolayers in six-well culture plates (2 × 106 per well) were washed twice and infected with 500 μl of 10-fold serial dilutions of viral stock suspension per well. Viral adsorption was allowed to proceed for 2 h, and the cells were then washed, covered with 2.4% (vol/vol) Avicel RC-581 (FMC BioPolymers, Brussels, Belgium) and 2× culture medium (2× DMEM supplemented with 40 mM Tricine, 2.2 g/liter NaHCO3, 50 mM MgCl2, 4 mM l-glutamine, 4% fetal bovine serum [FBS]) and incubated at 37°C under a humidified atmosphere containing 5% CO2. After a week, cells were fixed by incubation for 1 h with one volume of 4% paraformaldehyde. The overlays were removed from the fixed cells, which were washed with phosphate-buffered saline (PBS) and stained with crystal violet. Coxsackievirus B3 (CV-B3), an EV-B virus, was used as a positive control for plaque assays. CV-B3-infected cells were fixed and stained 3 days postinfection.
Growth kinetics.Growth kinetics assays were performed in 96-well culture plates seeded overnight with 1 × 104 Vero E6 or RD cells per well. The cells were washed and infected in three independent replicates at a multiplicity of infection of 1. Virus adsorption was allowed to proceed for 2 h. The free viral inoculum was then removed, and the cells were washed twice and incubated at 37°C. Infected cells were frozen at different times postinfection. Plates were subjected to three freeze-thaw cycles, and each supernatant was then collected and clarified by centrifugation. Virus titers were determined in microplates by calculating the 50% tissue culture infective dose (TCID50) according to the Kärber method. The number of viral genome copies from each supernatant was evaluated by reverse transcription-quantitative PCR (RT-qPCR), as previously described (71).
ACKNOWLEDGMENTS
We thank Pierre Charneau for providing the 293T/T7 cells and Julie Sappa for help with preparing the English manuscript.
This work was supported by the Institut Pasteur, a Program Transversal de Recherche PTR484 grant, and an Action Concertée des Instituts Pasteur ACIP A22-16 and Roux Howard Cantarini postdoctoral fellowship (R.V.) (https://www.pasteur.fr/fr) and by the Fondation Total grant S-CM15010-05B (https://www.foundation.total/fr).
The funders had no role in study design, data collection and interpretation, or in the decision to submit the work for publication.
R.V., R.R., J.-M.H., and F.D. initiated and designed this study; R.V., R.R., M.-L.J., M.B., S.R., S.A., J.R., B.B., and J.-L.B. collected data and performed experiments; R.V., J.-L.B., and F.D. wrote the first draft of the manuscript. All authors revised the manuscript.
FOOTNOTES
- Received 23 September 2018.
- Accepted 18 December 2018.
- Accepted manuscript posted online 2 January 2019.
Supplemental material for this article may be found at https://doi.org/10.1128/JVI.01667-18.
- Copyright © 2019 American Society for Microbiology.
