Molecular epidemiology of coxsackievirus A16 circulating in children in Beijing, China from 2010 to 2019

Background Coxsackievirus A16 (CVA16) is one of the major etiological agents of hand, foot and mouth disease (HFMD). This study aimed to investigate the molecular epidemiology and evolutionary characteristics of CVA16. Methods Throat swabs were collected from children with HFMD and suspected HFMD during 2010–2019. Enteroviruses (EVs) were detected and typed by real-time reverse transcription-polymerase chain reaction (RT-PCR) and RT-PCR. The genotype, evolutionary rate, the most recent common ancestor, population dynamics and selection pressure of CVA16 were analyzed based on viral protein gene (VP1) by bioinformatics software. Results A total of 4709 throat swabs were screened. EVs were detected in 3180 samples and 814 were CVA16 positive. More than 81% of CVA16-positive children were under 5 years old. The prevalence of CVA16 showed obvious periodic fluctuations with a high level during 2010–2012 followed by an apparent decline during 2013–2017. However, the activities of CVA16 increased gradually during 2018–2019. All the Beijing CVA16 strains belonged to sub-genotype B1, and B1b was the dominant strain. One B1c strain was detected in Beijing for the first time in 2016. The estimated mean evolutionary rate of VP1 gene was 4.49 × 10–3 substitution/site/year. Methionine gradually fixed at site-23 of VP1 since 2012. Two sites were detected under episodic positive selection, one of which (site-223) located in neutralizing linear epitope PEP71. Conclusions The dominant strains of CVA16 belonged to clade B1b and evolved in a fast evolutionary rate during 2010–2019 in Beijing. To provide more favorable data for HFMD prevention and control, it is necessary to keep attention on molecular epidemiological and evolutionary characteristics of CVA16. Supplementary Information The online version contains supplementary material available at 10.1007/s12519-021-00451-y.


Introduction
Enterovirus (EV) is a common pathogen, and is responsible for many infectious diseases in human, including hand, foot and mouth disease (HFMD), herpangina (HA), acute hemorrhagic conjunctivitis, respiratory infections, acute myocarditis, meningitis, encephalitis and acute flaccid paralysis [1]. HFMD was classified as a statutorily notifiable infectious disease in China in 2008 [2]. The clinical manifestations of HFMD are highly complex and heterogeneous, which makes it difficult for doctors to give an exact clinical diagnosis, especially at the early stage of disease. The main manifestations of HFMD are fever and rash on hands, feet, mouth, and buttocks, however, central nervous system complications and cardiopulmonary failure may occur in severe cases [3]. Therefore, early recognition of cases and identification of the pathogens are the key to prevention and control for HFMD.
Coxsackievirus A16 (CVA16) is a member of species Enterovirus A in genus Enterovirus, family Picornaviridae. As one of the main pathogens of HFMD, CVA16 was responsible for several HFMD outbreaks in the world, especially in Asia-Pacific region [4][5][6]. HFMD caused by CVA16 infection is generally mild and self-limiting. However, CVA16 can occasionally cause severe and fatal cases [7]. CVA16 contains a positive sense, single-stranded RNA genome [8]. The nucleotide sequence of viral protein gene (VP1), encoding the most important structural protein, is well correlated with EV serotype and genotype [9,10]. Up to now, global CVA16 strains can be divided into three genotypes, A, B and D based on the phylogenetic tree and genetic diversity of VP1 gene. Genotype B can be further divided into three sub-genotypes, B1, B2 and B3. Sub-genotype B2 was the predominant type before 2000, then replaced by B1 strains [11]. Sub-genotype B1 contains clade B1a, B1b and B1c. Two novel CVA16 strains isolated by Chen et al. in 2017 in Shenzhen, China, were designated as sub-genotype B3 by the phylogenetic reconstruction of VP1 [12]. Genotype D was first detected in Peru in 2009, then circulated in France from 2011 to 2014 [13]. In 2016, the first outbreak and spread of genotype D in China was reported in Shanghai [14].
Our laboratory has been conducting etiological surveillance for HFMD and EV-associated infectious diseases since 2007 [15,16]. During the monitoring period, more types of EVs apart from EV-A71 and CVA16 were further identified along with the continuous modification of primers and polymerase chain reaction (PCR) amplification conditions [17]. This study focused on the molecular epidemiology and evolutionary features of CVA16 circulating in children with HFMD and suspected HFMD in Beijing from 2010 to 2019.

Patients selection and samples collection
Patients involved in this study visited the Department of Infectious Diseases, Children's Hospital of Capital Institute of Pediatrics during the period from March 2010 to October 2019. Throat swabs were collected from patients who were under 18 years old with clinical diagnoses of HFMD, HA and rash and fever illness. Among 4709 patients, the ratio of male and female was 1.39:1. The mean age was 3.35 years (range from 9 days to 17 years 6 months).

Enterovirus detection and typing
The processes for EV detection and typing have been described previously [17]. Briefly, viral RNA was extracted from clinical samples using QIAamp Viral RNA Mini Kit (QIAGEN, Germany), and real-time reverse transcription (RT)-PCRs for pan-EV, EV-A71, CVA16, CVA6 and CVA10 were performed. Samples, which were only positive for pan-EV, were further typed using RT-PCR and sequencing.

VP1 gene amplification and sequencing
CVA16-positive samples were selected randomly and proportionally each year to amplify the complete VP1 gene using primers CVA16-VP1-F and CVA16-VP1-R (or CVA16-VP1-F1 and CVA16-VP1-R1) ( Table 1). The positive PCR products were sequenced by Sino Geno Max Co. Ltd. (Beijing, China). The VP1 sequences (891 bp) were edited using DNAStar v. 5.01 software and compared with publicly available sequences in GenBank using BLAST (http:// www. ncbi. nlm. nih. gov/ BLAST/). The similarity and divergence of nucleotide and deduced amino acid sequences were estimated by MEGA v. 6.05 software (P-distance model) [18].

Phylogenetic analysis
Phylogenetic tree was constructed by neighbor-joining (NJ) and maximum likelihood (ML) methods based on Kimura 2-parameter model using MEGA v. 6.05 software. The bootstrap analyses with 1000 repetitions were performed to estimate the reliability of the phylogenetic inference at each branch node.

Molecular evolution and population dynamics analysis
To assess whether there was sufficient temporal signal in the VP1 sequence dataset to proceed with molecular clock analysis, a regression analysis of root-to-tip genetic distance against sampling date based on the ML tree of Beijing strains was performed using TempEst v. 1.5.3 software [19]. The Bayesian molecular clock phylogeny and evolutionary rate of CVA16 VP1 genes were analyzed using Markov chain Monte Carlo method in BEAST v. 1.10.4 software [20]. The chain length was 100 million steps with sampling every 10,000 steps. The substitution model was selected

Selection pressure analysis
Natural selected sites were inferred by mixed-effects model for episodic diversifying selection (MEME), single likelihood ancestor counting (SLAC) and fast unconstrained Bayesian approximation (FUBAR) methods implemented in Datamonkey online server (http:// datam onkey. org/) [22][23][24]. MEME employs a mixed-effects maximum likelihood approach which is capable of identifying both episodic and pervasive selection at the level of an individual site. SLAC uses a combination of maximum likelihood and counting approaches with the most rigorous test result. FUBAR uses a Bayesian approach which can avoid misleading inference due to model misspecification. Both SLAC and FUBAR assume that the selection pressure for each site is constant along the entire phylogeny, estimating pervasive selection at the level of an individual site. Positively selected sites were determined by a P value of < 0.05 (MEME and SLAC methods) or a posterior probability of > 0.9 (FUBAR method).

Epidemiology
A total of 4709 throat swabs were collected from children with clinical diagnoses of HFMD, HA and rash and fever illness from March 2010 to October 2019. EVs were detected in 67.5% (n = 3180) of samples, and 17.3% (n = 814) were positive for CVA16. Figure 1a shows the annual number of collected cases and EV-positive cases. The number of screening and laboratory confirmed EV-positive cases after 2013 significantly increased because of the HFMD outbreak of CVA6 in China and the atypical clinical manifestation caused by CVA6 [25,26]. During 2010-2012, the detection rates of CVA16 were at high level (42.4%, 39.0%, 47.9%, respectively), while the detection rate of CVA6 was extremely low. In 2013, the detection rate of CVA16 dropped rapidly (18.4%), along with a sharp rise of CVA6 (48.2%) (Fig. 1b). In 2014, the detection rate of CVA16 increased again, and that of CVA6 dropped rapidly. With a fluctuating decline, the detection rate of CVA16 decreased to the lowest level in 2017 (4.6%), followed by an increase to 12.6% Among the CVA16-positive children, 789 cases with complete demographic information were selected to analyze age and sex distributions. Male to female ratio was 1.28:1. The median age was 3.0 years old (range from 1 month to 14 years 2 days). 81.6% (644/789) were younger than 5 years old. Most of them were in the age of 1-3 years old. The proportion of CVA16 infection declined with age among children older than 4 years old (Fig. 1d). Phylogenetic tree was constructed using NJ method based on the alignment of 164 unduplicated VP1 sequences generated in this study. A total of 95 reference sequences were obtained from GenBank database including the VP1 sequences of prototype G-10 strain, global CVA16 strains and the prototype strain of EV-A71 (BrCr). Five additional sequences collected in Beijing between 2010 and 2014 (obtained from GenBank database) were also included. The NJ tree in Supplementary Fig. 1 shows that all the 169 Beijing CVA16 strains belonged to B1 sub-genotype, and attached to 3 clades (B1a, B1b and B1c). Overall, 16 (9.0%) out of 177 Beijing strains (including the 8 duplicate sequences  (Table 2).

Phylogenetic analysis of VP1 sequence
Nucleotide and deduced amino acid divergences of VP1 were described in Table 3. All of the Beijing and reference strains including the prototype EV-A71 strain (BrCr) were analyzed. As a result, the mean nucleotide divergences between genotypes A, B and D were higher than 14.0%, the mean nucleotide divergences between sub-genotypes B1, B2 and B3 were greater than or equal to 10.0%, and less than 10.0% of the mean nucleotide divergence was existed between clades B1a, B1b and B1c. The divergences of the deduced VP1 amino acid sequences between genotypes A and B, A and D were obvious; however, the divergence between genotypes B and D was small, similar to that of intra-genotype B, indicating the closer evolutionary relationship between genotypes B and D. The consensus amino acid variations of VP1 in different branches were shown on the right in Supplementary Fig. 1. The amino acid at site-13 was different in the strains among different genotypes A, B and D. At site-14, a separate cluster of B1b strains had a consistent N14S mutation, becoming the same amino acid as genotype A strains. A wider cluster of B1b strains contained a consistent mutation (L23M) at site-23. The amino acid at site-23 of the B1b strains collected in Beijing from 2010 to 2011 was mainly leucine (L), whereas most of the strains emerged L23M mutation in 2012 (Table 4). Methionine (M) at site-23 was fixed in the majority of the B1b strains collected from 2013 to 2019. All strains with the amino acid mutation mentioned above belonged to cluster 3 in Supplementary Fig. 2a. Compared with other strains, the B1c strains had synchronous mutations of I235V and T240A.

Molecular evolution and population dynamics
For the 169 Beijing CVA16 strains collected during 2010-2019, a Bayesian MCC tree based on VP1 gene was generated to confirm the evolutionary relationship and to explore the timescale of CVA16. The correlation coefficient calculated by TempEst was 0.899, and the root-to-tip plot showed a positive correlation between genetic divergence and sampling time ( Supplementary Fig. 3), suggesting that the dataset was suitable for phylogenetic molecular clock analysis. As shown in Supplementary Fig. 2a, the time-scaled MCC tree contains three clades corresponding to clades B1a, B1b and B1c in Supplementary Fig. 1. Clade B1b can be further divided into 3 clusters (cluster 1, cluster 2 and cluster 3). The strains that emerged earlier (2010-2013) are mostly located at the root of the MCC tree. Most strains (cluster 3) that emerged later evolved from early strains, and further formed multiple new branches. Cluster 3 of clade B1b became the predominant strains since 2014, and formed a sustained epidemic through multiple transmission chains in the local area. However, there was no obvious ladder-like topological structure for clade B1a on MCC tree, indicating that B1a strains had a limited prevalence in Beijing, most of which may be imported from the other regions. Compared the epidemiological data ( Fig. 1b) with BSP, the

Selection pressure analysis of CVA16 strains in Beijing
Over the past decade, CVA16 continuously circulated in Beijing. To understand the natural selection pressure on VP1 gene, 169 unduplicated VP1 nucleotide sequences of Beijing CVA16 strains described above were analyzed by Datamonkey. The results of SLAC and FUBAR analyses showed no positive selection site on VP1 and most sites were under negative selection (204 and 275 sites, respectively). The result of SLAC analysis showed that the overall dN/dS ratio of VP1 gene is very low (0.02). These results above suggest that VP1 gene of Beijing CVA16 strains was subjected to very strong purifying selection. However, two sites under episodic positive selection were detected by MEME. The positively selected amino acid substitutions were D2V and G223N distributed in two strains (S854 sampled in 2010 and S5822 sampled in 2018) (marked in Supplementary Fig. 1). The mutation at site-223 was located on the linear neutralizing epitope on the surface of the virus capsid (PEP71, aa. 211-225 [27]), resulting in the change of non-polar hydrophobic amino acid (G) to polar neutral amino acid (N).

Discussion
This retrospective study investigated the molecular epidemiology and genetic evolutionary characteristics of CVA16 circulating in children diagnosed with HFMD and suspected cases in Beijing during 2010-2019. Over the past decade, the prevalence of CVA16 has changed significantly. CVA16 and EV-A71 were replaced by CVA6 in 2013 and during 2015-2019. On the whole, CVA16 had a fluctuating prevalence pattern in Beijing from 2010 to 2018, with the detection rate increasing every other year. The above prevalence pattern also existed in the other regions of China [28,29]. CVA6 and CVA10 also had the same trend, but the detection rate of CVA10 was always at a low level. However, different countries may have different epidemic patterns due to the differences in climate and medical-health conditions [30]. The detection rate of EV-A71 dropped obviously during 2017-2019, in line with another research about the application of EV-A71 vaccine in China [31]. However, the prevalence of CVA16 showed an upward trend from 2018 to 2019. It suggests that EV-A71 vaccine has limited cross-protection effect on CVA16. Therefore, it is necessary to use a multivalent vaccine that can prevent common types of EVs. Most children with CVA16 infection aged younger than 5 years, and younger children (1-3 years old) accounted for a larger proportion, consistent with HFMD-related studies in China and other Asian countries [32,33]. The global CVA16 strains can be divided into 3 genotypes, and genotype B can be divided into 3 sub-genotypes. In recent years, sub-genotype B1 was the prevalent strain worldwide [4], also contributed to all the CVA16 infections in this study. Clade B1a was first reported in Japan in 1995, then became the dominant strain in mainland China, Malaysia and Thailand [34][35][36]. B1b, the complex recombinant of CVA16 with CVA4 and EV-A71, co-circulated with B1a in China during 1999-2008 [37,38]. After 2010, B1b became the predominant strain in China [12,39,40]. Clade B1c strain was first emerged in Malaysia (2005Malaysia ( -2007, then detected in several countries, such as Russia, France, Indian and Japan, and induced an outbreak of HFMD in Indian in 2013 [34,41,42]. Recently, a B1c strain sampled in 2017 was reported in the west area of China [29]. In this study, clade B1b was the dominant strain during 2010-2019, circulated in Beijing and evolved continuously, while the prevalence of clade B1a and B1c strains were limited. Clade B1c and genotype D of CVA16 were newly emerged types in recent years, and transmitted in several countries [13,14,42]. Genotype D strains experienced multiple recombination during transmission, which might change the biological characteristics of the virus, then change the spreading ability and pathogenicity of the virus [43]. Therefore, it should be paid close attention to the molecular epidemiology of CVA16 to timely detect and control the potential HFMD outbreak caused by the new genotypes of CVA16. The divergence of VP1 nucleotide sequences greater than 15% is one of the reference basis for CVA16 genotyping [10,38]. However, in this research, the divergence between strains of genotype B and D is slightly lower than 15%, and the divergence of amino acid sequences was relatively low. It might be due to the recombinant origin of genotype D [13]. The considerable overlap between the nucleotide sequence divergence of sub-genotypes (B1-B3) and clades (B1a-B1c) suggested that it was difficult to distinguish sub-genotype and clade only by the divergence of VP1 nucleotide sequences.
The BSP based on VP1 gene of CVA16 can reflect the pathogen population dynamics. The BSP reconstructed using CVA16 VP1 gene worldwide sampled from 1980 to 2013 showed the relationship between the viral genetic diversity and the outbreak of HFMD caused by CVA16. Every sharp increase of CVA16 genetic diversity reflected the emergence of a new genotype, which resulted in a large-scale HFMD outbreak [41]. In mainland China, CVA16 genetic diversity increased from 2000 to 2009, then decreased gradually from 2010 to 2013. From 2014 to 2016, it continued to increase again, then decreased slightly from 2017 to 2018 [29]. The genetic diversity correlated with EPS. There was a similar result in this study. The EPS was at a relatively high level from 2014 to 2016, but the detection rate of CVA16 did not reach a correspondingly high level. It may be due to the increased focus on atypical HFMD cases caused by CVA6 and increased efforts to screening for EVs after 2013. It should be noted that the EPS and prevalence trend of CVA16 continued to increase from 2018 to 2019 in Beijing, therefore, implying the enhancement of CVA16 activity and the raising alertness about the potential outbreak.
Site-23 of VP1 is located at the N-terminus of the protein, which is related to the stabilization of the capsid structure and the process of RNA release [45,46]. The prevalence of CVA16 in Beijing was related to the spread of cluster 3 of clade B1b from 2014 to 2019, all of which had L23M mutation. This mutation was also found in an outbreak of CVA16 in Taizhou from 2011 to 2013 [47]. Whether it is beneficial for virus transmission remains to be analyzed in depth.
VP1 protein of EV contains major epitopes and is mainly affected by the selective pressure of the host's immune system. Most mutations were not beneficial for immune escape of the virus. One site under episodic positive selection pressure (site-223) was located in the linear neutralizing epitope PEP71 [27]. Whether this mutation (G223N) emerged in 2018 could increase the chance of immune escape of the virus and continue to accumulate in the virus population requires close monitoring. The site subject to positive selection may affect the protective effect of the vaccine [48]. Therefore, it is necessary to find suitable neutralizing epitopes to develop effective subunit vaccine, nucleic acid vaccine, etc. Other neutralizing epitopes (PEP32, PEP37, PEP55, PEP63 and PEP91) of VP1 protein, under purifying selection pressure, were more suitable as candidate epitopes for vaccine design than PEP71.
This study had two limitations. First, the detailed clinical information of CVA16-positive patients was incomplete, which limited the further analysis of the clinical symptoms. Second, only local CVA16 strains were involved in this study to analyze the evolutionary characteristics, which limited to understand the transmission and evolution worldwide.
In conclusion, this study summarized the molecular epidemiological and genetic evolutionary characteristics of Beijing CVA16 strains circulating between March 2010 and October 2019. Although CVA16 is not the predominant pathogen of HFMD in recent years, there could be an outbreak of CVA16 in the future due to the lack of available vaccine, accumulation of the susceptible host population, continuous evolution of the virus and import of new genotypes. Therefore, it is necessary to keep attention on molecular epidemiological and evolutionary characteristics of CVA16.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.