HFMD is classified as a legally reportable infectious disease in China. Standardized case questionnaires are used to collect demographic and epidemiologic data, which are reported in a timely manner to the web-based China Disease Prevention and Control Information System [2,28]. For the purposes of this study, these aforementioned data were extracted from this national database. In our analysis, an “epidemic season” for HFMD means that the number of cases reported in that month constituted at least 15% of the total cases reported throughout the year. The “epidemic peak” was identified as the month witnessing the highest influx of reported cases.

The baseline data of HFMD cases in Henan Children’s Hospital (Children’s Hospital Affiliated to Zhengzhou University) and the Third Affiliated Hospital of Zhengzhou University (Henan Maternal and Child Health Hospital of Henan Province), which had advantages in the diagnosis and treatment of HFMD, were also included in this study. From both hospitals, the pathogenetic data of HFMD cases from 2021 to 2023 were selected as representation. Biological samples including fecal samples and throat swabs were collected from HFMD cases for the study.

This study received ethical approval from the Committee for Ethical Review of Zhengzhou University (ZZUIRB2023-180), and written informed consent was obtained from all participants’ parents. The obtained data have been anonymized or deidentified.

Viral RNA was extracted from fecal samples of the abovementioned HFMD cases. Each fecal sample was stored at −80 °C; the sample was thawed at room temperature, and 0.1 g of the sample was placed in a 1.5-mL sterile, enzyme-free, grinding tube. Next, 1 mL of TRIzol reagent (Thermo Fisher Scientific) along with 2 sterile, enzyme-free, steel beads (3 mm) were added to the tube. The tube was then placed in a precooled tissue grinder and grinded 3 times for 30 seconds, each time at a frequency of 60 Hz. Chloroform was then added and the mixture was centrifuged, retaining the supernatant. Isopropanol was then added and allowed to stand for 10 minutes before being placed in the centrifuge again and discarding the supernatant. Then, 1 mL of 75% ethanol was added and mixed thoroughly, and another centrifugation was initiated. At the final stage, the supernatant was discarded, and the precipitate was dried and resuspended in 30 μL of RNase-free water. The concentration of RNA was measured using a NanoDrop spectrophotometer (Thermo Fisher Scientific).

Viral RNA was converted to complementary DNA using the Hifair II First Strand cDNA Synthesis Kit (Yeasen Biotechnology) according to the instructions. The VP1 section of the CV-A6 sequence was amplified using primers (forward primer: AAGGACACTGACGAGATTCAGC; reverse primer: AGTGGCGAGATGTCGGTTT; 1026 base pairs). Then, the amplified fragments were sent to Beijing Liuhe Huada Gene Technology Co., Ltd. for sequencing.

For the study, we retrieved the full-length sequence of CV-A6 from the GenBank database using “Coxsackievirus A6” as an index term with a cutoff date of January 1, 2024. From the GenBank database, all CV-A6 (near) whole-genome sequences (sequence length of 6600-7500 base pairs) were extracted. Any disreputable or subpar sequences were excluded. This includes laboratory-adapted strains, clones, strains with elevated passage numbers, sequences containing numerous unascertained bases, and mislabeled sequences that belonged to other serotypes. Ultimately, 858 near full-length sequences were gathered from GenBank. We constructed a phylogenetic evolutionary tree using the full-length CV-A6 sequences of 858 different strains to understand the phylogeny of CV-A6 strains in China. All data sets are listed in Multimedia Appendix 1. Next, we selected representative strains from different countries and years to understand the phylogeny based on the CV-A6 VP1 isolates and excluded inferior-quality sequences. The representative strains should cover more provinces of China to ensure a high geographic representation and cover each cluster of the phylogenetic trees to acquire a high phylogenic representation. After strict quality control, we compiled a data set of 313 VP1 sequences for phylogenetic and phylodynamic analyses. The samples cover 16 countries worldwide, and the time series were from 1970 to 2023. All data sets are listed in Multimedia Appendix 2. The complete VP1 nucleotide sequences of CV-A6 was 915 nucleotides long, which encode 305 amino acids.

We used MEGA 7.0 software [29] to compare the VP1 nucleotide sequence of CV-A6. Then, we calculated a maximum likelihood (ML) phylogenetic tree using IQ-TREE (version 1.6.12) [30] under the best-fit nucleotide substitution model as designated by Model Finder. For CV-A6, the SYM+G4 model was deemed the best fit, as outlined in Table 1 [31]. We used the ultrafast bootstrap method with 1000 replicates for our analyses. Analysis of temporal molecular evolutionary signals for the data set was conducted using TempEst (version 1.5) [32]. In brief, regression analyses were used to determine the relationship between sampling dates and root-to-tip genetic divergence obtained from the ML phylogeny. The slope of the regression line provides an estimate of the rate of evolution in substitutions per site per year, and the intercept with the time-axis constitutes an estimate of the age of the root. The time-scaled tree along with historical population dynamics were estimated using the BEAST package (version 1.10.4) [33]. Information regarding the time of the most recent common ancestor (tMRCA) and the rate of evolution (represented as substitutions per site per year) were derived from this estimation. To reconstruct the phylogeny, we used a Bayesian MCMC approach [33]. The nucleotide substitution model was tested through the use of Model Finder, with the Bayesian Sky grid model featuring a relaxed molecular clock being used. The MCMC lengths of 50,000,000 generations were completed with individual samplings convened every 5000 generations. Output from BEAST was analyzed within the Tracer program (version 1.7.1) [34] to ensure convergence through graphical checks and adequate quality control parameters of the posterior distribution (effective sample size>200). The maximum clade credibility (MCC) tree was ascertained through the use of Tree Annotator (version 1.10.4) and was depicted in FigTree (version 1.4.4) [35] in Multimedia Appendix 3. Then, the tMRCA and its 95% highest posterior density (HPD) were quantified in chronological years.

This study further analyzed the differences between the 34 CV-A6 strains and the prototype CV-A6 strain. The nucleotide sequences and deduced amino acid sequences of the above 34 strains were compared with the amino acid sequence of the CV-A6 prototype strain. Amino acid sites with mutations or deletions were analyzed using BioEdit software (version 7.2.5.0) [36].

From 2021 to 2023, a total of 85,060 HFMD cases were included in our study (n=23,440, 27.56% in 2021; n=26,914, 31.64% in 2022; and n=34,706, 40.8% in 2023). Our data showed that the epidemic season in Henan usually lasted from June to August, with peaks around June and July (Figure 1A). During this period, the monthly case reporting rate ranged from 20.7%(4854/23,440) to 35% (12,135/34,706) of the total annual number of cases. Next, we analyzed the pathogen composition of 2850 laboratory-confirmed cases from Henan Children’s Hospital and the Third Affiliated Hospital of Zhengzhou University. The analysis identified 8 EV types, but 1145 EV-positive specimens (40.84%, 1164/2850) could not be typed. Among them, CV-A6 was the most frequently identified EV type, accounting for 652 (22.88%) out of 2850 cases, followed by CV-A16 (n=624, 21.89%), EV-A71 (n=326, 11.44%), and CV-A10 (n=84, 2.95%). In 2021, CV-A16 was the main pathogen of HFMD with 23.48% (331/1410) of the cases. CV-A6 emerged as the primary pathogen for HFMD in 2022 and 2023, constituting 27.73% (203/732) and 37.01% (262/708) of cases, respectively. Most importantly, the detection rate of CV-A6 spiked from 13.26% (187/1410) in 2021 to 37.01% (262/708) in 2023 (Table 2). Out of all identified EV pathogen types, CV-A6 had the highest composition ratio at 22.88% (652/2850), and the composition ratios of other pathogen types are depicted in Figure 1B.

To analyze the genomic characteristics and patterns of genetic evolution of isolated strains, we downloaded the full-length sequence of CV-A6 from the NCBI database and performed phylogenetic analysis on the strains isolated in this study. Genogroups A, B, and C contained very few members due to the lack of complete VP1 nucleotide sequence data and limited virus circulation. The virulent strain isolated in this study was subtype D3 (Figure 2A). In this study, the full-length sequences in China were analyzed by subgenotype and year. It was found that since 2013, the prevalence of subgenotype D1 viruses had gradually declined, and a large number of D2 strains began to appear, accompanied by a small number of D3 subtypes. Since 2017, the D3 subtype has begun to proliferate, and by 2020, all CV-A6 strains in China were of the D3 subtype (Figure 2B).

In this study, a total of 514 full-length CV-A6 sequences from various regions within China were downloaded from the GenBank database. Only 1 CV-A6 isolate was obtained before 2008, and 44 isolates were reported between 2008 and 2012; strains isolated between 2013 and 2023 accounted for 91.24% (469/514; Figure 3A). These 514 isolates provide a representative coverage of 7 different Chinese geographic regions: East China (comprising Shandong, Jiangsu, Zhejiang, Fujian, Anhui, Shanghai, and Taiwan), South China (Guangdong, Guangxi, Hainan, Hong Kong, and Macao), Central China (Henan, Hubei, Hunan, and Jiangxi), North China (Inner Mongolia, Tianjin, Beijing, Hebei, and Shanxi), Northwest China (Xinjiang, Qinghai, Shaanxi, Ningxia, and Gansu), Southwest China (Sichuan, Yunnan, Guizhou, and Chongqing), and Northeast China (Heilongjiang, Jilin, and Liaoning; Figure 3B).

To gain insights into the evolutionary patterns of the isolated CV-A6 strains in our study and their relationships with previous CV-A6 strains and other viruses from different regions, we assembled a data set comprising 313 VP1 coding sequences. Using this data set, we constructed an initial ML tree for root-to-tip analysis. This analysis revealed a strong positive temporal signal with an R² value of 1, indicating a clear evolutionary trend over time (Figure 4A). We developed an MCC tree that estimated the divergence times of CV-A6 based on the 313 full-length VP1 sequences. The MCC tree exhibited consistency with the phylogenetic tree in terms of topology and typing results (Figure 4B). Furthermore, the mean substitution rate and tMRCA were generally calculated by temporal phylogenies for the evolution of CV-A6 full-length VP1 sequences. The results indicated a mean substitution rate of 4.5456×10⁻³ (95% HPD 3.8466×10⁻³ to 5.3796×10⁻³) substitution per site per year. The predicted tMRCA for the EV-A71 prototype strain was approximately 1824 (95% HPD 1702-1902). The ancestor of global CV-A6 isolates emerged approximately 108 years ago, dating the origin of CV-A6 to around 1915 (95% HPD 1870-1948). The specific details of the tMRCA of the genetic subtypes (A-D) of the CV-A6 strains, along with their respective 95% HPD ranges, are shown in Table 3. Furthermore, the earliest estimated time of CV-A6 transmission in China, as determined through BEAST analysis, dates back to 1989 (95% HPD 985-1991). This estimation aligns with the node date of the CV-A6 Chinese cluster as displayed in the MCC tree. In addition, our analysis predicted that the tMRCA for the isolated CV-A6 strains could trace back to 2019 (95% HPD 2018-2020).

The 305 amino acids encoded in the VP1 gene of 34 CV-A6 isolates were compared with Gdula (AY421764.1) using MegAlign software (DNASTAR) to determine the amino acid mutation levels. In this study, we excluded strains with identical amino acid mutation sites and retained 17 strains with different mutation sites (Figure 5A). Amino acid substitutions or deletions occurred at all 11 of these sites (loci 5, 8, 10, 14, 32, 98, 160, 174, 194, 261, and 279). In addition, some strains had amino acid substitutions at locus 137. Unresolved questions about CV-A6 mutations need to be investigated to understand their roles in CV-A6 evolution and to enrich the information on relevant mutations in EVs. Further monitoring and analysis of mutated loci and epidemiological investigations are needed to determine whether mutations occurring at these loci are associated with changes in virus transmission, pathogenicity, and virulence. From the analysis of evolutionary differences between the CV-A6 strains isolated in this study and the prototype CV-A6 strain, we found that the genetic differences between the isolates were small, while the differences between the isolates and the prototype strain (AY421764.1) were large, and the values were all greater than 0.199 (Figure 5B). The greater the genetic distance, the higher the probability of mutation and recombination in the virus and the CV-A6 virus evolving over time.