Characteristics of CVA16 B1c strains isolated for the first time in the Heilongjiang Province of China

Hand, foot, and mouth disease (HFMD) represents a globally prevalent infectious disease that is caused by enteroviruses. Enterovirus A71 (EV-A71), coxsackievirus A16 (CVA16), and coxsackievirus A6 (CVA6) are recognized as the predominant causative agents of HFMD. CVA16 is a member of the genus Enterovirus within the family Picornaviridae. B1a and B1b are the most prevalent subgenotypes, whereas the B1c subgenotype is relatively scarce. In this study, a comprehensive analysis was conducted on 15 of CVA16 B1c strains isolated from samples of patients diagnosed with HFMD in Jixi (Heilongjiang Province, China) in 2022. Subsequently, whole genome sequencing of these strains was carried out. Phylogenetic origin and potential recombination events were analyzed by aligning sequences of isolated of CVA16 B1c strains with related sequences deposited in GenBank. The CVA16 B1c isolates examined in this study exhibited a high degree of similarities. Specifically, the nucleotide similarity within the VP1 region ranged from 99.6 to 100%. The average nucleotide substitution rate of CVA16 B1c viruses worldwide was estimated to be 5.14 × 10⁻³ (4.13–6.27 × 10⁻³) substitution/site/year, and the most recent common ancestor could be traced back to 2003. The earliest CVA16 B1c strain isolated in China was traced back to 2011. Transmission pathway analysis suggested that Chinese strains may have originated in India. Recombination analysis showed that CVA16 B1c strains likely undergone recombination events with EV-A71 and CVA4. In conclusion, the analysis of a cluster of CVA16 B1c cases detected for the first time in Heilongjiang Province not only expanded the gene sequence library of CVA16 B1c strains in China but also offered an epidemiological basis for further investigations into the antigen–antibody interactions and pathogenicity of CVA16 B1c.

1 Introduction

Hand, foot and mouth disease (HFMD) is an infectious disease caused by various enteroviruses. It typically presents with mouth ulcers and rash on the hands and feet (Zhu et al., 2023). In 1957, an illness caused by coxsackievirus (mainly CVA16) was first reported in Toronto. Fever and mouth ulcers were characteristic manifestations of HFMD (Robinson et al., 1958). Previous studies have identified EV-A71 and CVA16 as the major causative pathogens of HFMD, which caused large-scale outbreaks in the United States and Asia–Pacific region (Zhu et al., 2023; González-Sanz et al., 2019; Noisumdaeng et al., 2019; Liu et al., 2000; Bible et al., 2008). CVA16 was first isolated in South Africa in 1951. According to the phylogenetic tree and genetic diversity of the VP1 gene, CVA16 strains can be divided into three genotypes: A, B, and D (Hu et al., 2021). Gene group B viruses can be further divided into B1a, B1b, B1c, and B2 clades (Zhang et al., 2010a). B2 was the dominant clade before 2000, but B1 has gradually replaced B2 as the most common clade after 1997 (Perera et al., 2007). Genotype D was first identified in Peru in 2009 and then spread to France from 2011 to 2014 (Hassel et al., 2017). In 2016, the SH-HP-16-51 strain was isolated in a patient with mild HFMD in Shanghai, China, and tested positive for the D gene (Wang et al., 2018a).

In 2005, the CVA16 B1c subtype was isolated in Malaysia for the first time. In recent years, CVA16 B1c has been detected successively in Malaysia (Chen et al., 2013; Perera et al., 2007; Chan et al., 2012), Japan (Mizuta et al., 2013), India (Ganorkar et al., 2017; Mamidi et al., 2024; Rao et al., 2017; Palani et al., 2016), and other countries. Among them, the HFMD study that broke out in the Andaman Islands of India in 2013 showed that CVA16 was the main pathogenic factor and was closely related to the B1c gene group (Palani et al., 2016). In China, CVA16 B1a and B1b strains co-circulated and spread continuously (Liu et al., 2000; Zhang et al., 2010a; Han et al., 2020). However, B1c was only detected in Shanghai in 2014 (Wang et al., 2018b), Xinjiang in 2017, Beijing in 2016 and 2022 (Hu et al., 2021; Liang et al., 2024), Guangdong in 2018 and 2023 (Xie et al., 2020; Zeng et al., 2025), lacking systematic epidemiological analysis. In this study, CVA16 B1c isolates detected in Jixi (Heilongjiang Province, China) in 2022 were analyzed using molecular epidemiology approaches. The evolution and propagation origins of B1c strains in China were explored using comparisons with similar sequences in GenBank. These results enrich the CVA16 B1c gene database in China and provide scientific data for further prevention and control of infections caused by these viruses.

2 Materials and methods

2.1 Ethics statement

All throat swab samples collected in this study were sent to the national laboratory after being processed at provincial and municipal levels. All specimens were collected after informed consent was obtained. The experimental protocol was approved by the Second Ethics Review Committee of the National Institute for Viral Disease Control and Prevention at the Chinese Center for Disease Control and Prevention.

2.2 Sample collection and virus isolation

In 2008, China established a national HFMD pathogen surveillance network. Representative samples from each province are sent to the national HFMD laboratory for verification. Clinical cases of HFMD are diagnosed according to the Guidelines for the Diagnosis and Treatment of HFMD (2018 edition) (Li et al., 2018). In this study, isolates from throat swab samples of children with mild HFMD, originally sent from the Heilongjiang Province in 2022, were used. Samples were collected from Jixi (n = 46), Harbin (n = 5), Jiamusi (n = 6), Yichun (n = 3), Qiqihar (n = 12), and Mudanjiang (n = 1). The samples were treated and inoculated into human rhabdomyosarcoma cells. The infected viral suspension was harvested after a complete cytopathic effect was observed. The samples were transported to the National Polio Laboratory for further isolation and identification within 24 h according to a rigorous cold chain procedure. The VP1 sequence was amplified using the universal EV primers. The genotype was determined using a criterion of VP1 nucleotide similarity >75% (amino acid similarity >85%) (Oberste et al., 1999; Oberste and Pallansch, 2005). MEGA7.0 software (v7.0.25; Sudhir Kumar, Arizona State University, Tempe, AZ, USA) (Kumar et al., 2016) was used to compare VP1 sequences, construct the phylogenetic tree through the “Models” option and reference sequence, and set the bootstrap value to 1,000 for the tree test. The gene subtype was determined based on the aggregation of the reference sequences in the evolutionary tree.

2.3 CVA16 B1c whole genome sequencing

The complete genome sequence was obtained using a “primer walk” strategy. 0001S48 was used as the upstream primer to amplify the start 5′ part of the sequence, and 7500A was used as the downstream primer to amplify the end 3′ part of the sequence (Yang et al., 2003). Sequencher software and online Oligonucleotide Properties Calculator¹ were used to design primers for segmented amplification. Some primers were taken from a previous study by Han et al. (2020) (Supplementary Table S1). A QIAGEN One-Step RT-PCR kit was used (QIAGEN, Hilden, Germany). Amplification was performed using the following steps: reverse transcription (50°C, 30 min, 94°C, 2 min), sequence amplification (94°C, 30°C, 50°C, 30°C, 72°C, 1 min 20 s) for 40 cycles, treatment at 72°C for 7 min, and preservation at 4°C. PCR products were identified using 1% agarose gel electrophoresis and sent to Tsingke Biotechnology Co., Ltd. (Beijing, China) for sequence determination. The sequencing results were collated and spliced using Sequencher (v5.4.5) to obtain the final complete genome sequence.

2.4 Analysis of VP1 sequences in CVA16 strains from Jixi City

MEGA7.0 software (v7.0.25) (Kumar et al., 2016) was used to verify the results of CVA16 genotyping. Nucleotide and amino acid genetic distances were calculated by “DISTANCE” program. The results of CVA16 genotyping were verified according to the division principle of 15–25% nucleotide difference of different genotypes within the same serotype proposed by Brown et al. (1999). Nucleotide and amino acid similarities of different gene subtypes of Jixi CVA16 strains were calculated using Sequence Identity Matrix program in BioEdit (v7.0.9.0) (Hall, 1999). ESPript3.0 online website was then used to analyze different Jixi CVA16 VP1 gene subtype area locus mutations (Robert and Gouet, 2014). The PROVEAN program was used to evaluate whether the variations had an impact on the function of the VP1 protein. Use −2.5 as the reference threshold (a score ≤ − 2.5 is considered harmful compilation, and a score > − 2.5 is considered neutral variation) (Choi et al., 2012).

2.5 Phylogenetic analysis

To elucidate the global evolutionary origin of CVA16 B1c strains, complete CVA16 VP1 sequences were obtained from GenBank and screened using TempEst software. Outliers were identified according to the vertical distance from the root to the node and residual, using the “root-to-tip” option to check the correlation (Rambaut et al., 2016). The “Models” package in MEGA software (v7.0.25) was used to select the optimal nucleotide replacement model. The molecular clock model was analyzed using a strict and an uncorrelated, relaxed clock combined with constant velocity, exponential growth, and logarithmic growth as the priority population growth modes. The Markov chain Monte Carlo method implemented in BEAST software (v1.10.4) (Suchard et al., 2018) was used, the chain length was set to 500,000,000, and the sampling frequency was set to 50,000. Convergence was tested using Tracer (v1.7.1) with an effective sample size > 200 as the criterion (Rambaut et al., 2018). A Bayesian maximum clade credibility (MCC) tree was constructed using TreeAnnotator (v1.10.4), and the top 10% of the sample trees were removed as burn-in. Phylogenetic trees were visualized using FigTree (v1.4.4). To analyze the dynamics of population diversity, Bayesian skyline plots were constructed based on the complete VP1 sequences.

2.6 Evolutionary analysis

To determine the source of CVA16 B1c strains in China, the global dynamic propagation of CVA16 B1c viruses was analyzed. An asymmetric alternative model with the Bayesian stochastic search variable selection option was selected in BEAST to infer the spread rate of the possible regions of different nodes and to reconstruct the ancestral origin of all strains. SpreadD3 (Bielejec et al., 2016) was used to calculate migration paths and Bayesian factor values (BF) between different regions. The BF value was used as the criterion to determine possible paths. BF > 100 indicated extremely strong evidence support, 10 < BF ≤ 100 indicated strong evidence support, and 3 < BF ≤ 10 indicated moderate evidence support. The numbers of Markov jumps and rewards used for state transitions between different regions were used to describe the migration states entering and leaving a region (Minin and Suchard, 2008).

2.7 Recombination analysis

The whole genome sequences of the selected CVA16 B1c strains divided into the 5′-UTR, P1, P2, and P3 regions were compared with those of the CVA16 prototype strain. The EV-A prototype strain sequence was obtained from GenBank, and phylogenetic trees of CVA16 B1c and EV-A prototype strains were constructed based on the 5’-UTR, P1, P2, and P3 regions. BLAST was used to screen for other group A sequences with more than 85% similarity to each region of CVA16 B1c. Sequences from India, Malaysia, France, and China were randomly selected, along with one sequence from this study. Potential Recombination event analysis was conducted using SimPlot software (version 3.5.1) and Recombination Detection Program (RDP4, version 4.100). Seven methods from RDP4 were adopted, including RDP, GENECONV, 3Seq, Chimaera, SiScan, MaxChi and BootScan. Recombination detected by at least three of the seven methods was judged to be possible, with p ≤ 0.05 (Martin et al., 2015). The Kimura’s two-parameter model was adopted in Simplot for the nucleotide substitution model, and the ratio of nucleotide transition and transversion was set to 2.0 (Lole et al., 1999). The recombinant results were verified using BioEdit by comparing the similarity between the reference sequence, recombinant donor, and prototype strain in different segments.

3 Results

3.1 Pathogen spectrum composition in HFMD samples isolated in the Heilongjiang Province in 2022

The VP1 sequence was determined in 73 positive specimens from six cities in the Heilongjiang Province. There were significant differences in the pathogen spectrum composition (Figures 1A,B). CVA16 was mainly detected in Jixi, where most specimens were collected. The highest pathogen diversity was detected in Qiqihar. CVA16, CVA6, CVA4, and CVA10 strains were isolated there, with CVA6 strains being the predominant pathogens. In the present study, isolates detected in Harbin and Jiamusi were represented by a single subtype, CVA16 and CVA6, respectively. The lowest number of samples was sent from Mudanjiang, with only one CVA16 positive case. Overall, CVA16 strains were isolated from the largest number of detected cases, accounting for 60% (44 cases, mainly CVA16 B1a). CVA6 strains were second most prevalent (30%; Figure 1C). Further analysis of CVA16 positive samples from Jixi showed that 15 B1c strains were detected (Figure 2). This accounted for 43% of CVA16 samples from Jixi (Figure 1D). For more details, please refer to Supplementary Table S2.

Figure 1

Figure 1. Spectrum of pathogenic coxsackievirus strains from HFMD samples isolated in the Heilongjiang Province (China) in 2022. (A, B) are graphs of the number of enterovirus pathogens detected in six regions of Heilongjiang Province. (C) represents the detected quantities of different enteroviruses and different subtypes of CVA16 in Heilongjiang Province; (D) represents the detected quantities of different subtypes of CVA16 in Jixi City.

Figure 2

Figure 2. Phylogenetic tree of CVA16 strains from Jixi (Heilongjiang Province, China).

3.2 Analysis of VP1 gene characteristics of CVA16 strains from Jixi

The genetic distances between all genotypes were calculated and compared. The inter-genotype genetic distances were greater than 15%, whereas the intra-genotype distances were less than 15%. The genetic distances within each subtype of B genotype were 0.047–0.075. The genetic distances between the subtypes were 0.101–0.145, in the range of 8–15%. The nucleotide similarities of Jixi CVA16 B1c isolates were 99.6–100%. The nucleotide similarities between Jixi CVA16 B1c and Heilongjiang Province CVA16 B1a and CVA16 B1b were 86–88.5%. Four characteristic loci of amino acid mutations were revealed in CVA16 B1c and other clades in Jixi strains: P3S, V25A, I235V, and T240A (Supplementary Figure S1). The PROVEAN results showed that the PROVEAN scores of these four variations were 1.458, 0.064, −0.063 and 0.304 respectively, all belonging to neutral mutations.

3.3 Evolutionary dynamics of CVA16 B1c

In total, 120 CVA16 B1c sequences from 2005 to 2022 were screened by TempEst (Supplementary Table S3), and an uncorrelated relaxed clock and logarithmic prior by marginal likelihood value were selected as the best molecular clock model for CVA16 B1c to construct the MCC tree. The average base substitution rate of CVA16 B1c worldwide was estimated at 5.14 × 10⁻³ substitution/site/year with a 95% highest posterior density (HPD) of 4.13 × 10⁻³– 6.27 × 10⁻³ substitution/site/year. The most recent common ancestor (MRCA) of global CVA16 B1c was dated to 2003 (95% HPD: 2000–2004; Figure 3A). All Chinese strains clustered with Indian strains. Chinese strains from 2022 were closely related to the Indian strains from 2018, with an approximate earliest origin in 2021. Chinese strains isolated from 2013 to 2017 were closely related to the Indian strains isolated in 2013. The earliest strain, CVA16 B1c, isolated in China, was traced back to 2011 and was a Guangdong strain. No other related sequences were obtained during the same period, and this was presumed to be an imported case. All French and Japanese strains isolated in 2011 were most closely related to the Indian strains from 2010. It is likely that the Malaysian strains clustered into a single clade because they originated earlier. Bayesian skyline analysis showed that the effective population size changed over time (Figure 3B). It remained relatively stable until 2011, after which three fluctuations were observed. The first wave occurred in 2011–2012, with a small increase in the population. The second wave witnessed a dynamic change of decline, increase and then decline from 2015 to 2018, and then tended to stabilize. In the third wave of 2021, the population showed a downward trend again.

Figure 3

Figure 3. Evolutionary analyses of CVA16 B1c strains. (A) Maximum Clade Credibility tree built using the Markov chain Monte Carlo method for the global VP1 coding region. Branches are colored according to the country of isolation. (B) Bayesian skyline plot based on VP1 sequences of CVA16 B1c strains reflects the effective population size in 2004–2022. (C) Histograms based on Markov jumps and rewards in different countries.

3.4 Analysis of the geographical spread of CVA16 B1c

Markov jumps and rewards analysis showed that India, Russia and Malaysia were mainly export-oriented, whereas other countries were mainly import-oriented (Figure 3C). To analyze the global dynamic propagation of CVA16 B1c further, the propagation network of CVA16 B1c strains was reconstructed. Geographical transmission analysis was performed for CVA16 B1c sequences from seven countries (14 from China, one from France, 88 from India, seven from Japan, one from Kazakhstan, seven from Malaysia, and two from Russia). Eight migration paths supported by BFs were identified (Figure 4). BF > 6,000 from India to China provided extremely strong support to the notion that Chinese CVA16 B1c strains were imported from India. The main migration paths could be: (I) from India to China, (II) from Russia to Kazakhstan, (III) from India to Japan, (IV) from Japan to France, (V) from Malaysia to Russia, (VI) from Malaysia to India, (VII) from Russia to India, and (VII) from France to Japan.

Figure 4

Figure 4. Analysis of the spatio-temporal dynamic propagation of CVA16 B1c strains. Only migration pathways with Bayesian factor (BF) values > 3 are shown. BF > 100 indicates extremely strong support, 10 < BF ≤ 100 indicates strong support, and 3 < BF ≤ 10 indicates moderate support.

3.5 Recombination analysis

The evolutionary trees of CVA16 5′-UTR, P1, P2, and P3 regions were constructed. The P1 region phylogenetic tree showed that all CVA16 B1c strains clustered together with CVA16 prototype strains (Figure 5B), which was consistent with the CVA16 serotyped results. However, phylogenetic trees based on the 5′-UTR, P2, and P3 regions showed different results. The CVA16 B1c strains clustered with CVA4 in the 5′-UTR region, and with EV-A71 in the P2 region (Figures 5A,C). The P3 region showed close similarity with EV-A71, CVA3, CVA12, CVA2, CVA6, and CVA10 (Figure 5D). Through RDP4 and Simplot analysis, it was found that the possible recombinant donors in the 5′-UTR region and the non-coding region were mainly CVA4 and EV-A71 (Figure 6 and Supplementary Figure S2). The similarity between different regions showed that the P1 region had the highest similarity with the CVA16 prototype strain. It has the highest similarity with EV-A71 in regions P2 and P3. The similarity with CVA4 was the highest in the 5’-UTR region (Supplementary Table S4).

Figure 5

Figure 5. Phylogenetic trees based on the 5′-UTR, P1, P2, and P3 regions of CVA16 B1c genomes and prototype EV-A strains. (A) 5′-UTR; (B) P1 region; (C) P2 region; (D) P3 region. Notes: Triangles represent CVA16 B1c sequences. Black square represents the prototype strain.

Figure 6

Figure 6. Bootscanning analysis of the whole CVA16 B1c strain genomes. The nucleotide transition and transversion ratios were set to 2.0, and the analysis was based on a 200-nucleotide sliding window with a 20-nucleotide step size.

4 Discussion

In this study, four characteristic loci were found in VP1 sequences of Jixi CVA16 B1c strains and strains of other clades: P3S, V25A, I235V, and T240A. Similarly, specific mutations at these sites were also found in the study of CVA16 B1c in Beijing and Guangdong (Hu et al., 2021; Liang et al., 2024; Zeng et al., 2025; Zhou et al., 2020). The VP1 region contains the major antigenic sites and most serotype-specific neutralization sites, and the substitution of amino acids may affect the function of proteins (Mateu, 1995). According to the prediction results of PROVEAN, P3S, V25A, I235V, and T240A are all neutral mutations. This may not cause a change in protein function, but it may also have an impact when it occurs simultaneously with the substitution of other regions. However, these mutations only exist in CVA16 B1c, and the same amino acid mutations have not been found in other genotypes of CVA16. At present, the impact of amino acid changes in B1c of CVA16 on the structure remains blank, and the biological significance reflected by such amino acid changes should be further evaluated. Whether these amino acid changes are related to virus antigenicity and pathogenicity needs further studies.

In the evolutionary analysis, the average nucleotide substitution rate of CVA16 B1c was estimated to be 5.14 × 10⁻³ substitution/site/year. Other studies on the evolutionary analysis of CVA16 have shown that the evolutionary rate of CVA16 is 3.36–6.66 × 10⁻³ substitution/site/year (Han et al., 2020; Zhao et al., 2016). It can be seen that CVA16 B1c does not show significant differences in terms of external evolutionary pressure. The analysis of global CVA16 B1c transmission suggested that Chinese B1c strains may have been imported from India. This conclusion is supported by a high Bayesian factor values. In this study, 15 CVA16 B1c strains were detected in the Heilongjiang Province of China, indicating a small outbreak. Apart from this study, no large-scale prevalence of B1c strains has been found in China, whereas long-term circulation of B1c strains have been reported in both India and Malaysia (Ganorkar et al., 2017; Chen et al., 2013; Chan et al., 2012). All CVA16 B1c strains in China were closely related to Indian strains, which further suggests the potential import risk of endemic strains in China and neighboring countries. Therefore, continuous and effective monitoring of HFMD worldwide should be emphasized to prevent the outbreak of CVA16 B1c in the future.

The symptoms of HFMD caused by CVA16 are generally mild: few severe cases and nervous system-related complications have been reported, and most of them were self-limiting (Zhu et al., 2010). However, some studies have shown that CVA16 can cause more serious complications and clinical symptoms (Wright et al., 1963; Wang et al., 2004), especially upon co-infection with EV-A71 (Chang et al., 1999; Cai et al., 2019). Recombination is one of the main mechanisms of enterovirus evolution. It was found in this study that CVA16 may have recombined with EV-A71 and CVA4. Other studies on CVA16 from Taiwan, Guangdong, central and northeastern China have suggested that the 5 ‘-UTR and non-structural protein coding regions of CVA16 subgenotypes B1 and B2 may recombine with CVA4 and EV-A71, respectively (Cheng et al., 2022; Han et al., 2014; Zhao et al., 2011; Lin et al., 2015). Previous studies have confirmed that a large outbreak of HFMD with fatal neurological complications in Fuyang City, China, in 2008 was mainly caused by the C4 sub-genotype of EV-A71, which had recombined with CVA16 in the 3D region (Zhang et al., 2010b). EV recombination may affect virulence and pathogenicity of the virus. Therefore, further attention should be paid to the recombination of CVA16 and EV-A71 strains.

Since the monovalent EV-A71 vaccine was launched in China the end of 2015, the detection rate of EV-A71 has decreased significantly. However, the lack of cross-protection against other EVs causing HFMD led to the detection rates of CVA6, CVA10 and CVA16 significantly exceeding that of EV-A71 (Zhang et al., 2022; Hong et al., 2022). At present, there is no specific vaccine for CVA16 to prevent HFMD on the market. Studies have indicated that there are EV-related multivalent vaccines in the clinical trial stage. These multivalent vaccines contain CVA16 and are developed based on the previously prevalent B1a and B1b subtypes (Zeng et al., 2025). CVA16 B1c has been occasionally detected in some provinces in recent years, such as Shanghai (Wang et al., 2018b), Beijing (Hu et al., 2021; Liang et al., 2024), and Guangdong (Zeng et al., 2025). Due to the antigenic differences among subgenotypes, it is necessary to develop an effective multivalent vaccine with cross-protective immune response. In addition to vaccines, disease prevention and control also requires further strengthening of comprehensive monitoring of EVs. The Asia-Pacific Network for Enterovirus Surveillance (APNES) was established in 2017 to assess the burden of disease through laboratory diagnosis and data collection (Chiu et al., 2020). The European Union has also established the European non-polio enterovirus network (ENPEN) to improve the diagnosis of enteroviruses, monitor enterovirus data and prevent the spread of enteroviruses (Harvala et al., 2018). For China, merely conducting HFMD monitoring cannot detect all enteroviruses in a timely manner. To supplement its limitations, a network of encephalitis and meningitis monitoring laboratories was established based on sentinel hospitals. At the same time, patients with acute flaccid paralysis (AFP) were monitored, and environmental surveillance of sewage water was carried out. Strengthen prevention and control measures to reduce the risk of epidemic outbreak.

In conclusion, this study is the first detection of CVA16 B1c in Heilongjiang province. Combined with the available CVA16 B1c data in the GenBank database, the isolated CVA16 B1c in Jixi City was analyzed. The global phylogeny and spread of CVA16 B1c were deeply understood through the Bayesian method, and recombination analysis was carried out. This study enriched the global genetic database of CVA16 B1c and provided an epidemiological basis for its further research.

Data availability statement

The accession numbers NMDCN0007NAQ - NMDCN0007NB8 of the National Microbiology Data Center (NMDC).

Ethics statement

The studies involving humans were approved by the Second Ethics Review Committee of the National Institute for Viral Disease Control and Prevention at the Chinese Center for Disease Control and Prevention. The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent for participation in this study was provided by the participants’ legal guardians/next of kin. Written informed consent was obtained from the minor(s)’ legal guardian/next of kin for the publication of any potentially identifiable images or data included in this article.

Author contributions

TY: Writing – original draft, Formal analysis, Data curation, Writing – review & editing. QS: Formal analysis, Data curation, Writing – review & editing, Supervision. HLi: Data curation, Writing – review & editing, Formal analysis. MY: Resources, Writing – review & editing. SC: Writing – review & editing, Resources. WS: Resources, Writing – review & editing. JXu: Writing – review & editing, Resources. DY: Supervision, Writing – review & editing, Methodology. SZ: Writing – review & editing, Methodology, Supervision. ZH: Methodology, Data curation, Resources, Writing – review & editing. JXi: Methodology, Data curation, Writing – review & editing, Formal analysis. HLu: Writing – review & editing, Methodology, Data curation. YL: Writing – review & editing, Data curation, Methodology. QY: Project administration, Supervision, Writing – review & editing, Methodology, Funding acquisition, Data curation, Formal analysis. YZ: Data curation, Methodology, Supervision, Project administration, Writing – review & editing, Funding acquisition, Formal analysis.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This study was supported by National Key Laboratory of Intelligent Tracking and Forecasting for Infectious Diseases, National Institute for Viral Disease Control and Prevention, Chinese Center for Disease Control and Prevention (Project Nos. ZDGWNLJS25-36) and Health and Family Planning Commission of Heilongjiang province (Project Nos. 2018561).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The authors declare that no Gen AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2025.1634547/full#supplementary-material

Footnotes

1. ^https://www.northwestern.edu/

References

Bible, J. M., Iturriza-Gomara, M., Megson, B., Brown, D., Pantelidis, P., Earl, P., et al. (2008). Molecular epidemiology of human enterovirus 71 in the United Kingdom from 1998 to 2006. J. Clin. Microbiol. 46, 3192–3200. doi: 10.1128/JCM.00628-08

PubMed Abstract | Crossref Full Text | Google Scholar

Bielejec, F., Baele, G., Vrancken, B., Suchard, M. A., Rambaut, A., and Lemey, P. (2016). SpreaD3: interactive visualization of spatiotemporal history and trait evolutionary processes. Mol. Biol. Evol. 33, 2167–2169. doi: 10.1093/molbev/msw082

PubMed Abstract | Crossref Full Text | Google Scholar

Brown, B. A., Oberste, M. S., Alexander, J. P., Kennett, M. L., and Pallansch, M. A. (1999). Molecular epidemiology and evolution of enterovirus 71 strains isolated from 1970 to 1998. J. Virol. 73, 9969–9975. doi: 10.1128/JVI.73.12.9969-9975.1999

PubMed Abstract | Crossref Full Text | Google Scholar

Cai, K., Wang, Y., Guo, Z., Yu, H., Li, H., Zhang, L., et al. (2019). Clinical characteristics and managements of severe hand, foot and mouth disease caused by enterovirus A71 and coxsackievirus A16 in Shanghai, China. BMC Infect. Dis. 19:285. doi: 10.1186/s12879-019-3878-6

PubMed Abstract | Crossref Full Text | Google Scholar

Chan, Y. F., Wee, K. L., Chiam, C. W., Khor, C. S., Chan, S. Y., Wan Nor Amalina, W. M. Z., et al. (2012). Comparative genetic analysis of VP4, VP1 and 3D gene regions of enterovirus 71 and coxsackievirus A16 circulating in Malaysia between 1997-2008. Trop. Biomed. 29, 451–466.

Google Scholar

Chang, L. Y., LIN, T. Y., Huang, Y. C., TSAO, K. C., SHIH, S. R., KUO, M. L., et al. (1999). Comparison of enterovirus 71 and coxsackievirus A16 clinical illnesses during the Taiwan enterovirus epidemic, 1998. Pediatr. Infect. Dis. J. 18, 1092–1096. doi: 10.1097/00006454-199912000-00013

Crossref Full Text | Google Scholar

Chen, X., Tan, X., Li, J., Jin, Y., Gong, L., Hong, M., et al. (2013). Molecular epidemiology of Coxsackievirus A16: Intratype and prevalent intertype recombination identified. PLoS One 8:e82861. doi: 10.1371/journal.pone.0082861

PubMed Abstract | Crossref Full Text | Google Scholar

Cheng, D., Chiu, Y. W., Huang, S. W., Lien, Y. Y., Chang, C. L., Tsai, H. P., et al. (2022). Genetic and cross neutralization analyses of Coxsackievirus A16 circulating in Taiwan from 1998 to 2021 suggest dominant genotype B1 can serve as vaccine candidate. Viruses 14:2306. doi: 10.3390/v14102306

PubMed Abstract | Crossref Full Text | Google Scholar

Chiu, M. L., Luo, S. T., Chen, Y. Y., Chung, W. Y., Duong, V., Dussart, P., et al. (2020). Establishment of Asia-Pacific network for enterovirus surveillance. Vaccine 38, 1–9. doi: 10.1016/j.vaccine.2019.09.111

PubMed Abstract | Crossref Full Text | Google Scholar

Choi, Y., Sims, G. E., Murphy, S., Miller, J. R., and Chan, A. P. (2012). Predicting the functional effect of amino acid substitutions and indels. PLoS One 7:e46688. doi: 10.1371/journal.pone.0046688

PubMed Abstract | Crossref Full Text | Google Scholar

Ganorkar, N. N., Patil, P. R., Tikute, S. S., and Gopalkrishna, V. (2017). Genetic characterization of enterovirus strains identified in hand, foot and mouth disease (HFMD): emergence of B1c, C1 subgenotypes, E2 sublineage of CVA16, EV71 and CVA6 strains in India. Infect. Genet. Evol. 54, 192–199. doi: 10.1016/j.meegid.2017.05.024

PubMed Abstract | Crossref Full Text | Google Scholar

González-Sanz, R., Casas-Alba, D., Launes, C., Muñoz-Almagro, C., Ruiz-García, M. M., Alonso, M., et al. (2019). Molecular epidemiology of an enterovirus A71 outbreak associated with severe neurological disease, Spain, 2016. Euro Surveill. 24:1800089. doi: 10.2807/1560-7917.ES.2019.24.7.1800089

PubMed Abstract | Crossref Full Text | Google Scholar

Hall, T. A. (1999). Bioedit: a user-friendly biological sequence alignment editor and analysis program for windows 95/98/NT. Nuclc Acids Symp. Ser. 41, 95–98. doi: 10.1021/bk-1999-0734.ch008

Crossref Full Text | Google Scholar

Han, Z., Song, Y., Xiao, J., Jiang, L., Huang, W., Wei, H., et al. (2020). Genomic epidemiology of coxsackievirus A16 in mainland of China, 2000-18. Virus Evol. 6:veaa084. doi: 10.1093/ve/veaa084

PubMed Abstract | Crossref Full Text | Google Scholar

Han, J. F., Yu, N., Pan, Y. X., He, S. J., Xu, L. J., Cao, R. Y., et al. (2014). Phenotypic and genomic characterization of human coxsackievirus A16 strains with distinct virulence in mice. Virus Res. 179, 212–219. doi: 10.1016/j.virusres.2013.10.020

PubMed Abstract | Crossref Full Text | Google Scholar

Harvala, H., Broberg, E., Benschop, K., Berginc, N., Ladhani, S., Susi, P., et al. (2018). Recommendations for enterovirus diagnostics and characterisation within and beyond Europe. J. Clin. Virol. 101, 11–17. doi: 10.1016/j.jcv.2018.01.008

PubMed Abstract | Crossref Full Text | Google Scholar

Hassel, C., Mirand, A., Farkas, A., Diedrich, S., Huemer, H. P., Peigue-Lafeuille, H., et al. (2017). Phylogeography of Coxsackievirus A16 reveals global transmission pathways and recent emergence and spread of a recombinant Genogroup. J. Virol. 91, e00630–e00617. doi: 10.1128/JVI.00630-17

PubMed Abstract | Crossref Full Text | Google Scholar

Hong, J., Liu, F., Qi, H., Tu, W., Ward, M. P., Ren, M., et al. (2022). Changing epidemiology of hand, foot, and mouth disease in China, 2013−2019: a population-based study. Lancet Reg. Health West. Pac. 20:100370. doi: 10.1016/j.lanwpc.2021.100370

PubMed Abstract | Crossref Full Text | Google Scholar

Hu, Y. F., Jia, L. P., Yu, F. Y., Liu, L. Y., Song, Q. W., Dong, H. J., et al. (2021). Molecular epidemiology of coxsackievirus A16 circulating in children in Beijing, China from 2010 to 2019. World J. Pediatr. 17, 508–516. doi: 10.1007/s12519-021-00451-y

PubMed Abstract | Crossref Full Text | Google Scholar

Kumar, S., Stecher, G., and Tamura, K. (2016). MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol. Biol. Evol. 33, 1870–1874. doi: 10.1093/molbev/msw054

PubMed Abstract | Crossref Full Text | Google Scholar

Li, X. W., Ni, X., Qian, S. Y., Wang, Q., Jiang, R. M., Xu, W. B., et al. (2018). Chinese guidelines for the diagnosis and treatment of hand, foot and mouth disease (2018 edition). World J. Pediatr. 14, 437–447. doi: 10.1007/s12519-018-0189-8

PubMed Abstract | Crossref Full Text | Google Scholar

Liang, Z., Lin, C., Huo, D., Yang, Y., Feng, Z., Cui, S., et al. (2024). First detection of multiple cases related to CV-A16 strain of B1c clade in Beijing in 2022. J. Med. Virol. 96:e29796. doi: 10.1002/jmv.29796

PubMed Abstract | Crossref Full Text | Google Scholar

Lin, C. H., Wang, Y. B., Chen, S. H., Hsiung, C. A., and Lin, C. Y. (2015). Precise genotyping and recombination detection of enterovirus. BMC Genomics 16:S8. doi: 10.1186/1471-2164-16-S12-S8

PubMed Abstract | Crossref Full Text | Google Scholar

Liu, C. C., Tseng, H. W., Wang, S. M., Wang, J. R., and Su, I. J. (2000). An outbreak of enterovirus 71 infection in Taiwan, 1998: epidemiologic and clinical manifestations. J. Clin. Virol. 17, 23–30. doi: 10.1016/S1386-6532(00)00068-8

PubMed Abstract | Crossref Full Text | Google Scholar

Lole, K. S., Bollinger, R. C., Paranjape, R. S., Gadkari, D., Kulkarni, S. S., Novak, N. G., et al. (1999). Full-length human immunodeficiency virus type 1 genomes from subtype C-infected Seroconverters in India, with evidence of Intersubtype recombination. J. Virol. 73, 152–160. doi: 10.1128/JVI.73.1.152-160.1999

PubMed Abstract | Crossref Full Text | Google Scholar

Mamidi, P., Mishra, B., Panda, S., Ray, A., Jena, D., Biswas, V. K., et al. (2024). Identification of a novel Indian clade of Coxsackievirus A16 in hand, foot and mouth disease outbreak in 2022. Future Virol. 19, 161–172. doi: 10.1080/17460794.2024.2350247

Crossref Full Text | Google Scholar

Martin, D. P., Murrell, B., Golden, M., Khoosal, A., and Muhire, B. (2015). RDP4: detection and analysis of recombination patterns in virus genomes. Virus Evol. 1:vev003. doi: 10.1093/ve/vev003

PubMed Abstract | Crossref Full Text | Google Scholar

Mateu, M. G. (1995). Antibody recognition of picornaviruses and escape from neutralization: a structural view. Virus Res. 38, 1–24. doi: 10.1016/0168-1702(95)00048-u

PubMed Abstract | Crossref Full Text | Google Scholar

Minin, V. N., and Suchard, M. A. (2008). Counting labeled transitions in continuous-time Markov models of evolution. J. Math. Biol. 56, 391–412. doi: 10.1007/s00285-007-0120-8

PubMed Abstract | Crossref Full Text | Google Scholar

Mizuta, K., Abiko, C., Aoki, Y., Ikeda, T., Matsuzaki, Y., Hongo, S., et al. (2013). Molecular epidemiology of Coxsackievirus A16 strains isolated from children in Yamagata, Japan between 1988 and 2011. Microbiol. Immunol. 57, 400–405. doi: 10.1111/1348-0421.12041

PubMed Abstract | Crossref Full Text | Google Scholar

Noisumdaeng, P., Korkusol, A., Prasertsopon, J., Sangsiriwut, K., Chokephaibulkit, K., Mungaomklang, A., et al. (2019). Longitudinal study on enterovirus A71 and coxsackievirus A16 genotype/subgenotype replacements in hand, foot and mouth disease patients in Thailand, 2000–2017. Int. J. Infect. Dis. 80, 84–91. doi: 10.1016/j.ijid.2018.12.020

PubMed Abstract | Crossref Full Text | Google Scholar

Oberste, M. S., Maher, K., Kilpatrick, D. R., and Pallansch, M. A. (1999). Molecular evolution of the human enteroviruses: correlation of serotype with VP1 sequence and application to picornavirus classification. J. Virol. 73, 1941–1948. doi: 10.1128/JVI.73.3.1941-1948.1999

PubMed Abstract | Crossref Full Text | Google Scholar

Oberste, M. S., and Pallansch, M. A. (2005). Enterovirus molecular detection and typing. Rev. Res. Med. Microbiol. 16, 163–171. doi: 10.1097/01.revmedmi.0000184741.90926.35

Crossref Full Text | Google Scholar

Palani, S., Nagarajan, M., Biswas, A. K., Maile, A., and Paluru, V. (2016). B1c genetic subtype of coxsackievirus A16 associated with hand, foot and mouth disease in Andaman Islands, India. Trans. R. Soc. Trop. Med. Hyg. 110, 421–423. doi: 10.1093/trstmh/trw042

PubMed Abstract | Crossref Full Text | Google Scholar

Perera, D., Yusof, M. A., Podin, Y., Ooi, M. H., Thao, N. T. T., Wong, K. K., et al. (2007). Molecular phylogeny of modern coxsackievirus A16. Arch. Virol. 152, 1201–1208. doi: 10.1007/s00705-006-0934-5

PubMed Abstract | Crossref Full Text | Google Scholar

Rambaut, A., Drummond, A. J., Xie, D., Baele, G., and Suchard, M. A. (2018). Posterior summarization in Bayesian Phylogenetics using tracer 1.7. Syst. Biol. 67, 901–904. doi: 10.1093/sysbio/syy032

PubMed Abstract | Crossref Full Text | Google Scholar

Rambaut, A., Lam, T. T., Carvalho, L. M., and Pybus, O. G. (2016). Exploring the temporal structure of heterochronous sequences using TempEst (formerly path-O-gen). Virus Evol. 2:vew007. doi: 10.1093/ve/vew007

PubMed Abstract | Crossref Full Text | Google Scholar

Rao, D. C., Naidu, J. R., Maiya, P. P., Babu, A., and Bailly, J. L. (2017). Large-scale HFMD epidemics caused by Coxsackievirus A16 in Bangalore, India during 2013 and 2015. Infect. Genet. Evol. 55, 228–235. doi: 10.1016/j.meegid.2017.08.030

PubMed Abstract | Crossref Full Text | Google Scholar

Robert, X., and Gouet, P. (2014). Deciphering key features in protein structures with the new ENDscript server. Nucleic Acids Res. 42, W320–W324. doi: 10.1093/nar/gku316

PubMed Abstract | Crossref Full Text | Google Scholar

Robinson, C. R., Doane, F. W., and Rhodes, A. J. (1958). Report of an outbreak of febrile illness with pharyngeal lesions and exanthem: Toronto, summer 1957; isolation of group A Coxsackie virus. Can. Med. Assoc. J. 79, 615–621.

Google Scholar

Suchard, M. A., Lemey, P., Baele, G., Ayres, D. L., Drummond, A. J., and Rambaut, A. (2018). Bayesian phylogenetic and phylodynamic data integration using BEAST 1.10. Virus Evol. 4:vey016. doi: 10.1093/ve/vey016

PubMed Abstract | Crossref Full Text | Google Scholar

Wang, C. Y., Li, L. F., Wu, M. H., Lee, C. Y., and Huang, L. M. (2004). Fatal coxsackievirus A16 infection. Pediatr. Infect. Dis. J. 23, 275–276. doi: 10.1097/01.inf.0000115950.63906.78

Crossref Full Text | Google Scholar

Wang, J., Teng, Z., Chu, W., Fang, F., Cui, X., Guo, X., et al. (2018a). The emergence and spread of one Coxsackievirus A16 Genogroup D novel recombinant strain that caused a clustering HFMD outbreak in Shanghai, China, 2016. Emerg. Microbes Infect. 7, 1–3. doi: 10.1038/s41426-018-0134-x

PubMed Abstract | Crossref Full Text | Google Scholar

Wang, J., Teng, Z., Cui, X., Li, C., Pan, H., Zheng, Y., et al. (2018b). Epidemiological and serological surveillance of hand-foot-and-mouth disease in Shanghai, China, 2012-2016. Emerg. Microbes Infect. 7:8. doi: 10.1038/s41426-017-0011-z

PubMed Abstract | Crossref Full Text | Google Scholar

Wright, H. T., Landing, B. H., Lennette, E. H., and McAllister, R. M. (1963). Fatal infection in an infant associated with Coxsackie virus group a, type 16. N. Engl. J. Med. 268, 1041–1044. doi: 10.1056/NEJM196305092681904

PubMed Abstract | Crossref Full Text | Google Scholar

Xie, J., Yang, X. H., Hu, S. Q., Zhan, W. L., Zhang, C. B., Liu, H., et al. (2020). Co-circulation of coxsackieviruses A-6, A-10, and A-16 causes hand, foot, and mouth disease in Guangzhou city, China. BMC Infect. Dis. 20:271. doi: 10.1186/s12879-020-04992-x

PubMed Abstract | Crossref Full Text | Google Scholar

Yang, C. F., Naguib, T., Yang, S. J., Nasr, E., Jorba, J., Ahmed, N., et al. (2003). Circulation of endemic type 2 vaccine-derived poliovirus in Egypt from 1983 to 1993. J. Virol. 77, 8366–8377. doi: 10.1128/JVI.77.15.8366-8377.2003

PubMed Abstract | Crossref Full Text | Google Scholar

Zeng, H., Zeng, B., Yi, L., Qu, L., Cao, J., Yang, F., et al. (2025). The emergence of Coxsackievirus A16 subgenotype B1c: a key driver of the hand, foot, and mouth disease epidemic in Guangdong, China. Viruses 17:219. doi: 10.3390/v17020219

PubMed Abstract | Crossref Full Text | Google Scholar

Zhang, Y., Wang, D., Yan, D., Zhu, S., Liu, J., Wang, H., et al. (2010a). Molecular evidence of persistent epidemic and evolution of subgenotype B1 Coxsackievirus A16-associated hand, foot, and mouth disease in China. J. Clin. Microbiol. 48, 619–622. doi: 10.1128/JCM.02338-09

PubMed Abstract | Crossref Full Text | Google Scholar

Zhang, X., Zhang, Y., Li, H., and Liu, L. (2022). Hand-foot-and-mouth disease-associated enterovirus and the development of multivalent HFMD vaccines. Int. J. Mol. Sci. 24:169. doi: 10.3390/ijms24010169

PubMed Abstract | Crossref Full Text | Google Scholar

Zhang, Y., Zhu, Z., Yang, W., Ren, J., Tan, X., Wang, Y., et al. (2010b). An emerging recombinant human enterovirus 71 responsible for the 2008 outbreak of hand foot and mouth disease in Fuyang city of China. Virol. J. 7:94. doi: 10.1186/1743-422X-7-94

PubMed Abstract | Crossref Full Text | Google Scholar

Zhao, K., Han, X., Wang, G., Hu, W., Zhang, W., and Yu, X. F. (2011). Circulating Coxsackievirus A16 identified as recombinant type A human enterovirus, China. Emerg. Infect. Dis. 17, 1537–1540. doi: 10.3201/eid1708.101719

PubMed Abstract | Crossref Full Text | Google Scholar

Zhao, G., Zhang, X., Wang, C., Wang, G., and Li, F. (2016). Characterization of VP1 sequence of Coxsackievirus A16 isolates by Bayesian evolutionary method. Virol. J. 13:130. doi: 10.1186/s12985-016-0578-3

PubMed Abstract | Crossref Full Text | Google Scholar

Zhou, Y., Tan, L. V., Luo, K., Van Tan, L., Liao, Q., Wang, L., et al. (2020). Genetic variation of multiple serotypes of enteroviruses associated with hand, foot and mouth disease in southern China. Virol. Sin. 36:61. doi: 10.1007/s12250-020-00266-7

PubMed Abstract | Crossref Full Text | Google Scholar

Zhu, P., Ji, W., Li, D., Li, Z., Chen, Y., Dai, B., et al. (2023). Current status of hand-foot-and-mouth disease. J. Biomed. Sci. 30:15. doi: 10.1186/s12929-023-00908-4

PubMed Abstract | Crossref Full Text | Google Scholar

Zhu, Z., Zhu, S., Guo, X., Wang, J., Wang, D., Yan, D., et al. (2010). Retrospective seroepidemiology indicated that human enterovirus 71 and coxsackievirus A16 circulated wildly in central and southern China before large-scale outbreaks from 2008. Virol. J. 7:300. doi: 10.1186/1743-422X-7-300

PubMed Abstract | Crossref Full Text | Google Scholar