The emerging sub-genotype C2 of CoxsackievirusA10 Associated with Hand, Foot and Mouth Disease extensively circulating in mainland of China

Coxsackievirus A10 (CV-A10) associated with Hand, foot, and mouth disease (HFMD) cases emerged increasingly in recent years. In this study, the samples from nation-wide HFMD surveillance, including 27 out of 31 provinces in China were investigated, and the continuous and extensive virological surveillance, covered 13 years, were conducted to provide a comprehensive molecular characterization analysis of CV-A10. 855 CV-A10 viruses (33 severe cases included), were isolated from HFMD children patients during 2009 to 2016 in China. 164 representative sequences from these viruses, together with 117 CV-A10 sequences downloaded from GenBank based on entire VP1 were recruited in this study. Two new genotypes (F and G) and two sub-genotypes (C1 and C2) were identified. Among 264 Chinese sequences, 9 of them were genotype B, 8 of them were C1, and the other (247) were C2, the predominant sub-genotype in China since 2012. Chinese C2 viruses showed obvious temporal characteristics and can be divided into 3 clusters (cluster 1~3). Cluster 3 viruses was circulating extensively during 2014 and 2016 with more severe cases. It is very necessary and important to continuously conduct the extensive virological surveillance for CV-A10, and further evolutionary studies will provide more evidence on its evolution and virulence.

occurred in 2008 4,5 . In the same year, CV-A10 was reported as the most prevalent virus causing the outbreak of HFMD in the Singapore 6 . In mainland of China, following EV-71 and CV-A16, CV-A10 was the third most common virus detected in HFMD during 2009-2011 7 . Most diseases associated with CV-A10 were wild and self-limiting, but there also had severe 8 and death cases 9 been reported in domestic and overseas.
Several studies on CV-A10 genotyping were based on 5′-UTR 10 , VP4 11 , or partial VP1 3,4,12,13 . However, VP1 plays a critical role in mediating binding receptor and the complete VP1 has been used widely to identify EV serotypes [14][15][16][17][18] as its coding region contains many important neutralizing antigenic sites. According to another recent study, the cell surface molecule KREMEN1 was verified as an entry receptor for CV-A10 19 , KREMEN1 overexpression enhances CV-A10 binding to the cell surface and increases susceptibility to infection, indicating that KREMEN1 is a rate-limiting factor for CV-A10 infection.
Although CV-A10 occupied certain proportion in HFMD pathogenic spectrum, its prevalence, genetic character and the criterion for genotyping are still not clear. In this study, the samples from nation-wide HFMD surveillance, including 27 out of 31 provinces in China were investigated, and continuous and extensive virological surveillance were conducted to cover more than 13 years to provide a comprehensive epidemiological and molecular characterization analysis of CV-A10, and to establish the standardization for genotyping.  Table 1). The number of CV-A10 cases were increasing, and the proportion of CV-A10 associated with severe cases were on the rise as well during 2012 and 2016, as shown in Table 1. 855 CV-A10 viruses from 26 out of 31 provinces obtained in this study, together with 189 Chinese CV-A10 sequences downloaded from GenBank during 2004-2016 were included in this study for further analysis.  Table S1). The data showed that CV-A10 has been extensively circulating around the mainland of China since 2012, covered 7 geographic regions and 27 out of 31 provinces of mainland of China (one additional province involved in GenBank) during 2004-2016. While before 2012, CV-A10 was distributing in limited regions and provinces with much less number of cases.

Temporal and geographic distribution of CV-A10 isolates circulated in
Thirty-three out of 855 cases were diagnosed as severe cases: patients with neurological complications or cardiopulmonary complications according to the National Guideline for Diagnosis and Management of HFMD cases issued by National Health and Family Planning Commission of China in 2010 (http://www.nhfpc.gov.cn/zwgkzt/ wsbysj/201004/46884.shtml). Most of patients (716/855, 83.84%) were 1~4 years old, the median age of these patients was 2.00 years (range 0.1-24 years), while 95.12% of severe cases were less than 3 years old. Nearly 3/4 (631/855, 73.81%) of cases occurred during April and July (Fig. 1B), when was spring and summer time in China.
Seven genotypes of CV-A10 were assigned based on entire VP1. There were 284 VP1 sequences available from GenBank before November 15 th , 2017. All these VP1 sequences were from 9 countries during 1950 and 2016, including the prototype virus Kowalik strain. The representative sequences were selected to include the sequences covering all the location/country/provinces and time range, meanwhile the sequences with high homologies or with significant errors in sequences, were not included in this study. In all, total 281sequences were included to be performed the phylogenetic tree, including 164 out of 855 VP1 sequences in this study and 117 out of 284 sequences from GenBank.
A phylogenetic tree was constructed based on these representative viruses of CV-A10, and the sequences were assigned to seven genotypes (A, B, C, D, E, F and G) with at least 14.97% nucleotide divergence between different genotypes (Fig. 2). And two new genotypes (F and G) were assigned in this study, together with previously assigned genotypes A-E 18,20 . Genotype F was consisted of two Indian sequences in 2005, although, these two sequences were not full length of VP1, with 7 bases or 22 bases deletion in 3′ terminal of VP1 region in individual Indian sequences, respectively. One Taiwan China sequence isolated in 2008 was grouped separately, assigned genotype G, and it seems that it circulated only in Taiwan China (Fig. 3). Pairwise nucleotide and amino acid sequences identities in VP1 region between different genotype of viruses in this study were showed in following table (Table 2).
Year All

CV-A10
CV-A10/ Other EV Representative viruses of CV-A10 on partial VP1 in the world. As complete VP1 sequences of CV-A10 available were limited, in order to investigate CV-A10 virus's origin and spread in the world, and to verify the genotypes on entire VP1, the representative viruses of CV-A10 on partial VP1 were selected and performed the phylogenetic analysis (Fig. 3). Based on 239 bps of VP1, we found that CV-A10 can be divided into 9 genotypes, with the mean group distance 14.5-24.6%. Compared with the entire VP1 region tree, those selective partial VP1 sequences were still in the same cluster, except for two potential genotypes (H and I) isolated in earlier year but been ignored. Genotype H was composed of sequences before 2001 that emerged in Japan, Finland and Germany, and genotype I was only composed with one Egypt strain of 2007, which was 18.70% nucleotide difference with the prototype strain Kowalik. has been evolved as the absolutely predominant genotype since 2009, and can be further divided into C1 and C2 sub-genotypes. Sub-genotype C1 was circulating during 2008-2012, while C2 was the predominant sub-genotype after 2010. The nucleotide divergence between sub-genotypes C1 and C2 was 5.8%. All the other Chinese viruses grouped in C2 and can be categorized into 3 clusters: C2-cluster1, C2-cluster2, and C2-cluster3.
The group mean distance among these 3clusters were range from 3.7% to 6.2%. Chinese viruses showed obvious temporal characteristics ( We have compared the nucleotide and amino acid mutations among the sub-genotype C1, C2, and C2 cluster 1~clus-ter3, described in the Table 3 (Table 3). In addition, one synonymous mutation on nucleotide site 813 was also observed: most of sequences from C2 cluster 2 and cluster 3 had the mutation of cysteine (C) to glycine (G), while the C1 and cluster 1 were barely found. Homologous  Supplementary Table S1.
comparison analysis on nucleotide sequences showed that there was no difference between mild cases and severe cases in amino acid coding region. However, we found all the severe cases located in C2 except for one Ningxia strain in 2010, and nearly half of severe viruses (17/33) located in cluster3, which was continuous circulating during 2014 and 2016.  The average evolutionary divergence of genotype C sequences between different time groups was calculated using MEGA 5.03 with p-distance model ( Table 4). The results revealed that the nucleotide divergences within individual year from 2008-2016 were increasing gradually, except for 2009, as limited viruses available in 2009. Analysis of all these viruses did not show obvious regional distribution characteristics (Fig. 3B).

Discussion
HFMD is a common disease among children, especially for children below five years. Although EV-71 and CV-A16 were the predominant pathogens caused HFMD in recent years, the proportion of other EV has been increasing 7,20-22 , among which, CV-A10 was reported to be associated with sporadic HFMD cases and outbreak events in several countries 3,5,6,12 and topped the list of other EV in many provinces of China 23,24 . In this study, we found an apparently increasing proportion of CV-A10 among HFMD cases in mainland of China since during 2012 and 2016, except for 2013 and 2015, when CV-A6 was responsible for most of HFMD outbreaks in China [25][26][27] . This reveals that CV-A10 is gradually becoming one of the most frequent EV serotypes during the epidemic interval of EV-71 and CV-A16. As most of the studies were performed based on province level with a small sample size. This study, we collected the nationwide representative samples and covered over a 13 years' time span in order to provide a comprehensive epidemiology and molecular characterization analysis of CV-A10.
All age groups are susceptible to be infected with CV-A10, but most cases have occurred in children below 4 years old (with median age of 2), which was similar with the infection occurred in Spain 4,5 , France 3 and Singapore 6,28 . However, the median age was 6 years old in the HFMD outbreak caused by CV-A6 and CV-A10 of Finland 12 , and it was initially ignored because of most cases were adult infections, which suggested that a low herd immunity to these viruses. Most of reports about CV-A10 in China occurred mainly in spring and summer 10,29-32 , like the EV infection, but the 84% of HMFD cases in Finland were reported in October -December (autumn and winter).
The first strain of CV-A10 was isolated in New York in 1950 33 , prototype of CV-A10, Kowalik strain. Since then, CV-A10 was uncommon virus in Europe and United States 34-38 until several outbreaks was reported in Europe during 2008-2010. China national notifiable disease reported system showed that the proportion of other EV associated with HFMD was increasing in recent years 22 . In this study, we further identified the CV-A10 was the major pathogen of other EV responsible for HFMD in China. The number of CV-A10 circulating in China was becoming more and more prevalent (more patients and more provinces involved) and had a trend of increasing severe cases proportion from 2012-2016 (Table 1 and Fig. 1). This is a very important pathogenic information to guide the EV71 vaccination [39][40][41] and HFMD control and treatment.
Seven genotypes were assigned based on entire VP1 in this study, which was consistent with other studies 18,20,42 except for two new genotypes F and G. In addition, we rebuilt a phylogenetic tree on partial VP1 to investigate CV-A10 virus's origin and spread in the world, and to verify the genotype classification based on method of   complete VP1. In general, the genotyping results were consistent with each other, and two potential genotypes H and I were identified in separate group according to partial VP1. The comprehensive search from available public database showed H and I genotype viruses were obtained only from Japan, Finland, Germany and Egypt before 2007, and no more viruses were detected after that. Both of these genotypes might be disappeared after 2007, or they could not be detected because of surveillance gap. However, it is necessary to obtain the VP1 complete sequences of these suspicious genotype viruses to verify its genotyping classification and to better understand its spread in the world. As VP1 complete coding region contains many important neutralizing antigenic sites and there were generally consistent classification of genotyping between two sequencing window, we recommended VP1 to be as the standard target sequencing window for CV-A10. According to both complete VP1 and partial VP1 sequences, we found the similar situation: different genotype viruses have different geographic distribution pattern. Genotype B and G were mainly circulating in Chinese mainland and its surrounding areas, such as Taiwan China. While most of Chinese viruses, European viruses and Indian surrounding viruses, which solely clustered in genotype C, D, and F, respectively, showing geographic distribution circumscribed characterization. This is similar to other EV, like EV-A71 1 , but different with respiratory viruses such as measles 43,44 and human respiratory syncytial viruses 45,46 , which spread globally and different genotypes distributed cross over with each other in phylogenetic tree, not much geographic distribution circumscribed characterization.
There was an obvious genotype shift occurred in CV-A10 epidemic during 2004-2016 in mainland of China. Genotype B had been primary genotype circulating during 2004-2009, following by genotype C becoming predominant genotype during 2010-2016. It might not rule out the undetected genotype C before 2008, because of HFMD surveillance was initiated. Later on, based on the extensive and continuous HFMD surveillance, we concluded genotype B had disappeared completely. Similarly, sub-genotype C1 was mainly circulating during 2008-2012 and limited in few areas, but C2 gradually replaced C1 and became the predominant sub-genotype especially after 2010, which was similar to the alternation between C4a and C4b of EV-71 in China 15,17 . As the spread of sub-genotype C2 viruses, it has been evoluted and divided into 3 clusters. Viruses that isolated between 2008~2013 usually located in the cluster 1 and cluster 2, and 2014-2016 viruses were clustered in cluster 3. Comparing with cluster 1 and 2 viruses, cluster 3 viruses distributed in much more province with more prevalence during 2014 and 2016 (Fig. 4).
According to the mutations analysis among these sub-genotypes and clusters, we found that C2-cluster3 had three significant mutations and resulted in two nonsynonymous mutations and one synonymous mutation in VP1 region. Studies have shown that VP1 plays a critical role in mediating binding receptor, these mutations might cause some changes on the structure and function of VP1 protein, which might make the viruses having stronger transmissibility, infectivity and virulence.
The average evolutionary divergence of genotype C sequences between different year revealed that the nucleotide divergences within individual year from 2008-2016 were increasing gradually (Table 3), and this indicated that CV-A10 evolved constantly over years, its gene polymorphism gradually increased in these years. Based on Bayesian theory 17 , the scale of the pandemic was associated with the gene polymorphism: as the gene polymorphism increased, so do the scale of the pandemic. This is consistent with the findings in this study.
Most diseases associated with CV-A10 were mild and self-limiting, but the data from Singapore HFMD outbreak of in 2008 showed: comparing with EV-71, the virulence of CV-A10 and CV-A6 might be lower, but their transmissibility was stronger 6 . In the study of Cao et al., CV-A10 was demonstrated to be associated with severe complications defined by the same criteria, although with less effect than EV-71 7 . In this study, we found that the proportion of severe cases caused by CV-A10 had been gradually increasing since 2012. In addition, nearly half of severe viruses (17/33) located in cluster3.
All these evidence indicated C2-cluster3 viruses might have been evoluted with stronger transmissibility and virulence during the fitness with the immunity of host after persistent circulation in population. It is very necessary and very important to conduct the continuous and extensive virological surveillance for CV-A10, and further evolutionary studies, to better understand its evolution, transmissibility and virulence, will provide more evidence based scientific data to guide disease control and treatment.
One step RT-PCR was performed using the One Step PT-PCR kit Ver.2 (TaKaRa, #RR057A). The reaction system as followed: 3 μl of viral RNA was added to reaction mixture (total volume 25 μl), containing 12.5 μl2 × Buffer, 7.5 μldeionized H 2 O, 1 μl Enzyme Mix, and 0.5 μl each of the specific primers. PCR profile were 50 °C for 30 min, then 94 °C for 3 min; 32 cycles at 94 °C for 30 s, 50 °C for 30 s and 72 °C for 1 min and 20 s; and final extension step at 72 °C for 10 min. The products were analyzed by 1.5% agarose gel electrophoresis, and positive products were purified using the QIAquick Gel extraction kit (Qiagen, Valencia, CA). All the amplicons were sequenced by using both upper and down primers on an ABI Prism 3100 genetic analyzer 29 Table S2). Phylogenetic analysis using neighbor joining (NJ) and maximum likelihood (ML) was performed. A phylogenetic tree was constructed with the Kimura-2 parameter evolutionary models, and the reliability of it was tested by 1000 bootstrap replicates. Bootstrap values greater than 80% were considered statistically significant for grouping. Nucleotide Accession Number. The entireVP1 nucleotide sequence that represent different year and genotype/cluster of this study were deposited in the GenBank database under the accession number KF999730-KF999786 and MG838781-MG838884. (Supplementary Table S3).
Ethics Statement. This study did not involve human participants or human experimentation. Only specimen (stool samples, throat swab samples) collected from HFMD patients for public health purposes at the urging of the Ministry of Health, P. R. of China. Written informed consent for the use of their clinical samples was obtained from the parents of the children whose samples were analyzed. This study was approved by the second session of the Ethics Review Committee of the National Institute for Viral Disease Control and Prevention (NIVDC), Chinese Center for Disease Control and Prevention, all experimental protocols were approved by NIVDC, and the methods were carried out in accordance with the approved guidelines.