Molecular analysis of Coxsackievirus A24 variant isolates from three outbreaks of acute hemorrhagic conjunctivitis in 1988, 1994 and 2007 in Beijing, China.

Coxsackievirus A24 variant (CVA24v) is a major pathogen that causes continued outbreaks and pandemics of acute hemorrhagic conjunctivitis (AHC). In China, the first confirmed outbreak of CVA24v-related AHC occurred in Beijing in 1988, followed by another two significant outbreaks respectively in 1994 and 2007, which coincides with the three-stage dynamic distribution of AHC in the world after 1970s. To illustrate the genetic characteristics of CVA24v in different periods, a total of 23 strains were isolated from those three outbreaks and the whole genome of those isolations were sequenced and analyzed. Compared with the prototype strain, the 23 strains shared four nucleotide deletions in the 5' UTR except the 0744 strain isolated in 2007. And at the 98th site, one nucleotide insertion was found in all the strains collected from 2007. From 1994 to 2007, amino acid polarity in the VP1 region at the 25th and the 32nd site were changed. Both the 3C and VP1 phylogenetic tree indicated that isolates from 1988 and 1994 belonged to Genotype III (GIII), and 2007 strains to Genotype IV (GIV). According to the Bayesian analysis based on complete genome sequence, the most recent common ancestors for the isolates in 1988, 1994 and 2007 were respectively estimated around October 1987, February 1993 and December 2004. The evolutionary rate of the CVA24v was estimated to be 7.45 ​× ​10-3 substitutions/site/year. Our study indicated that the early epidemic of CVA24v in Chinese mainland was the GIII. Point mutations and amino acid changes in different genotypes of CVA24v may generate intensity differences of the AHC outbreak. CVA24v has been evolving constantly with a relatively rapid rate.

Introduction

Acute hemorrhagic conjunctivitis (AHC) is a highly contagious acute viral eye disease. Repeated outbreaks of AHC have been frequently reported on a global scale. Coxsackievirus A24 variant (CVA24v) has been identified as the major causative agent of AHC outbreaks (Burr et al., 2017; Fonseca et al., 2012; Miyamura et al., 1988). As an antigenic variant of CVA24, CVA24v was first isolated from the AHC outbreak in Singapore in 1970 (Lim and Yin-Murphy, 1971). In the subsequent decade, the AHC epidemic caused by CVA24v had been mainly confined to Southeast Asia and India (Christopher et al., 1982; Miyamura et al., 1988). However, this confinement was broken after the outbreak of AHC in Singapore in 1985, after which the epidemic caused by CVA24v quickly spread to Asia, America and Africa (Sawyer et al., 1989; Tavares et al., 2006; Lévêque et al., 2007). The dynamic process of the AHC caused by CVA24v in the world after 1970s has three characteristically distinct periods: The rapid expansion period with a large number of cases reported in the 1980s; The relatively silent period with moderate outbreaks occurred in 1990s; The re-emerging period with extremely rapid spread after 2000 (Yen et al., 2016).

In China, the first outbreak of AHC was reported in 1971, but the pathogen of that outbreak was unknown. Except for the remote areas, that outbreak almost spread to all provinces and cities across the country, including Hong Kong and Taiwan. Another outbreak of AHC occurred in Hong Kong in 1975 (Chang et al., 1977). Then in 1986, AHC outbreaks spread to Beijing and Guangzhou (Mu et al., 1992). In 1988, an outbreak of AHC in Guangzhou was reported, with the estimated actual number of cases reaching more than 490,000 within 70 days. AHC was also reported in Qingdao in 1994 and 1997 (Meng et al., 1999). In 2002, the local epidemic of AHC occurred again in coastal areas such as Zhejiang and Fujian (Li et al., 2013). In Yunnan and Jiangsu, there were also reported AHC in 2007 and 2010, respectively.

Although AHC occurs frequently, those outbreaks were caused by EV70 or certain serotypes of adenoviruses. It was not until 1988 in Beijing AHC outbreak, CVA24v was firstly identified in China (Mu et al., 1990). During that large outbreak, nearly 400,000 cases were reported, with an incidence of more than 10% (Shi and He, 1999). Therefore, in 1989, a monitoring system of AHC was established throughout Beijing. According to the system, during the period from 1989 to 2006, except for a moderate outbreak of AHC caused by CVA24v with 18101 cases reported (incidence of 1.71%) in 1994, the intensity of AHC epidemic in Beijing had been relatively low (Shi and He, 1999; Wang et al., 2014). However, since a nationwide AHC epidemic caused by CVA24v occurred in China in 2007, the incidence has risen sharply, which is demonstrated by a total of 134 outbreaks covering 14 provinces including Beijing and a total of 74,263 case reports, representing an increase by 540.58% compared with 11,593 cases in 2006 (Wang et al., 2014). In that massive outbreak of AHC in 2007, the reported incidence rate of Beijing was 3.6%, ranking the 10th in China. And all the 18 districts and counties in Beijing had reported cases (Ding et al., 2009). The three outbreaks of AHC in 1988, 1994 and 2007 in Beijing happened to coincide with the three characteristically distinct stages of the AHC outbreak in the world, which respectively are the rapid expansion period in the 1980s, the relatively silent period in 1990s and the re-emerging period after 2000. This would be a golden opportunity for us to explore the evolution of CVA24v in different periods.

As the first phylogenetic analysis of CVA24v genomes was based on the non-structural 3C region (Ishiko et al., 1992), therefore, CVA24v has been genetically classified into four genotypes based on the 3C region. Genotype I (GI) was detected in Singapore and Hong Kong from 1970 to 1971, Genotype II (GII) in Singapore and Thailand in 1975 and Genotype III (GIII) in France, Asian and African countries from 1985 to 1994, whereas Genotype IV (GIV) was detected after 2000 (Nidaira et al., 2014). However, since VP1 region is regarded as the target gene for enterovirus genotyping, phylogenetic analyses based on the VP1 region of CVA24v has surged, especially post-2000. Subsequently, CVA24v has also been classified into four genotypes based on the entire VP1 region. GI was detected in Singapore in 1970, GII in Brazil and Jamaica in 1987 and GIII in the Dominican Republic in 1993 and in the USA in 1998, whereas GIV were detected after 2002 (De et al., 2012; Nidaira et al., 2014). Controversy about whether the 3C or VP1 region of CV-A24v was more informative has lasted for a long time because of the limited VP1 sequences (Fonseca et al., 2020).

In our study, we added 23 whole genome sequences of CVA24v, which enriched the gene pool of CVA24v sequences to contribute to more accurate phylogenetic and genotyping analyses. On the other hand, 14 strains were isolated from two of the three AHC outbreaks caused by CVA24v occurred before 2000, which enables us to identify the genotypes from early epidemics in Chinese mainland. To sum up, by this study, we can confirm the genotypes in Chinese mainland before 2000 as well as explore the evolutionary characteristics of CVA24v in different periods.

Materials and methods

Specimen collection and virus isolation

During three Beijing AHC outbreaks in 1988, 1994 and 2007, conjunctival swab specimen was collected within 3 days from outpatient clinically diagnosed with AHC. Identification of CVA24v was completed by Beijing Center for Disease Control and Prevention (CDC). Among all the samples that have been identified as CVA24v, parts of them were sent to the National Polio Laboratory and stored under −70 ​°C. We took them for further virus amplification. Two hundred microliter of the sample solution was inoculated onto HEp-2 ​cells (provided by the WHO Global Poliovirus Specialized Laboratory in the US and purchased from the American Type Culture Collection, Manassas, VA, USA), then observed for 7 days for cytopathic effects (CPE). The culture temperature was 36 ​°C with 5% CO2. After complete CPE was observed, infected cell cultures would be harvested.

Sequencing and genetic analysis

The viral nucleic acid was extracted with QIAamp Viral Mini Kit (QIAGEN, Valencia, CA, USA) according to the instruction for use. The full-length genomes of the 23 CVA24v were amplified by specific primers (De et al., 2012). The products were purified with QIAquick Gel Extraction Kit (QIAGEN, Valencia, CA, USA), then labelled with BigDye Terminator v3.0 Cycle Sequencing Kit (Applied Biosystems Foster City, USA). After purification by dextran gel G-50 (Pharmacia, Sweden), the labelled product was sequenced on an ABI3100 sequencing analyzer (Applied Biosystems, Hitachi, Japan).

The sequence analysis, correction and collation were carried out by Sequencher 5.4.5 software (Gene Code, Ann Arbor, Michigan, USA). The nucleotide sequence alignment and homology analysis were conducted by MEGA7.0 software (Sudhir Kumar, Arizona State University, Arizona, USA). The phylogenetic trees of VP1 and 3C for genotyping were constructed by the maximum likelihood (ML) method in MEGA7.0 software with the reliability evaluated by 1000 bootstrap (Kumar et al., 2016). The reference strains of VP1 and 3C were retrieved and downloaded from GenBank. The ML tree of the VP1 region was constructed based on 23 strains from Beijing and 39 worldwide strains from GenBank (Supplementary Table S2). And the ML tree of the 3C region was based on 23 strains from this study and 59 from GenBank (Supplementary Table S3).

Phylogenetic analysis and estimation of evolutionary rate

To estimate the molecular evolutionary rate and relationship among CVA24v isolates, phylogenetic analysis was conducted based on the full-length genome sequences of the 23 strains of Beijing and the 42 whole genome sequences collected from GenBank (Supplementary Table S4). Evolutionary rate was estimated using the Bayesian Markov chain Monte Carlo (MCMC) method in BEAST software (version1.7.5) (Drummond et al., 2012). The dataset was analyzed under a lognormal relaxed uncorrelated clock using the General Time Reversible (GTR) ​+ ​Gamma (G) ​+ ​Invariable sites (I) model supported by the jModelTest software (version 2) (Darriba et al., 2012). The MCMC chain was run for 200,000,000 steps and sampled every 20,000 steps. Model evaluation was performed by Tracer v1.6 software (Rambaut et al., 2018). The maximum clade credibility (MCC) tree was viewed and edited by the Figtree.v1.4.4 software (Drummond et al., 2012).

Results

Viral isolation and genetic analysis

Twenty-three samples that have been identified as CVA24v by Beijing CDC were selected for virus amplification in HEp-2 ​cell line in the National Polio Laboratory. CPE of the 23 samples were all observed in seven days. Then seven, seven and nine infected cell cultures from 1988, 1994 and 2007 Beijing AHC outbreak were harvested respectively and then the strains were isolated for further research. The complete genome sequences of the 23 Beijing CVA24v strains were determined. They are 7456–7457 ​nt in length with a 5′ UTR of about 750 ​nt and a 3′ UTR of about 69 ​nt. A long open reading frame of 6642 ​nt, corresponding to 2214 aa, encoded the polyprotein precursor. The background information of the isolated strains was listed in Supplementary Table S1. The nucleotide sequence homologies of the 23 strains were 88.9%～99.9% corresponding to 96.3%～100% of amino acid homology. The nucleotide and amino acid sequence homologies compared with the prototype strain are shown in Table 1.

Table 1

Comparison of nucleotide and amino acid similarities between the prototype strain and the Beijing strains.

Year	Number of strain	% identity with D90457.1/Singapore/1970
		Nucleotide	Amino acid
1988	7	87.5–87.6	96.8–96.9
1994	7	87.1–87.2	96.5–96.8
2007	9	85.7–86.1	96.3–96.7

Point mutations and amino acid changes

The 23 strains shared four nucleotide deletions from the 116th site to the 119th site in the 5′ UTR except the strain of 0744 isolated in 2007 compared with the prototype strain (EH24/70, GenBank number: D90457). Furthermore, at the 98th site, one nucleotide insertion was found in all the strains of 2007. No nucleotide insertion and deletion were detected in the 3′ UTR. All the 23 strains shared 43 amino acid mutations in polyprotein precursors compared with the prototype strain. The amino acid changes and frequencies in different regions of the genome are shown in Fig. 1. Mutated amino acid residues were distributed mainly in the VP1 region, in which 11 of the 43 mutations shared by all the Beijing isolates compared to the prototype sequence. The amino acid changes in the VP1 region compared with the prototype are shown in Table 2. The total length of the VP1 region and 3C region was 915 and 549 nucleotides respectively, encoding 305 and 183 amino acids. Two common amino acid changes were observed in the 3C region compared with the prototype as shown in Table 3.

Fig. 1

Summary of amino acid changes and frequencies in different regions of the genome of CVA24v among the prototype strain (D90457.1/Singapore/1970) and strains from 1988 (n ​= ​7), 1994 (n ​= ​7) and 2007 (n ​= ​9).

Table 2

Amino acid mutations in the VP1 region.

Strains	Amino acid site in the VP1 region
	11	16	25	32	51	56	89	100	103	146	147	151	188	196	250	256	297	301	305
D90457.1-VP1/Singapore/1970	S	L	L	S	V	I	M	E	K	T	S	Y	T	I	F	I	I	N	F
8809-VP1_1988	T	.	S	P	A	V	I	.	.	A	.	H	.	M	Y	T	T	.	L
8811-VP1_1988	T	.	S	P	A	V	I	.	.	A	.	H	.	M	Y	T	T	.	L
8812-VP1_1988	T	.	S	P	A	V	I	.	.	A	N	H	.	M	Y	T	T	.	L
8814-VP1_1988	T	.	S	P	A	V	I	.	.	A	.	H	.	M	Y	T	T	.	L
8816-VP1_1988	T	.	S	P	A	V	I	.	.	A	.	H	.	M	Y	T	T	.	L
8818-VP1_1988	T	.	S	P	A	V	I	.	.	A	.	H	.	M	Y	T	T	.	L
8819-VP1_1988	T	.	S	P	A	V	I	.	.	A	.	H	.	M	Y	T	T	.	L
9406-VP1_1994	T	F	S	P	A	V	I	.	.	A	.	H	.	M	Y	T	T	.	L
9409-VP1_1994	T	.	S	P	A	V	I	.	.	A	.	H	.	M	Y	T	T	.	L
9410-VP1_1994	T	F	S	P	A	V	I	.	.	A	.	H	.	M	Y	T	T	.	L
9414-VP1_1994	T	F	P	P	A	V	I	.	.	A	.	H	.	M	Y	T	T	.	L
9415-VP1_1994	T	F	S	P	A	V	I	.	.	A	.	H	S	M	Y	T	T	.	L
9417-VP1_1994	T	F	S	P	A	V	I	.	.	A	.	H	.	M	Y	T	T	.	L
9418-VP1_1994	T	F	S	P	A	V	I	.	.	A	.	H	.	M	Y	T	T	.	L
0740-VP1_2007	T	.	H	L	A	V	I	D	R	A	.	H	.	M	Y	T	T	.	L
0769-VP1_2007	T	.	H	L	A	V	I	D	R	A	.	H	.	M	Y	T	T	.	L
0771-VP1_2007	T	.	H	L	A	V	I	D	R	A	.	H	.	M	Y	T	T	.	L
0775-VP1_2007	T	.	H	L	A	V	I	D	R	A	.	H	.	M	Y	T	T	.	L
0778-VP1_2007	T	.	H	L	A	V	I	D	R	A	.	H	.	M	Y	T	T	.	L
0779-VP1_2007	T	.	H	L	A	V	I	D	R	A	.	H	.	M	Y	T	T	.	L
0783-VP1_2007	T	.	H	L	A	V	I	D	R	A	.	H	.	M	Y	T	T	.	L
0784-VP1_2007	T	.	H	L	A	V	I	D	R	A	.	H	.	M	Y	T	T	.	L
0744-VP1_2007	T	.	H	L	A	V	I	D	R	A	.	H	.	M	Y	T	T	D	L

Table 3

Amino acid mutations in the 3C region.

Strains	Amino acid site in the 3C region
	15	31	42	52	68	92	114	139	151	160
D90457.1–3C/Singapore/1970	I	H	A	K	T	T	V	N	I	M
8809-3C_1988	V	.	.	.	.	A	.	H	.	.
8811-3C_1988	V	.	.	.	.	A	.	H	.	.
8812-3C_1988	V	.	.	.	.	A	.	H	.	.
8814-3C_1988	V	.	.	.	.	A	.	H	.	.
8816-3C_1988	V	.	.	.	.	A	.	H	.	.
8818-3C_1988	V	.	.	.	.	A	.	H	.	.
8819-3C_1988	V	.	.	.	.	A	.	H	.	.
9406-3C_1994	V	.	.	.	I	A	.	H	.	.
9409-3C_1994	V	.	.	.	I	A	.	H	.	.
9410-3C_1994	V	.	.	.	I	A	.	H	.	.
9414-3C_1994	V	Y	.	.	I	A	.	H	.	.
9415-3C_1994	V	Y	.	.	I	A	.	H	.	.
9417-3C_1994	V	.	.	.	I	A	.	H	.	.
9418-3C_1994	V	.	.	.	I	A	.	H	.	.
0740-3C_2007	V	Y	.	.	V	.	I	H	.	I
0769-3C_2007	V	Y	.	.	V	.	I	H	.	I
0771-3C_2007	V	Y	.	.	V	.	I	H	.	I
0775-3C_2007	V	Y	V	.	V	.	I	H	.	I
0778-3C_2007	V	Y	.	.	V	.	I	H	.	I
0779-3C_2007	V	Y	.	.	V	.	I	H	.	I
0783-3C_2007	V	Y	.	.	V	.	I	H	.	I
0784-3C_2007	V	Y	.	.	V	.	I	H	.	I
0744-3C_2007	V	Y	V	R	I	A	I	H	V	I

To further investigate the genetic relationship of the CVA24v causing three outbreaks in Beijing, genotyping was performed based on the 3C region and the VP1 region. We retrieved and downloaded all CVA24v sequences from GenBank and removed part of highly similar sequences isolated from the same outbreak. Then we obtained 59, 39 and 42 sequences respectively of 3C, VP1 and whole genome (Supplementary Table S2, S3 and S4). The ML tree of the 3C region was constructed based on 23 strains from Beijing and 59 worldwide strains from GenBank (Fig. 2). And the ML tree of the VP1 region was based on 23 strains from this study and 39 from GenBank (Fig. 3). Both of the phylogenetic trees showed that the strains from 1988 and 1994 belonged to GIII, and the strains from 2007 belonged to GIV. To assess time-scale evolution, phylogenetic analysis was performed using the Bayesian MCMC method. The MCC tree was constructed based on the whole genome sequences of the 23 strains from Beijing and 42 worldwide strains from GenBank (Fig. 4).

Fig. 2

Phylogenetic analysis based on entire 3C region. 3C phylogenetic tree was constructed by the maximum likelihood method in the MEGA 7.0 program with 59 worldwide 3C strains from GenBank and 23 Beijing strains in this study. The bootstrap values in 1000 replicates for major lineages are shown as percentages. Only bootstrap values＞80% was shown at the nodes. Where and when each strain was isolated are reported next to its GenBank accession number. ▲, ●, ♦ represent the 1988, 1994 and 2007 isolates, respectively.

Fig. 3

Phylogenetic analysis based on entire VP1 region. VP1 phylogenetic tree was constructed by the maximum likelihood method in the MEGA 7.0 program with 39 worldwide VP1 strains from GenBank and 23 Beijing strains in this study. The bootstrap values in 1000 replicates for major lineages are shown as percentages. Only bootstrap values＞80% was shown at the nodes. Where and when each strain was isolated are reported next to its GenBank accession number. ▲, ●, ♦ represent the 1988, 1994 and 2007 isolates, respectively.

Fig. 4

Phylogenetic analysis based on whole genome sequences. The MCC phylogenetic tree was constructed based on the whole genome sequences of 23 Beijing (in red) and 42 worldwide CV-A24v strains. For each branch, the color indicates the evolutionary rate. 95% HPDs of TMRCAs are indicated by the bars at nodes. Where and when each strain was isolated are reported next to its GenBank accession number.

The time for the most recent common ancestor (tMRCA) (95% highest probability density, 95% HPD) of the Beijing strains was estimated based on divergence time analysis with whole genome sequences. The tMRCAs of the isolates in 1988, 1994 and 2007 were respectively estimated to be from Oct. 1987 (Sep. 1987–Dec. 1987), Feb. 1993 (Oct. 1992–Jun. 1993) and Dec. 2004 (Feb. 2004–Sep. 2005). Furthermore, the rate of molecular evolution of the whole genome sequence based on the 23 Beijing isolates and 42 worldwide strains was estimated 7.45 (6.26–8.55) ​× ​10−3 substitutions/site/year. The rate of molecular evolution of the 42 worldwide whole genome sequences was 6.28 (4.93–7.45) ​× ​10−3 substitutions/site/year while the evolutionary rate of the 23 Beijing isolates in our study was 8.12 (7.49–8.70) ​× ​10−3 substitutions/site/year.

Discussion

Since 1970s, reports on pandemic of AHC, a highly contagious eye disease with no vaccines or antiviral drugs currently available, were mainly related to CVA24v. During the preparation for the 11th Asian Games in 1990, Beijing carried out a number of targeted thematic studies, including the epidemic characteristics and the control and prevention measures of AHC. And the research results were applied to the health and disease prevention for the Asian Games. Since then, AHC has been listed as one of the Class C infectious diseases in China (Li et al., 1993). All AHC cases diagnosed by physicians should be registered (Yan et al., 2010). The Beijing monitoring system of AHC captured the early outbreaks with sufficient information and dynamic data, which allowed us to conduct study on AHC epidemics caused by CVA24v in different times.

The epidemiology of the three AHC outbreaks in Beijing is consistent with that reported worldwide. The research based on this opportunity shows that these three different stages correspond to the prevalence of different genotypes of CVA24v at the level of molecular evolution. Due to high similarity among the VP1 sequences of CVA24v, the criteria used for enterovirus genotyping “a cluster of genetically related viruses with less than 15% nt divergence” is inadequate (Yen et al., 2016). Debate about whether the 3C or VP1 region of CVA24v was more informative has lasted for a long time (De et al., 2012; Nidaira et al., 2014). The genotype II is their point of contention. Some of them held the view that the GII appeared around 1975, while others classified the two isolates from 1987 as GII (Chu et al., 2009; De et al., 2012; Nidaira et al., 2014). In our study, we confirmed that the results of phylogenetic analysis based on both 3C and VP1 region are consistent with the complete genome sequence. The worldwide CVA24v can be divided into four genotypes. The genotype of CVA24v in the rapid expansion period in the 1980s and the relatively silent period in the 1990s was all GIII, while the re-emerging period after 2000 was GIV. The lack of VP1 sequences of early times, especially around 1975, lead to differences of the genotyping.

Although 23 strains in this study were isolated from the same place, the intensity of the three AHC outbreaks varied. The 2007 outbreak of AHC is clearly more widespread. In terms of genetic sequence, the 5' UTR and 3' UTR of the enterovirus genome contribute to viral RNA replication. The 5' UTR contains an oriL (nt 1–90) and an IRES (nt 91–745) that guide initiation of translation in a way that independent of the cap structure (Lin et al., 2009). In this study, we found a nucleotide insertion at the 98th site of the 5' UTR of all 2007 isolates, whereas this insertion was not found in the genomes of other isolates in 1988 and 1994. Since the inserted nucleotide was in the IRES of CVA24v, we hypothesize that this insertion may have an effect on the initiation of translation of CVA24v. The VP1, VP2 and VP3 region of CVA24v are located on the surface of the viral capsid and exposed to immune pressure. The major neutralization sites are located in VP1, VP2 and VP3 proteins and VP1 contains most neutralization epitopes. In this study, we found changes of the amino acid polarity in VP1 region at the 25th site (L→S→H) and the 32nd site (S→P→L) among the prototype sequence, the sequences of 1988 and 1994 strains (GIII) and the sequences of 2007 strains (GIV). The change of amino acids from hydrophobic to hydrophilic may affect the secondary structure and epitopes of capsid proteins and eventually lead to escape mutations that promote the outbreak of AHC in a relatively short period of time. Results of the study on CVA24v by De et al. also supported our view of that point mutations and amino acid changes in different genotypes may affected the translation of CVA24v thereby generated the intensity differences of the AHC outbreak (De et al., 2012).

We used the MCMC method to assess the molecular evolutionary rate and the most recent common ancestors (95% HPD) of the Beijing strains based on complete genome sequences. According to our analysis, the tMRCA of the isolates in 1988, 1994 and 2007 were respectively estimated to be from Oct. 1987, Feb. 1993 and Dec. 2004. The evolutionary rate of the CVA24v strains was estimated to be 7.45 (6.26–8.55) ​× ​10−3 substitutions/site/year. The previous reported evolutionary rate based on the VP1 and 3C region was respectively 6.37 (5.17–7.60) ​× ​10−3 substitutions/site/year and 6.15 (4.38–8.07) ​× ​10−3 substitutions/site/year (Nidaira et al., 2014). Actually the 42 worldwide CVA24v strains collected from GenBank evolved at a similar evolutionary rate of 6.28 (4.93–7.45) ​× ​10−3 substitutions/site/year with the report when the Beijing strains were not considered. The evolutionary rate of 23 Beijing CVA24v isolates in our study was 8.12 (7.49–8.70) ​× ​10−3 substitutions/site/year. With such a so high evolutionary rate, it is difficult for people to have sustained immunity against CVA24v. The loss of herd immunity against the virus may lead to widespread transmission. Therefore, the establishment of a new genotype often leads to an outbreak.

In conclusion, our study indicated that phylogenetic analysis based on 3C and VP1 region are consistent with the complete genome sequence. The strains from early epidemic of CVA24v in Chinese mainland are in GIII. CVA24v has been evolving constantly with a relatively rapid rate. Point mutations and amino acid changes in different genotypes of CVA24v may generate intensity differences of the AHC outbreak. Further studies on the genetic characterization and molecular evolution of CVA24v are necessary to prevent and control the outbreaks of AHC in China.

Data availability

All data relevant to the study are included in the article or uploaded as supplementary information. The VP1, 3C and complete genome sequences of 23 Beijing strains from this study are available in GenBank with the following accession numbers: MT741560–MT741582, MT741537–MT741559 and MZ171074–MZ171096.

Ethics statement

The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Ethical Committee of National Institute for Viral Diseases Control and Prevention, Center for Disease Control and Prevention, China. Informed consent was waived because the research is a retrospective study conducted with viral isolates obtained in routine diagnostics without accessing patient medical data.

Author contributions

Junhan Li: conceptualization, data curation, writing—original draft. Fang Huang: resources, supervision. Yong Zhang: methodology, funding acquisition, project administration. Tianjiao Ji: validation. Shuangli Zhu: resources. Dongyan Wang: resources. Zhenzhi Han: software. Jinbo Xiao: software. Fenfen Si: validation. Wenbo Xu: funding acquisition, supervision. Dongmei Yan: conceptualization, project administration, writing—review & editing.

Conflicts of interest

The authors declare that there are no conflicts of interest.