Introduction

Hand, foot, and mouth disease (HFMD) is a common infectious disease in young children that often caused by Enterovirus A species (EV-A) (Zhang et al., 2011). The EV-A species includes 25 serotypes viruses such as the coxsackievirus A group (CVA2–CVA8, CVA10, CVA12, CVA14, and CVA16), human enteroviruses (EV-A71, EV-A76, EV-A89–EV-A91, EV-A114, EV-A119–EV-A121), and non-human enteroviruses (EV-A92, EV-A122–EV-A125) (Simmonds et al., 2020). HFMD is mostly a kind of mild and self-limited disease, but some cases infected by EV-A71 can rapidly develop serious neurological complications and even death (Chen et al., 2019). Since 2007, EV-A71 has caused several large outbreaks of HFMD in the Western Pacific Region, especially in mainland China. Although EV-A71 vaccine has been put on the market since 2016, EV-A71 is still an important pathogen of severe and fatal HFMD (Zhao, 2019). No effective and specific antiviral drugs are available. Severe HFMD caused by EV-A71 remains an important public health problem for children.

As a member of EV-A in the genus Enterovirus of the family Picornaviridae, EV-A71 contains a single-stranded, positive-sense RNA genome of nearly 7.4 kb in length (Lin et al., 2017). The genome of EV-A71 consists of 5′- and 3′-untranslated regions (UTR) and one large open reading frame (ORF). The ORF encodes a large polyprotein divided into 3 regions, P1, P2, and P3. The P1 region encodes 4 structural proteins VP1–VP4, which constitutes the virus capsid, while the P2 and P3 regions encode 7 non-structural proteins 2A–2C and 3A–3D, which participate in polyprotein processing, viral replication, and interactions with the host (Wang and Li, 2019).

Based on phylogenetic analysis of VP1 gene, EV-A71 is currently classified into 8 genotypes A–H (Deshpande et al., 2003; Saxena et al., 2015; Bessaud et al., 2014; Majumdar et al., 2018). Genotype A contained a sole prototype strain BrCr isolated in California, USA, in 1969 (Schmidt et al., 1974) and some strains emerged in China in 2008–2010 (Vakulenko et al., 2019). The emergence of genotype A EV-A71 strains in China was considered as a release of the prototype strain from a laboratory into circulation (Vakulenko et al., 2019). Genotypes B and C are subdivided into subgenotypes B0–B5 and C0–C5, respectively (Noisumdaeng et al., 2018). In the Asia-Pacific region, subgenotypes C1, C2, B3, B4, and C4 prevailed in Australia, subgenotypes B3, B4, B5, C1, and C2 prevailed in Malaysia and Singapore, subgenotypes B2, C1, C2, B4, C4, and B5 prevailed in Japan, and subgenotypes B4, C2, C4, B5, and C5 prevailed in Taiwan (Mizuta et al., 2014). Different from these areas with the alternative circulation of several predominant subgenotypes, C4 was the sole predominant subgenotype in mainland China (Zhang et al., 2013). New genotypes were identified in India (genotypes D and G) (Deshpande et al., 2003; Saxena et al., 2015), sub-Saharan Africa (genotype E) (Bessaud et al., 2014), Madagascar (genotype F) (Bessaud et al., 2014), and Pakistan (genotype H) (Majumdar et al., 2018). EV-A71 shared high genetic diversity in developing countries where the circulating genotypes might be underestimated. Except for C4, genotype A and subgenotypes B5, C2, and C3 had emerged in mainland China (Zhang et al., 2013; Tao et al., 2012; Li et al., 2018; Zhang et al., 2009). This study aimed to know whether there were new genotypes or subgenotypes in mainland China by analyzing the phylogenetics of EV-A71 based on a large dataset of VP1 sequences available from GenBank database. And the genomic characteristics of EV-A71 strains belonging to different subgenotypes were analyzed by recombination analyses.

Materials and methods

Viral sequences

This study involved no human participants or human experimentation. The complete VP1 gene sequences of 3037 Chinese EV-A71 strains covering the time period 1987–2017 were downloaded from the GenBank database. The VP1 sequences of 10 international EV-A71 and CVA16 prototype strains were used as the reference sequences. Three thousand forty-eight VP1 sequences were used in phylogenetic analysis (Table S1). Moreover, the genome sequences of 6 EV-A71 strains whose isolated date was the latest in each subgenotype (C4a1, C4a2, C4b, B5, C2, and C6) were selected as the query sequences for recombination analysis. And another 13 Chinese EV-A71 strains, 10 international EV-A71 strains, and 11 CVA prototype strains (CVA2–8, CVA10, CVA12, CVA14, and CVA16) were used for comparison in recombination analysis. Forty genome sequences were employed in recombination analysis (Table S2).

Phylogenetic analysis

The VP1 gene sequences were converted into a merged file in Clustal format by SeqVerter™. Then, the merged file was converted into MEGA format by MEGA-X. Phylogenetic tree was constructed by MEGA-X using the neighbor-joining method (bootstrap test 1000 replicates) with Kimura 2-parameter model.

Recombination analysis

The merged sequence file was first aligned by MEGA-X with ClustalW method to generate a nucleotide alignment. The similarity plot and bootscan analysis were performed by SimPlot3.5.1. The Kimura 2-parameter distance model and the neighbor-joining tree model were selected in all bootscanning with a window size of 200 nucleotides moving along the alignment in 20-bp increments.

Results

Phylogenetic analysis of the entire VP1 gene of EV-A71 in mainland China

The phylogenetic tree showed that except for 3011 strains clustered into subgenotype C4, 17 strains clustered into genotype A, 4 strains clustered into subgenotype B5, 2 strains clustered into subgenotype C2, 2 orphan strains respectively clustered into subgenotype C0 and C3, and one EV-A71 strain DL71 (GenBank accession number KF982854) clustered into a new branch in genotype C (Fig. 1). The strain DL71 was isolated from Dalian City in Liaoning Province in 2012. It did not fall into any of the six known subgenotypes C0–C5 of EV-A71 but formed a separate genetic branch containing no other strains. The phylogenetic trees based on VP4 and 2A sequences also indicated that DL71 belonged to a previously unknown subgenotype in genotype C (Fig. S1). The new branch was defined as subgenotype C6. The phylogenetic analyses indicated that 7 genotypes or subgenotypes including A, B5, C0, C2, C3, C4, and C6 of EV-A71 circulated in mainland China.

Fig. 1
figure 1

Phylogenetic analysis of VP1 gene sequences of EV-A71 with reference strains. The evolutionary branches were colored differently by the viral isolated time. Clade C4a2 of subgenotype C4 was compressed to save space. The scale bar indicated the number of nucleotide substitutions per site

Full size image

Time and space distribution of subgenotypes of EV-A71 in mainland China

To describe the prevalence of EV-A71 subgenotypes in mainland China, the information including place and time of isolation and subgenotypes of 3037 EV-A71 strains was recorded (Fig. 2). In 1987–1997, EV-A71 exhibited an extremely low prevalence in mainland China. Since 1998, C4 EV-A71 strains emerged, which became the predominant subgenotype. Three shifts in the predominant subgenotype between clades C4b and C4a2 were observed. C4b was the predominant subgenotype in 1998–2002 and 2004. The first shift from C4b to C4a2 was observed in 2003. The second shift from C4a2 to C4b was observed in 2004, followed by the third shift from C4b to C4a2 in 2005. Since then, C4a2 has been the predominant subgenotype. Moreover, different from the sporadic emerging subgenotypes C0, C2, C3, and B5, the genotype A EV-A71 strains co-circulated with both C4b and C4a1 in Central China in 2009, and the subgenotype C6 EV-A71 strain co-circulated with the re-emerged C4a1 in Northeast China in 2012.

Fig. 2
figure 2

Temporal and geographical distribution of EV-A71 subgenotypes in mainland China during 1987–2017. a Temporal distribution of EV-A71 subgenotypes was analyzed by isolated time of EV-A71 strains. b Geographical distribution of EV-A71 subgenotypes was analyzed by the isolated places of EV-A71 strains

Full size image

Recombination analysis of Chinese C4 EV-A71 strains and other EV-A strains

To characterize the genomic features of the predominant subgenotype C4 EV-A71 strains, the complete genome sequences of 3 EV-A71 strains, BZ200805 (HQ694983, C4b, 2008), SHZH03 (AY465356, C4a1, 2003), and 160-50 (MG875331, C4a2, 2016), were selected as the query sequences for similarity and bootscan analysis (Fig. 3). Similarity plot analyses indicated that C4 EV-A71 strains (including C4b, C4a1, and C4a2) exhibited the highest degree of similarity to all C EV-A71 strains in the P1 region; however, in the P2 and P3 regions, they contained an unidentified sequence that was apparently not related to C EV-A71 strains (Fig. 3a, c, and e). Two different similarity patterns for C4 EV-A71 strains were observed in the P2 and P3 regions. One pattern showed that C4 EV-A71 strains resembled CVA4, CVA5, CVA14, and CVA16 in the region from 3700 to 6000 bp, and in the region after 6000 bp, they still resembled CVA4, CVA14, and CVA16 but away from CVA5. The other pattern showed that C4 EV-A71 strains shared high similarity to all B EV-A71 strains in the region from 3700 to 5650 bp, while in the region after 5650 bp, they still shared high similarity to B3 EV-A71 strains but away from B1, B2, B4, and B5 EV-A71 strains.

Fig. 3
figure 3

Recombination analysis of the complete genome of C4 EV-A71 strains in mainland China. The genome sequences of BZ200805, subgenotype C4b (a, b), SHZH03, subgenotype C4a1 (c, d), and 160-50, subgenotype C4a2 (e, f) were used as query sequences in the similarity plot (left panels) and bootscan analysis (right panels). The structural map of EV-A71 genome was displayed upon each panel. The possible crossover breakpoints marked with red vertical line was showed in the analysis

Full size image

Moreover, the bootscan analyses suggested possible recombination events with at least two recombinant breakpoints on C4 EV-A71 corresponding to the 2A/2B junction and 3C/3D junction (Fig. 3b, d, and f). C4 EV-A71 strains contained alternative regions that were closely related to B EV-A71 strains, CVA5, CVA4, CVA14, and CVA16 strains in the region between the 2A/2B and 3C/3D junctions. And in the region after the 3C/3D junction, C4 EV-A71 strains were more closely related to B3 EV-A71, CVA4, CVA14, and CVA16 strains.

Recombination analysis of Chinese B5, C2, C6 EV-A71 strains and other EV-A strains

The genomic characteristics of EV-A71 strains belonging to B5, C2, and C6 were also analyzed. The complete genome sequences of EV-A71 strains CQ2014-86 (KU647000, B5, 2014), 30-2/2015/BJ (MG214681, C2, 2015), and DL71 (KF982854, C6, 2012) were selected as the query sequences for the similarity and bootscan analysis (Fig. 4).

Fig. 4
figure 4

Emergence of EV-A71 subgenotypes B5, C2, and C6 by recombination with other EV-A strains. The whole genome sequences of CQ2014-86, subgenotype B5 (a, b), 30-2/2015/BJ, subgenotype C2 (c, d), and DL71, subgenotype C6 (e, f) were used in the similarity plot (left panels) and bootscan analysis (right panels). The structural map of EV-A71 genome was displayed upon each panel. The possible crossover breakpoints marked with red vertical line was showed in the analysis

Full size image

For B5 strain CQ2014-86 (Fig. 4a and b), the sequence before the 2A/2B junction (nearby 3750 bp) resembled C EV-A71 strains, while the fragment (3750–5662 bp) was more similar to C4 EV-A71, CVA4, CVA5, CVA14, and CVA16 strains than C1-C3 and C5 EV-A71 strains. The fragment (5662–6500 bp) of CQ2014-86 shared higher similarity to CVA3, CVA5, CVA6, CVA10, and CVA12 strains than C EV-A71 strains. The sequence after 6500 bp of CQ2014-86 was more similar to C1–C3 and C5 EV-A71 strains than C4 EV-A71 strains. Two recombination events possibly occurred in B5 EV-A71 strain: one occurred between C4 EV-A71 and CVA3, CVA5, CVA6, CVA10, and CVA12 strains with breakpoint located in the 3C gene, and the other occurred between C1–C3 and C5 EV-A71 and CVA4, CVA5, CVA14, and CVA16 strains with two breakpoints located in the 2A/2B junction and 3C gene.

For C2 strain 30-2/2015/BJ (Fig. 4c and d), the results provided strong evidence that it might be a recombinant virus between C4 EV-A71 and CVA8 strains with two breakpoints located in the 5′UTR (550 bp) and 2A/2B junction (3750 bp). The sequence (550–3750 bp) of 30-2/2015/BJ shared the highest similarity to C4 EV-A71 strains, while the 5′UTR (1–550 bp) and the sequence after 3750 bp of 30-2/2015/BJ shared the highest similarity to CVA8 strain.

For C6 strain DL71 (Fig. 4e and f), the similarity plot analysis exhibited an irregular helix-like recombination on DL71 formed by the winding of the genome sequences of C4 and C2 EV-A71 strains with multiple breakpoints. In the P1 region, DL71 shared higher similarity to C4 and C2 EV-A71 strains (>80%) than B EV-A71 strains (65–85%) and other EV-A strains (<75%). In the P2 and P3 regions, 4 fragments (4022–4390 bp, 4650–5830 bp, 6190–6560 bp, and 6890–7240 bp) of DL71 shared higher similarity to C4 EV-A71 strains (>80%) than C2 EV-A71 strain (<80%), while other 5 fragments (3700–4022-bp, 4390–4650-bp, 5830–6190-bp, 6560–6890-bp, and 7240-bp 3′UTR) of DL71 shared higher similarity to C2 EV-A71 strain (>80%) than C4 EV-A71 strains (<80%).

Discussion

In this study, the genetic diversity of EV-A71 strains in mainland China was well demonstrated by the molecular epidemiology and recombination analyses. Besides the predominant subgenotypes C4, some other sporadic subgenotypes A, B5, C0, C2, and C3 and the new subgenotype C6 emerged in mainland China. The first C2 EV-A71 strain 96200/SD/CHN/96 in mainland China isolated from an acute flaccid paralysis case shared high identity in VP1 with the second C2 EV-A71 strain 30-2/2015/BJ (Tao et al., 2012; Li et al., 2018). The later C2 strain was isolated from a child who had traveled back from Thailand, where C2 was commonly detected (Li et al., 2018). The C2 strains in mainland China might be imported from other countries or areas. In mainland China, the first B5 strain was detected in 2009, and then 3 strains were isolated in 2012, 2014, and 2016. Subgenotype B5 strains were first detected in Taiwan in 1998 (Wu et al., 2013) and then circulated widely in the Asia-Pacific region. In 2009, subgenotype B5 strains simultaneously circulated in Japan, Malaysia, and Taiwan. Frequent international travel might promote the import of sporadic viruses into mainland China. Moreover, a new finding in this study was the emergence of C6 EV-A71. One orphan strain DL71 was clustered into the new subgenotype C6 different from the previously known subgenotypes C0–C5 in genotype C. The results suggested that several subgenotypes of EV-A71 had existed in China. Thus, the emergence and importation of new genotypes or subgenotypes of EV-A71 should be continuously monitored (Li et al., 2018; Noisumdaeng et al., 2018).

Genetic recombination plays important roles in the genetic variation and evolution of viruses. Several studies have provided evidence of recombinant events in EV-A71, including intratypic recombination between different subgenotypes of EV-A71 and intertypic recombination between EV-A71 and other EV-A strains (Zhang et al., 2011; Zhang et al., 2013; Yoke-Fun and AbuBakar, 2006; Wang et al., 2012). However, the query sequences of EV-A71 selected in recombination analysis in previous studies were mostly the consensus sequences or the oldest strains within different subgenotypes (Zhang et al., 2013), and few sequences of newly isolated EV-A71 strains were used as the query sequences. In the present study, the complete genome sequences of the latest EV-A71 strains within each subgenotype were used as the query sequences. Our comprehensive recombination analysis exhibited two patterns of recombination in C4 EV-A71 strains (Fig. 3). We found that intertypic recombination of C4 sequences occurred between the structural genes of C EV-A71 and the non-structural genes of CVA4, CVA5, CVA14, and CVA16, which was consistent with the findings in previous reports (Zhang et al., 2011; Zhang et al., 2013). However, the evidence of intratypic recombination between EV-A71 genotypes C and B was not sufficient (Zhang et al., 2013). In addition to intertypic recombination, our study also showed evidence of intratypic recombination between EV-A71 genotypes C and B. The fragment (3700–5650 bp) of C4 EV-A71 resembled B EV-A71, and this recombination pattern was also found in other C4 EV-A71 strains isolated in Thailand (Noisumdaeng et al., 2018; Mauleekoonphairoj et al., 2015), Taiwan (Huang et al., 2008), and Hong Kong (Yip et al., 2013). Moreover, we found interestingly that except for the fragment (3700–5650 bp) of C4 being closer to B3, the sequence after 5650 bp of C4 EV-A71 still shared high similarity to B3 but away from B1, B2, B4, and B5. A phenomenon in a previous study that the P2 and P3 regions of B3 EV-A71 had high similarity to CVA4, CVA14, and CVA16 (Yoke-Fun and AbuBakar, 2006) also supported our results. Thus, there was strong evidence to support the conclusion of the intertypic recombination between genotype C EV-A71 strains and other EV-A strains and the intratypic recombination between EV-A71 genotypes C and B. Several studies renamed this “double-recombinant” C4 EV-A71 strain genotype D (Yip et al., 2010; Yip et al., 2013).

The genomic characteristics of EV-A71 strains belonging to B5, C2, and C6 in mainland China were also analyzed (Fig. 4). The recombination manner of C2 EV-A71 strain 30-2/2015/BJ in mainland China was similar to other C2 EV-A71 strains isolated in Malaysia and Taiwan (Yoke-Fun and AbuBakar, 2006), which resembled C4 EV-A71 strains in the P1 region and resembled the prototype strain of CVA8 in the 5′UTR, P2, and P3 regions. The B5 EV-A71 strain CQ2014-86 was closer to C EV-A71 strains in the P1 region but contained an unidentified sequence in the P2 and P3 regions. Given the global distribution of EV-A71 subgenotypes C2 and B5 (Mizuta et al., 2014; Wu et al., 2013), the sporadic C2 and B5 strains that emerged in mainland China were most likely to have been imported from other countries and areas. An interesting finding was that DL71 (Dalian City, 2012) was a recombinant between the dominant subgenotype C4 and the sporadic subgenotype C2 EV-A71 strains with a helix-like recombination manner. The details of recombination events in DL71 needed further studies.

The HFMD outbreaks and EV-A71 circulation in mainland China provided conditions for C6 EV-A71 emergence. Firstly, frequent outbreaks of HFMD occurred in Dalian City with the highest number of cases in 2012 during the period of 2010–2013 (Gao et al., 2016; Wang et al., 2019). Secondly, as one parent strain of C6, C4 was the single predominant subgenotype in mainland China with much genetic diversity (including C4a1, C4a2, and C4b). Thirdly, as another parent strain of C6, C2 has been widely detected in many countries for several years (World Health Organization [WHO], 2011). The importation of C2 from neighboring countries led to its sporadic circulation in mainland China. Thus, the prevalence of HFMD, the genetic diversity of the predominant subgenotype C4, and the sporadic circulation of C2 might together promote the recombination and emergence of C6 EV-A71. Moreover, no evidence was reported about the clinical symptoms and disease severity of DL71. Based on reviewing the potential virulence determinants of EV-A71 in previous reports (Huang et al., 2015; Kung et al., 2010; Li et al., 2016; Li et al., 2011; Liu et al., 2014b; Ma et al., 2018; Meng and Kwang, 2014; Wen et al., 2013; Yuan et al., 2015), the nucleotides and amino acids located at these virulence sites in DL71 were analyzed (Table S3). In DL71 genome, two nucleotides and four amino acids might be with high virulence at six sites, while six nucleotides and twelve amino acids might be with low virulence at eighteen sites. In fact, the pathogenicity of EV-A71 was very complex and influenced by many factors including virus virulence and host immunity. It was difficult to judge the pathogenicity of C6 EV-A71. This was one limitation in the study.

Our present findings that recombination breakpoints in the EV-A71 genome were frequently located in the non-structural regions 5′UTR, P2, and P3 but rarely in the structural region P1 were consistent with previous reports (Zhang et al., 2011; Yip et al., 2010; Noisumdaeng et al., 2018; Zhang et al., 2013; Yoke-Fun and AbuBakar, 2006; Huang et al., 2008; Zhang et al., 2010; Wang et al., 2012; Li et al., 2012; Yip et al., 2013; Liu et al., 2014a; Mauleekoonphairoj et al., 2015). On one hand, the P1 region of EV-A71 maintained a high degree of similarity within EV-A71 strains, but away from other EV-A strains, possibly because of the need for structural constraint to ensure the proper assembly of viral capsid (Wang et al., 2012). While recombination occurred in the structural genes, it might produce nonviable or unstable progeny viruses that could not survive (Stanway et al., 1986). On the other hand, the reason for the non-structural regions P2 and P3 as the hot spots for recombination in EV-A71 was that the P2 and P3 regions of EV-A shared more homology than the structural region (Wang et al., 2012). The acquisition of non-structural genes encoding replicative components such as protease, helicase, and polymerase from B EV-A71 strains and other EV-A strains (CVA4, CVA5, CVA14, and CVA16) in C4 EV-A71 strains might provide evolutionary advantages for the “double-recombinant” strains (Yip et al., 2013). In addition, the inconsistent topology of phylogenetic trees suggested that it was not sufficient to genotype the recombinant viruses by VP1 alone. Given that these recombinant viruses might be a threat to public health in mainland China, it was still worthwhile to closely monitor the prevalence of EV-A71 and the emergence of new subgenotypes by performing extensive studies on the viral complete genome.