Published on in Vol 10 (2024)

Preprints (earlier versions) of this paper are available at https://preprints.jmir.org/preprint/59604, first published .
Molecular Evolutionary Dynamics of Coxsackievirus A6 Causing Hand, Foot, and Mouth Disease From 2021 to 2023 in China: Genomic Epidemiology Study

Molecular Evolutionary Dynamics of Coxsackievirus A6 Causing Hand, Foot, and Mouth Disease From 2021 to 2023 in China: Genomic Epidemiology Study

Molecular Evolutionary Dynamics of Coxsackievirus A6 Causing Hand, Foot, and Mouth Disease From 2021 to 2023 in China: Genomic Epidemiology Study

1Department of Infectious Diseases, Children’s Hospital Affiliated to Zhengzhou University, , Zhengzhou, , China

2College of Public Health, Zhengzhou University, , Zhengzhou, , China

3Henan International Joint Laboratory of Children’s Infectious Diseases, Children's Hospital Affiliated to Zhengzhou University, , Zhengzhou, , China

4NHC Key Laboratory of Birth Defects Prevention, Henan Key Laboratory of Population Defects Prevention, , Zhengzhou, , China

Corresponding Author:

Yuefei Jin, PhD


Background: Hand, foot, and mouth disease (HFMD) is a global public health concern, notably within the Asia-Pacific region. Recently, the primary pathogen causing HFMD outbreaks across numerous countries, including China, is coxsackievirus (CV) A6, one of the most prevalent enteroviruses in the world. It is a new variant that has undergone genetic recombination and evolution, which might not only induce modifications in the clinical manifestations of HFMD but also heighten its pathogenicity because of nucleotide mutation accumulation.

Objective: The study assessed the epidemiological characteristics of HFMD in China and characterized the molecular epidemiology of the major pathogen (CV-A6) causing HFMD. We attempted to establish the association between disease progression and viral genetic evolution through a molecular epidemiological study.

Methods: Surveillance data from the Chinese Center for Disease Control and Prevention from 2021 to 2023 were used to analyze the epidemiological seasons and peaks of HFMD in Henan, China, and capture the results of HFMD pathogen typing. We analyzed the evolutionary characteristics of all full-length CV-A6 sequences in the NCBI database and the isolated sequences in Henan. To characterize the molecular evolution of CV-A6, time-scaled tree and historical population dynamics regarding CV-A6 sequences were estimated. Additionally, we analyzed the isolated strains for mutated or missing amino acid sites compared to the prototype CV-A6 strain.

Results: The 2021-2023 epidemic seasons for HFMD in Henan usually lasted from June to August, with peaks around June and July. The monthly case reporting rate during the peak period ranged from 20.7% (4854/23,440) to 35% (12,135/34,706) of the total annual number of cases. Analysis of the pathogen composition of 2850 laboratory-confirmed cases identified 8 enterovirus serotypes, among which CV-A6 accounted for the highest proportion (652/2850, 22.88%). CV-A6 emerged as the major pathogen for HFMD in 2022 (203/732, 27.73%) and 2023 (262/708, 37.01%). We analyzed all CV-A6 full-length sequences in the NCBI database and the evolutionary features of viruses isolated in Henan. In China, the D3 subtype gradually appeared from 2011, and by 2019, all CV-A6 virus strains belonged to the D3 subtype. The VP1 sequences analyzed in Henan showed that its subtypes were consistent with the national subtypes. Furthermore, we analyzed the molecular evolutionary features of CV-A6 using Bayesian phylogeny and found that the most recent common ancestor of CV-A6 D3 dates back to 2006 in China, earlier than the 2011 HFMD outbreak. Moreover, the strains isolated in 2023 had mutations at several amino acid sites compared to the original strain.

Conclusions: The CV-A6 virus may have been introduced and circulating covertly within China prior to the large-scale HFMD outbreak. Our laboratory testing data confirmed the fluctuation and periodic patterns of CV-A6 prevalence. Our study provides valuable insights into understanding the evolutionary dynamics of CV-A6.

JMIR Public Health Surveill 2024;10:e59604

doi:10.2196/59604

Keywords



Hand, foot, and mouth disease (HFMD) is an acute infectious disease that mainly affects children aged <5 years; its symptoms include fever and rashes on the hands, feet, oral mucosa, and buttocks [1]. In China, the annual peak season for HFMD outbreaks usually lasts from April to June. Some areas also experience a surge of HFMD from October to November, notably in Southern China [2,3]. HFMD has increasingly become a public health concern and has a substantial effect on the disease burden and economic status of the affected regions [4,5]. HFMD is primarily caused by enterovirus (EV) A, and the most common pathogens are EV-A71 and coxsackievirus (CV) A16 [6]. Recently, CV-A6 has gained global concern, with a large number of patients reported in Finland [7], Singapore [8], Japan [9], Spain [10], and Thailand [11]. Since the end of 2012, CV-A6 has emerged as the main pathogen of HFMD in China [12]. By 2013, Guangdong, Jilin, and Beijing had recorded cases of HFMD caused by CV-A6 [13-15]. However, most studies have focused on major serotypes such as EV-A71 and CV-A16, with a limited amount of information on other EV serotypes. This is largely due to these other EV serotypes generally causing lesser disseminated case. CV-A6 is one of the most prevalent EVs in the world in recent years, and it is a new variant that has undergone genetic recombination and evolution. It is specifically associated with genotype D [16]. The CV-A6 D3 subtype associated with HFMD was first described in 2008 following an outbreak in Finland and subsequently spread to other European countries [16,17]. After 2018, the D3 subtype emerged in Asia and was described from samples collected during outbreaks in Philippines [18], Vietnam [19], and other countries [20]. HFMD cases caused by CV-A6 tend to have clinical characteristics that differ from those of typical HFMD, such as atypical herpes [21], and can lead to severe neurological or cardiopulmonary complications, such as encephalitis, aseptic meningitis, encephalomyelitis, acute flaccid paralysis, and fatal cardiopulmonary complications [22-24]. Recombination and evolution could potentially alter not only the clinical manifestations of HFMD but also heighten the disease’s pathogenicity due to cumulative nucleotide mutations [25]. Most importantly, CV-A6 appears to influence a more extensive demographic of children, inducing a more severe and prolonged ailment compared to other EVs [26]. However, the evolution and transmission of CV-A6 in China have been relatively underexplored in recent years [27]. In this study, 34 HFMD-associated CV-A6 strains were isolated from Henan Children’s Hospital and the Third Affiliated Hospital of Zhengzhou University, which greatly enriched the global CV-A6 VP1 sequence database. We downloaded all CV-A6 full-length sequences before January 2024 from the NCBI database to characterize the evolution of the isolated CV-A6 strains and to analyze the temporal and geographic distribution of different CV-A6 subtypes using Bayesian Markov chain Monte Carlo (MCMC) methodology in China. Furthermore, this study characterized the amino acid mutations of the VP1 gene. This study aimed to provide insight into the molecular epidemiology of CV-A6 and establish a link between disease development and viral evolution.


Collection of Epidemiological Data

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.

Data Sources

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.

Ethical Consideration

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.

RNA Extraction

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 Sequencing

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.

Data Set Construction

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.

Phylodynamic Analysis

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.

Table 1. The best-fit nucleotide substitution model of coxsackievirus A6 determined by Model Finder.
No.ModelLogarithm of site likelihooddfAICaAICcbBICc
1JC19,581.40614939,460.81239,488.8239,460.812
2JC+G418,154.72115036,609.44136,637.84237,429.203
3TN+F18,121.215436,550.39936,580.40536,580.405
4TN+F+G416,528.70415533,367.40833,397.82334,214.496
5TNe18,120.7115136,543.42136,572.21937,368.648
6TNe+G416,535.9515233,375.89933,405.09734,206.592
7K2P18,178.31115036,656.62136,685.02237,476.383
8K2P+G416,555.57615133,413.15233,441.9534,238.38
9K2P18,178.31115036,656.62136,685.02237,476.383
10K2P+G416,555.57615133,413.15233,441.9534,238.38
11F81+F19,582.94715239,469.89339,499.09140,300.586
12F81+F+G418,149.83515336,605.67136,635.27137,441.828
13HKY+F18,166.57715336,639.15536,668.75537,475.312
14HKY+F+G416,544.79815433,397.59733,427.60334,239.219
15SYM18,096.17315436,500.34536,530.35237,341.968
16SYM+G416,515.59915533,341.19733,371.61234,188.285
17TIM+F18,119.78715536,549.57536,579.9937,396.663
18TIM+F+G416,528.38115633,368.76233,399.58934,221.315
19TVM+F18,132.14715636,576.29436,607.12137,428.847
20TVM+F+G416,512.71715733,339.43533,370.67634,197.453
21TVMe18,153.43515336,612.86936,642.4737,449.027
22TVMe+G416,536.23315433,380.46533,410.47234,222.088
23GTR+F18,086.98115736,487.96236,487.96237,345.98
24GTR+F+G4d16,495.34915833,306.698e33,338.357e34,170.181e

aAIC: Akaike information criterion.

bAICc: corrected Akaike information criterion.

cBIC: Bayesian information criterion.

dBest-fit model: GTR+F+G4 was chosen according to BIC.

eBest fit.

Base Substitution and Amino Acid Mutation Analysis

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].


The Epidemiology and Etiology of HFMD in Henan, China, from 2021 to 2023

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.

Figure 1. Pathogen spectrum of HFMD in Henan, China, from 2021 to 2023: (A) the monthly incidence of HFMD and (B) composition ratio of various EVs in Henan, from 2021 to 2023. CV: coxsackievirus; EV: enterovirus; HFMD: hand, foot, and mouth disease.
Table 2. Enterovirus (EV) serotypes of patients with hand, foot, and mouth disease (n=2850) in Henan, China, from 2021 to 2023.
YearEV-A71, n (%)CVa-A16, n (%)CV-A6, n (%)CV-A10, n (%)Other EV serotypes, n (%)Total, n (%)
2021 (n=1410)141 (10)331 (23.48)187 (13.26)67 (4.75)684 (48.51)1410 (100)
2022 (n=732)76 (10.38)180 (24.59)203 (27.73)4 (0.55)269 (36.75)732 (100)
2023 (n=708)109 (15.40)113 (15.96)262 (37.01)13 (1.83)211 (29.80)708 (100)
Total (n=2850)326 (11.44)624 (21.89)652 (22.88)84 (2.95)1164 (40.84)2850 (100)

aCV: coxsackievirus.

Subtyping Analysis of CV-A6

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).

Figure 2. Subtyping analysis of CV-A6. (A) Genetic evolutionary characteristics of enteroviruses isolated in Henan, China. Maximum likelihood phylogenetic trees for enteroviruses were constructed with CV-A6 full-length sequences. The red box marked the strains that have been isolated this study. The colors of the branches represent countries, and the colors on the strains represent subtypes. (B) Annual percentage distribution of subgenotypes D1, D2, and D3 associated with HFMD from 2009 to 2023 in mainland China. CV: coxsackievirus; HFMD: hand, foot, and mouth disease.
Figure 3. The temporal and geographical distribution characteristics of CV-A6 in China. (A) Annual and regional distribution of 514 CV-A6 strains in China. (B) The distribution of CV-A6 strains associated to HFMD in 7 significant regions of China, where both region and sequence quantity are denoted and differentiated by the provided color key. The colors on the map represent different geographical areas, and the colors on the pie charts represent the pathogens of HFMD. CV: coxsackievirus; HFMD: hand, foot, and mouth disease.

Molecular Evolution Analysis of CV-A6

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 R2 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−3 (95% HPD 3.8466×10−3 to 5.3796×10−3) 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).

Figure 4. Molecular evolution analysis of CV-A6. (A) Temporal signal analysis of root-to-tip divergence regression versus date (R2=1). (B) The phylogenetic tree was constructed with maximum clade credibility for 313 VP1 sequences of CV-A6 strains. The purple bars at the nodes indicate the 95% HPDs of tMRCAs. Red vertical lines represent the 34 CV-A6 strains isolated in this study. CV: coxsackievirus; HPD: highest posterior density; tMRCA: time of the most recent common ancestor.
Table 3. The mean and 95% HPDa values of the tMRCAb and evolutionary rate of coxsackievirus A6.
Genetic subtypestMRCA (year), mean (95% HPD)Evolutionary rate (×10–3 substitutions/site/year), mean (95% HPD)
Group A1940 (1929‐1951)3.31 (0.24‐11.66)
Group B1973 (1959‐1984)2.39 (0.46‐7.7)
Group C1978 (1971‐2003)2.81 (0.27‐10.11)
Group D11994 (1992‐1997)2.47 (0.15‐8.53)
Group D22003 (1998‐2007)6.46 (0.98‐3.18)
Group D32000 (1999‐2003)5.30 (0.80‐15.94)

aHPD: highest posterior density.

btMRCA: time of the most recent common ancestor.

Amino Acid Variation Site Analysis of CV-A6

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.

Figure 5. Amino acid mutation sites and genetic differences analysis of the VP1 gene. (A) Analysis of amino acid mutations in the VP1 gene of the isolated strains of CV-A6. (B) Estimates of evolutionary divergence between sequences of isolated CV-A6 strains and the prototype CV-A6 strain (AY421764.1). Values are indicated by color shading. Analyses were conducted using the maximum composite likelihood model. The analysis involved 35 nucleotide sequences. CV: coxsackievirus.

HFMD outbreaks have proven to be a substantial public health concern in the Asia-Pacific region, particularly in China, over the past decade [2,37]. Consequently, attention from researchers toward CV-A6 has been growing. The sporadic new cases of CV-A6 infections over the past decade have indicated a possible worldwide resurgence of EV pathogens [7,17,38,39]. In this study, the epidemic season in Henan usually lasted from June to August, with peaks around June and July, from 2021 to 2023. Next, we analyzed the pathogen composition of 2850 laboratory-confirmed cases in Henan Children’s Hospital and the Third Affiliated Hospital of Zhengzhou University. Among them, CV-A6 was the most important pathogen causing HFMD (652/2850, 22.88%). We found that the D3 subtype gradually appeared from 2011 onward, and all CV-A6 virus strains belonged to the D3 subtype in 2019. We analyzed the molecular evolutionary features of CV-A6 using Bayesian phylogeny and found that the tMRCA of D3 subtype of CV-A6 in China was inferred to be 2005 (95% HPD 2002‐2007), suggesting that the D3 subtype had been circulating in China for many years. In addition, the strain isolated in 2023 had mutations at several amino acid sites compared to the original strain. Some mutations in the evolution of EVs are important for their reproduction, pathogenesis, and transmission [40-42]. These findings can provide ideas for the research of the pathogenic mechanism and allow researchers to look for ways to inhibit viral replication or transmission of new targets.

We note that there has been a significant increase in HFMD outbreaks in Henan Province during the summer months, with peak infections occurring between June and July. This aligns with the epidemic seasons observed in countries such as Thailand and India [20,43]. The specific cause for this has not been clarified. However, climate and geographical location are thought to contribute to the seasonal patterns of HFMD to some degree. High temperature and high humidity have been linked to an increased incidence of HFMD. Records also indicate a correlation between atmospheric pressure, wind speed, and the occurrence of HFMD [44]. Previous epidemiological studies on CV-A6 have revealed a transition from the D2 to D3 branch of CV-A6 evolution occurring in China and France [45,46]. Since 2010, there have been large-scale HFMD outbreaks in several countries, including Brazil [47], China [48], India [49], Japan [50], and Vietnam [51], which have been attributed to CV-A6 D3 strains. Our study also indicates that these CV-A6 D3 strains share closer phylogenetic relationships with viruses that circulated in Thailand, France, Japan, and Vietnam from 2008 to 2021, suggesting the presence of a virus spreading over long distances. Analyzing the phylogenetic tree of the complete VP1 sequence, we found evidence of coevolution between CV-A6 strains prevalent in Henan and CV-A6 strains sourced from other parts of China. This indicates that there may be multiple transmission chains for the strains studied. Given the vast topography and climatic diversity across China, the transmission and replication capacity of CV-A6 may vary by region. It is also worth noting that CV-A6’s subgenotype D3 may be more transmissible, infectious, and virulent, which could potentially explain the sustained international circulation of CV-A6. From the 313 CV-A6 strains included in this study, we constructed a molecular clock evolutionary tree for the VP1 gene. On the basis of previous studies, this study divided CV-A6 into 4 subtypes (A-D), of which D was divided into D1, D2, and D3 subtypes, and described the common ancestor time of each subtype in detail by Bayesian analysis method, which was missing in previous studies [52]. Among them, the first isolation of the D3 subtype of CV-A6 in China occurred in 2011, and the tMRCA was inferred to be 2005 (95% HPD 2002‐2007), suggesting that D3 experienced frequent transmission in China during the 14 years of the epidemic. In addition, the tMRCA of the 2023 isolates was 2019, suggesting that the CV-A6 D3 subtype had been secretly circulating in China for many years since the occurrence of this strain in 2023. This also means that the virus is in the process of continuing to evolve and could lead to large outbreaks in immunity-naive populations. Here, CV-A6 was seen to diverge into 4 evolutionary branches, with the CV-A6 strains situated in the D branch showing more active evolution. This analysis of the evolutionary source implies the importance of bolstering monitoring efforts for strains in the D branch. These strains not only comprise the Henan strains but also include CV-A6 isolates from regions neighboring mainland China. These strains have notably evolved continuously over a relatively short period since 1990, spreading across other regions and forming numerous chains of transmission. The molecular clock evolutionary tree of CV-A6 gives rise to a hypothesis that during the evolutionary genetic process, there were 3 major sources of global CV-A6 transmission: Vietnam; Japan; and various provinces of China, such as Henan, Shandong, Yunnan, etc. Currently, CV-A6 endemic regions are predominantly found in China and its neighboring countries. Mutation and recombination are common processes in EVs, driving significant genotypic and phenotypic variability [53]. These genetic alterations, which can result in antigenic shifts, could potentially spark outbreaks. Nonsynonymous mutations in VP1, leading to changes in amino acids, are frequently found among EVs. Our study has found that the average rates of VP1 nucleotide substitutions fall between 3.8466×10−3 and 5.3796×10−3 substitutions per site per year in CV-A6, which is consistent with the evolutionary rate estimated from expanded CV-A6 surveys (>300 VP1 sequences) and rates estimated previously [54]. It is imperative to note that genetic alterations at crucial sites in VP1 could give rise to immune-escape mutants. These mutants could impact the efficacy of vaccines as they allow the virus to circumvent immune defenses, possibly leading to poorer control over the disease in populations. Amino acid variants A5T, V30A, S137N, and I242V were also present in the CV-A6 epidemic strain that caused HFMD in Vietnam from 2011 to 2015 [19]. Furthermore, the structure of 137N is a bend, and the substitution from S to N at the 137th position will introduce a new N-glycosylation site [55]; however, the function of this locus needs to be further investigated. The CV-A6 may be mutating toward increased pathogenicity, which may be one of the reasons for the increased prevalence of CV-A6 in Zhengzhou city. Besides, it has been suggested that amino acid mutation sites (96-98, 102, 106, 151, 160, 165, 174, and 216) of the CV-A6 strain may be associated with major neutralizing antigenic epitopes [56]. The amino acid mutation sites (N10S, Q98L, S194T, F261L, and S279T) in this study are also consistent with the epidemiology reported by laboratories in the Philippines [18]. Moreover, we identified previously unreported variants at loci 8, 14, and 32 and other locations in this study, which may be the main drivers of the variants in the VP1 gene of CV-A6. Some mutations in the evolution of EVs are important for their reproduction, pathogenesis, and transmission [40-42]. Furthermore, the emergence of new genotypes or subtypes, intragenotype and intergenotype substitutions, and the introduction of new variants arising from spontaneous mutations and genetic recombination all pose challenges for vaccine development, selection, and vaccination policies. There is a need for laboratory-based surveillance and epidemiological surveillance for genetic change and evolution studies.

A number of important limitations need to be considered regarding these findings. First, as with any common self-limited disease, HFMD surveillance captures only the tip of the clinical iceberg, and most cases go undetected because they are asymptomatic, or because patients do not seek formal care or are not diagnosed and reported. Second, because of the lack of real-time and effective pathogen typing in some areas, the typing of atypical HFMD was insufficient. We may have underestimated the actual depth of harm caused by HFMD. We cannot rule out this monitoring bias. Finally, due to geographical restrictions, this study only selected cases from 2 representative hospitals in Henan Province for virus isolation, with a small sample size. Follow-up studies can further increase the sample size analysis and fully understand the epidemic trend of HFMD.

In this study, we used 34 HFMD-associated CV-A6 strains isolated from Henan Children’s Hospital and The Third Affiliated Hospital of Zhengzhou University, significantly enriching the global CV-A6 VP1 sequence database. We explored the molecular epidemic characteristics of CV-A6 in mainland China and around the world, and we predicted the amino acid mutation sites associated with CV-A6, which deepened our understanding of CV-A6 and provided valuable information for the molecular epidemiology of CV-A6. This not only enhanced our understanding of CV-A6 but also provided valuable insights into its molecular epidemiology. As the pathogen spectrum of HFMD grows more complex, the proportion of previously dominant strains, namely CV-A16 and EV-A71, have decreased. At the same time, there were some rare, easily ignored, and vulnerable epidemic serotypes of EVs that are cocirculating. In light of these findings, a variety of emerging technologies should be used to improve sample detection and enhance surveillance and reporting of other non–EV-A71 and non–CV-A6 EVs. This will enable a more comprehensive understanding of global HFMD epidemiologic trends and provide stronger options for managing and reducing future outbreaks.

Acknowledgments

All phases of this study were supported by the National Natural Science Foundation of China (YJ; 82372229 and 82002147); Open Research Fund of National Health Commission Key Laboratory of Birth Defects Prevention & Henan Key Laboratory of Population Defects Prevention (YJ; ZD202301); Open Project of Henan Province Engineering Research Center of Diagnosis and Treatment of Pediatric Infection and Critical Care (YJ; ERC202302); and 2024 Graduate Independent Innovation Project of Zhengzhou University (YC; 20240331).

Data Availability

All data generated or analyzed during this study are included in the published paper.

Conflicts of Interest

None declared.

Multimedia Appendix 1

Coxsackievirus A6 sequence numbers for phylogenetic evolutionary analysis.

XLSX File, 22 KB

Multimedia Appendix 2

Coxsackievirus A6 VP1 gene sequence numbers for phylogenetic and phylodynamic analyses.

XLSX File, 15 KB

Multimedia Appendix 3

Maximum clade credibility tree depicted in FigTree.

ZIP File, 6138 KB

  1. Yu F, Zhu R, Jia L, et al. Sub-genotype change and recombination of coxsackievirus A6s may be the cause of it being the predominant pathogen for HFMD in children in Beijing, as revealed by analysis of complete genome sequences. Int J Infect Dis. Oct 2020;99:156-162. [CrossRef] [Medline]
  2. Xing W, Liao Q, Viboud C, et al. Hand, foot, and mouth disease in China, 2008-12: an epidemiological study. Lancet Infect Dis. Apr 2014;14(4):308-318. [CrossRef] [Medline]
  3. Rui J, Luo K, Chen Q, et al. Early warning of hand, foot, and mouth disease transmission: a modeling study in mainland, china. PLoS Negl Trop Dis. Mar 24, 2021;15(3):e0009233. [CrossRef] [Medline]
  4. Feng H, Duan G, Zhang R, Zhang W. Time series analysis of hand-foot-mouth disease hospitalization in Zhengzhou: establishment of forecasting models using climate variables as predictors. PLoS One. Jan 31, 2014;9(1):e87916. [CrossRef] [Medline]
  5. Lee BY, Wateska AR, Bailey RR, Tai JHY, Bacon KM, Smith KJ. Forecasting the economic value of an enterovirus 71 (EV71) vaccine. Vaccine. Nov 16, 2010;28(49):7731-7736. [CrossRef] [Medline]
  6. Lu QB, Zhang XA, Wo Y, et al. Circulation of coxsackievirus A10 and A6 in hand-foot-mouth disease in China, 2009-2011. PLoS One. Dec 18, 2012;7(12):e52073. [CrossRef] [Medline]
  7. Osterback R, Vuorinen T, Linna M, Susi P, Hyypiä T, Waris M. Coxsackievirus A6 and hand, foot, and mouth disease, Finland. Emerg Infect Dis. Sep 2009;15(9):1485-1488. [CrossRef] [Medline]
  8. Wu Y, Yeo A, Phoon MC, et al. The largest outbreak of hand; foot and mouth disease in Singapore in 2008: the role of enterovirus 71 and coxsackievirus A strains. Int J Infect Dis. Dec 2010;14(12):e1076-e1081. [CrossRef] [Medline]
  9. Fujimoto T, Iizuka S, Enomoto M, et al. Hand, foot, and mouth disease caused by coxsackievirus A6, Japan, 2011. Emerg Infect Dis. Feb 2012;18(2):337-339. [CrossRef] [Medline]
  10. Montes M, Artieda J, Piñeiro LD, Gastesi M, Diez-Nieves I, Cilla G. Hand, foot, and mouth disease outbreak and coxsackievirus A6, northern Spain, 2011. Emerg Infect Dis. Apr 2013;19(4):676-678. [CrossRef] [Medline]
  11. Puenpa J, Chieochansin T, Linsuwanon P, et al. Hand, foot, and mouth disease caused by coxsackievirus A6, Thailand, 2012. Emerg Infect Dis. Apr 2013;19(4):641-643. [CrossRef] [Medline]
  12. Feng X, Guan W, Guo Y, et al. A novel recombinant lineage’s contribution to the outbreak of coxsackievirus A6-associated hand, foot and mouth disease in Shanghai, China, 2012-2013. Sci Rep. Jun 30, 2015;5(1):11700. [CrossRef] [Medline]
  13. Han JF, Xu S, Zhang Y, et al. Hand, foot, and mouth disease outbreak caused by coxsackievirus A6, China, 2013. J Infect. Sep 2014;69(3):303-305. [CrossRef] [Medline]
  14. Hongyan G, Chengjie M, Qiaozhi Y, et al. Hand, foot and mouth disease caused by coxsackievirus A6, Beijing, 2013. Pediatr Infect Dis J. Dec 2014;33(12):1302-1303. [CrossRef] [Medline]
  15. Lu J, Zeng H, Zheng H, et al. Hand, foot and mouth disease in Guangdong, China, in 2013: new trends in the continuing epidemic. Clin Microbiol Infect. Jul 2014;20(7):O442-O445. [CrossRef] [Medline]
  16. Bian L, Wang Y, Yao X, Mao Q, Xu M, Liang Z. Coxsackievirus A6: a new emerging pathogen causing hand, foot and mouth disease outbreaks worldwide. Expert Rev Anti Infect Ther. Jun 25, 2015;13(9):1061-1071. [CrossRef] [Medline]
  17. Song Y, Zhang Y, Ji T, et al. Persistent circulation of coxsackievirus A6 of genotype D3 in mainland of China between 2008 and 2015. Sci Rep. Jul 14, 2017;7(1):5491. [CrossRef] [Medline]
  18. Foronda JLM, Jiao MMAD, Climacosa FMM, Oshitani H, Apostol LNG. Epidemiological and molecular characterization of coxsackievirus A6 causing hand, foot, and mouth disease in the Philippines, 2012-2017. Infect Genet Evol. Oct 2023;114:105498. [CrossRef] [Medline]
  19. Anh NT, Nhu LNT, Van HMT, et al. Emerging coxsackievirus A6 causing hand, foot and mouth disease, Vietnam. Emerg Infect Dis. Apr 2018;24(4):654-662. [CrossRef] [Medline]
  20. Gopalkrishna V, Ganorkar N. Epidemiological and molecular characteristics of circulating CVA16, CVA6 strains and genotype distribution in hand, foot and mouth disease cases in 2017 to 2018 from Western India. J Med Virol. Jun 2021;93(6):3572-3580. [CrossRef] [Medline]
  21. Aswathyraj S, Arunkumar G, Alidjinou EK, Hober D. Hand, foot and mouth disease (HFMD): emerging epidemiology and the need for a vaccine strategy. Med Microbiol Immunol. Oct 2016;205(5):397-407. [CrossRef] [Medline]
  22. Li K, Li X, Si W, Liang H, Xia HM, Xu Y. Identifying risk factors for neurological complications and monitoring long-term neurological sequelae: protocol for the Guangzhou prospective cohort study on hand-foot-and-mouth disease. BMJ Open. Feb 24, 2019;9(2):e027224. [CrossRef] [Medline]
  23. Chen Y, Wang F, Jin Y. Coxsackievirus A6-induced central nervous system injury and limb paralysis. Pediatr Infect Dis J. May 7, 2024. [CrossRef] [Medline]
  24. Chen Y, Chen S, Duan G, et al. A recombinant CVA6 infection aggravates left ventricular dysfunction and myocardial injury in a child with cardiomyopathy. Pediatr Infect Dis J. May 1, 2024;43(5):e191-e193. [CrossRef] [Medline]
  25. Kimmis BD, Downing C, Tyring S. Hand-foot-and-mouth disease caused by coxsackievirus A6 on the rise. Cutis. Nov 2018;102(5):353-356. [Medline]
  26. Zhao TS, Du J, Sun DP, et al. A review and meta-analysis of the epidemiology and clinical presentation of coxsackievirus A6 causing hand-foot-mouth disease in China and global implications. Rev Med Virol. Mar 2020;30(2):e2087. [CrossRef] [Medline]
  27. He YQ, Chen L, Xu WB, et al. Emergence, circulation, and spatiotemporal phylogenetic analysis of coxsackievirus A6- and coxsackievirus A10-associated hand, foot, and mouth disease infections from 2008 to 2012 in Shenzhen, China. J Clin Microbiol. Nov 2013;51(11):3560-3566. [CrossRef] [Medline]
  28. Tan X, Huang X, Zhu S, et al. The persistent circulation of enterovirus 71 in People’s Republic of China: causing emerging nationwide epidemics since 2008. PLoS One. Sep 28, 2011;6(9):e25662. [CrossRef] [Medline]
  29. Kumar S, Stecher G, Tamura K. MEGA7: Molecular Evolutionary Genetics Analysis version 7.0 for bigger datasets. Mol Biol Evol. Jul 2016;33(7):1870-1874. [CrossRef] [Medline]
  30. IQ-TREE version 1.6.12. IQ-TREE. URL: http://www.iqtree.org/release/v1.6.12 [Accessed 2024-07-18]
  31. Trifinopoulos J, Nguyen LT, von Haeseler A, Minh BQ. W-IQ-TREE: a fast online phylogenetic tool for maximum likelihood analysis. Nucleic Acids Res. Jul 8, 2016;44(W1):W232-W235. [CrossRef] [Medline]
  32. Rambaut A, Lam TT, Max Carvalho L, Pybus OG. Exploring the temporal structure of heterochronous sequences using TempEst (formerly Path-O-Gen). Virus Evol. Apr 9, 2016;2(1):vew007. [CrossRef] [Medline]
  33. Suchard MA, Lemey P, Baele G, Ayres DL, Drummond AJ, Rambaut A. Bayesian phylogenetic and phylodynamic data integration using BEAST 1.10. Virus Evol. Jun 8, 2018;4(1):vey016. [CrossRef] [Medline]
  34. Yu L, Guo Q, Wei H, et al. Molecular epidemiology and evolution of coxsackievirus A14. Viruses. Nov 26, 2023;15(12):2323. [CrossRef] [Medline]
  35. Rambaut A. FigTree. Molecular evolution, phylogenetics and epidemiology. URL: http://tree.bio.ed.ac.uk/software/figtree [Accessed 2024-07-18]
  36. BioEdit 7.2. Software Informer. URL: https://bioedit.software.informer.com/7.2 [Accessed 2024-07-18]
  37. Zhang Y, Wang J, Guo W, et al. Emergence and transmission pathways of rapidly evolving evolutionary branch C4a strains of human enterovirus 71 in the Central Plain of China. PLoS One. Nov 18, 2011;6(11):e27895. [CrossRef] [Medline]
  38. Flett K, Youngster I, Huang J, et al. Hand, foot, and mouth disease caused by coxsackievirus A6. Emerg Infect Dis. Oct 2012;18(10):1702-1704. [CrossRef] [Medline]
  39. Gaunt E, Harvala H, Österback R, et al. Genetic characterization of human coxsackievirus A6 variants associated with atypical hand, foot and mouth disease: a potential role of recombination in emergence and pathogenicity. J Gen Virol. May 2015;96(Pt 5):1067-1079. [CrossRef] [Medline]
  40. Gnädig NF, Beaucourt S, Campagnola G, et al. Coxsackievirus B3 mutator strains are attenuated in vivo. Proc Natl Acad Sci U S A. Aug 21, 2012;109(34):E2294-E2303. [CrossRef] [Medline]
  41. Lauring AS, Andino R. Quasispecies theory and the behavior of RNA viruses. PLoS Pathog. Jul 22, 2010;6(7):e1001005. [CrossRef] [Medline]
  42. Pfeiffer JK, Kirkegaard K. Increased fidelity reduces poliovirus fitness and virulence under selective pressure in mice. PLoS Pathog. Oct 2005;1(2):e11. [CrossRef] [Medline]
  43. Upala P, Apidechkul T, Suttana W, Kullawong N, Tamornpark R, Inta C. Molecular epidemiology and clinical features of hand, foot and mouth disease in northern Thailand in 2016: a prospective cohort study. BMC Infect Dis. Dec 6, 2018;18(1):630. [CrossRef] [Medline]
  44. Li T, Yang Z, DI B, Wang M. Hand-foot-and-mouth disease and weather factors in Guangzhou, southern China. Epidemiol Infect. Aug 2014;142(8):1741-1750. [CrossRef] [Medline]
  45. Puenpa J, Vongpunsawad S, Österback R, et al. Molecular epidemiology and the evolution of human coxsackievirus A6. J Gen Virol. Dec 2016;97(12):3225-3231. [CrossRef] [Medline]
  46. Tomba Ngangas S, Bisseux M, Jugie G, et al. Coxsackievirus A6 recombinant subclades D3/A and D3/H were predominant in hand-foot-and-mouth disease outbreaks in the paediatric population, France, 2010-2018. Viruses. May 17, 2022;14(5):1078. [CrossRef] [Medline]
  47. Carmona RCC, Machado BC, Reis FC, et al. Hand, foot, and mouth disease outbreak by coxsackievirus A6 during COVID-19 pandemic in 2021, São Paulo, Brazil. J Clin Virol. Sep 2022;154:105245. [CrossRef] [Medline]
  48. Zhang M, Chen X, Wang W, Li Q, Xie Z. Genetic characteristics of coxsackievirus A6 from children with hand, foot and mouth disease in Beijing, China, 2017-2019. Infect Genet Evol. Dec 2022;106:105378. [CrossRef] [Medline]
  49. Saxena VK, Pawar SD, Qureshi T, et al. Isolation and molecular characterization of coxsackievirus A6 and coxsackievirus A16 from a case of recurrent hand, foot and mouth disease in Mumbai, Maharashtra, India, 2018. Virusdisease. Mar 2020;31(1):56-60. [CrossRef] [Medline]
  50. Mizuta K, Tanaka S, Komabayashi K, et al. Phylogenetic and antigenic analyses of coxsackievirus A6 isolates in Yamagata, Japan between 2001 and 2017. Vaccine. Feb 14, 2019;37(8):1109-1117. [CrossRef] [Medline]
  51. Hoa-Tran TN, Dao ATH, Nguyen AT, et al. Coxsackieviruses A6 and A16 associated with hand, foot, and mouth disease in Vietnam, 2008-2017: essential information for rational vaccine design. Vaccine. Dec 14, 2020;38(52):8273-8285. [CrossRef] [Medline]
  52. Noisumdaeng P, Puthavathana P. Molecular evolutionary dynamics of enterovirus A71, coxsackievirus A16 and coxsackievirus A6 causing hand, foot and mouth disease in Thailand, 2000-2022. Sci Rep. Oct 13, 2023;13(1):17359. [CrossRef] [Medline]
  53. Muslin C, Mac Kain A, Bessaud M, Blondel B, Delpeyroux F. Recombination in enteroviruses, a multi-step modular evolutionary process. Viruses. Sep 14, 2019;11(9):859. [CrossRef] [Medline]
  54. Tan X, Li L, Zhang B, et al. Molecular epidemiology of coxsackievirus A6 associated with outbreaks of hand, foot, and mouth disease in Tianjin, China, in 2013. Arch Virol. Apr 2015;160(4):1097-1104. [CrossRef] [Medline]
  55. Xu L, Zheng Q, Li S, et al. Atomic structures of coxsackievirus A6 and its complex with a neutralizing antibody. Nat Commun. Sep 11, 2017;8(1):505. [CrossRef] [Medline]
  56. Zhou Y, Van Tan L, Luo K, et al. Genetic variation of multiple serotypes of enteroviruses associated with hand, foot and mouth disease in Southern China. Virol Sin. Feb 2021;36(1):61-74. [CrossRef] [Medline]


CV: coxsackievirus
EV: enterovirus
HFMD: hand, foot, and mouth disease
HPD: highest posterior density
MCC: maximum clade credibility
MCMC: Markov chain Monte Carlo
ML: maximum likelihood
tMRCA: time of the most recent common ancestor


Edited by Amaryllis Mavragani; submitted 17.04.24; peer-reviewed by Hehe Zhao, Sanjaykumar Tikute, Tianjiao Ji; final revised version received 19.06.24; accepted 20.06.24; published 31.07.24.

Copyright

© Yu Chen, Shouhang Chen, Yuanfang Shen, Zhi Li, Xiaolong Li, Yaodong Zhang, Xiaolong Zhang, Fang Wang, Yuefei Jin. Originally published in JMIR Public Health and Surveillance (https://publichealth.jmir.org), 31.7.2024.

This is an open-access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work, first published in JMIR Public Health and Surveillance, is properly cited. The complete bibliographic information, a link to the original publication on https://publichealth.jmir.org, as well as this copyright and license information must be included.