Advertisement

Evolutionary phylogeography reveals novel genotypes of coxsackievirus A24 variant and updates the spatiotemporal dynamics in the population with acute hemorrhagic conjunctivitis

  • Peng Chen
    Affiliations
    Department of Healthcare-associated Infection Management, School and Hospital of Stomatology, Cheeloo College of Medicine, Shandong University & Shandong Key Laboratory of Oral Tissue Regeneration & Shandong Engineering Laboratory for Dental Materials and Oral Tissue Regeneration & Shandong Provincial Clinical Research Center for Oral Diseases, Jinan, China
    Search for articles by this author
  • Xiao-Juan Lin
    Affiliations
    Shandong Provincial Key Laboratory of Infectious Disease Control and Prevention; Shandong Center for Disease Control and Prevention, Jinan, China
    Search for articles by this author
  • Feng Ji
    Affiliations
    Shandong Provincial Key Laboratory of Infectious Disease Control and Prevention; Shandong Center for Disease Control and Prevention, Jinan, China
    Search for articles by this author
  • Yan Li
    Affiliations
    Shandong Provincial Key Laboratory of Infectious Disease Control and Prevention; Shandong Center for Disease Control and Prevention, Jinan, China
    Search for articles by this author
  • Su-Ting Wang
    Affiliations
    Shandong Provincial Key Laboratory of Infectious Disease Control and Prevention; Shandong Center for Disease Control and Prevention, Jinan, China
    Search for articles by this author
  • Yao Liu
    Affiliations
    Shandong Provincial Key Laboratory of Infectious Disease Control and Prevention; Shandong Center for Disease Control and Prevention, Jinan, China
    Search for articles by this author
  • Ze-Xin Tao
    Correspondence
    Ze-Xin Tao, Shandong Center for Disease Control and Prevention, No. 16992 Jingshi Road, Jinan 250014, People's Republic of China; Tel: +86 531 82679693; Fax: +86 531 82679693
    Affiliations
    Shandong Provincial Key Laboratory of Infectious Disease Control and Prevention; Shandong Center for Disease Control and Prevention, Jinan, China
    Search for articles by this author
  • Ai-Qiang Xu
    Correspondence
    Corresponding authors: Ai-Qiang Xu, Shandong Center for Disease Control and Prevention, No. 16992 Jingshi Road, Jinan 250014, People's Republic of China, Tel: +86 531 82679606; Fax: +86 531 82620031
    Affiliations
    Shandong Provincial Key Laboratory of Infectious Disease Control and Prevention; Shandong Center for Disease Control and Prevention, Jinan, China
    Search for articles by this author
Open AccessPublished:October 11, 2022DOI:https://doi.org/10.1016/j.ijid.2022.10.007

      Highlights

      • Eight coxsackievirus A24 variant (CVA24v) genotypes were divided, and two novel genotypes were reported first.
      • Five genotypes of CVA24v circulated previously in mainland China.
      • Special ocular tropism might impact the variation of CVA24v.
      • Residue 146T could contribute to the CVA24v outbreaks and its pathogenicity.
      • CVA24v possesses a potential neuroinvasion capability besides causing acute hemorrhagic conjunctivitis.

      Abstract

      Objectives

      The coxsackievirus A24 variant (CVA24v) has raised a remarkable concern because of its main etiological role in acute hemorrhagic conjunctivitis.

      Methods

      We conducted a retrospective study to summarize CVA24v isolated from acute hemorrhagic conjunctivitis outbreaks and acute flaccid paralysis surveillance in Shandong province, China during 1988-2020. Phylogenetic and phylogeographic methods based on the VP1 coding region were used to determine the CVA24v origin, spatiotemporal dynamics, and evolution. Also, the positive selection sites in the VP1 gene were identified and exhibited in the tertiary structure.

      Results

      The global CVA24vs were classified into eight genotypes (GⅠ-GⅧ). Here, 12 CVA24v isolates were detected, of which five strains were typed as two novel genotypes (GⅦ and GⅧ) and reported first in the world. The time to the most recent common ancestor of the global CVA24v was estimated around March 1965 and evolved with 5.573 × 10−3 substitutions/site/year. Four residues under positive selection were detected, and residue 146T might be adapted in the CVA24v pandemic. Phylogeographic analysis indicated that China was the main source sink for CVA24v dispersion in a long-lasting global pattern.

      Conclusion

      Our study updated the epidemiological characteristics of CVA24v and enabled a better understanding of the molecular mechanism underlying different genotypes. The results provided new insights into the CVA24v origin, spatiotemporal dynamics, and possibly, the determinants of viral tropism and pathogenicity.

      Graphical abstract

      Keywords

      Introduction

      Acute hemorrhagic conjunctivitis (AHC) is a highly contagious eye disease, with sudden onset of ocular pain, swelling, epiphora, and extensive subconjunctival hemorrhaging as its main manifestations. Enteroviruses (EVs) were implicated as the major causative agents involved in AHC (
      • Akçay E
      • Çarhan A
      • Hondur G
      • Tufan ZK
      • Duru N
      • Kılıç S
      • et al.
      Molecular identification of viral agents associated with acute conjunctivitis: a prospective controlled study.
      ). In particular, the coxsackievirus A24 variant (CVA24v), belonging to the EV-C species within the genus Enterovirus, raised a remarkable concern after its first discovery during an epidemic in Singapore in 1970 (
      • Mirkovic RR
      • Schmidt NJ
      • Yin-Murphy M
      • Melnick JL.
      Enterovirus etiology of the 1970 Singapore epidemic of acute conjunctivitis.
      ) because the virus was recognized as the major etiological agent of AHC with pandemic potential (
      • Enfissi A
      • Joffret ML
      • Delaune D
      • Delpeyroux F
      • Rousset D
      • Bessaud M.
      Coxsackievirus A24 variant associated with acute haemorrhagic conjunctivitis cases, French Guiana, 2017.
      ;
      • Fragoso-Fonseca DE
      • Escobar-Escamilla N
      • Rodríguez-Maldonado AP
      • Barrera-Badillo G
      • Garcés-Ayala F
      • Mendieta-Condado E
      • et al.
      Complete genome sequence of a Coxsackievirus type A24 variant causing an outbreak of acute haemorrhagic conjunctivitis in southeastern Mexico in 2017.
      ;
      • Sousa Jr, IP
      • Burlandy FM
      • Ferreira JL
      • Alves JCS
      • Sousa-Júnior EC
      • Tavares FN
      • et al.
      Re-emergence of a coxsackievirus A24 variant causing acute hemorrhagic conjunctivitis in Brazil from 2017 to 2018.
      ;
      • Zhang L
      • Zhao N
      • Huang X
      • Jin X
      • Geng X
      • Chan TC
      • et al.
      Molecular epidemiology of acute hemorrhagic conjunctivitis caused by coxsackie A type 24 variant in China, 2004–2014.
      ).
      Initially, CVA24v circulated only in Southeast Asia. Thereafter, it spread all over the world, including Africa (
      • Burr SE
      • Sillah A
      • Joof H
      • Bailey RL
      • Holland MJ.
      An outbreak of acute haemorrhagic conjunctivitis associated with coxsackievirus A24 variant in the Gambia, West Africa.
      ), Asia (
      • Zhang L
      • Zhao N
      • Huang X
      • Jin X
      • Geng X
      • Chan TC
      • et al.
      Molecular epidemiology of acute hemorrhagic conjunctivitis caused by coxsackie A type 24 variant in China, 2004–2014.
      ), Europe (
      • Enfissi A
      • Joffret ML
      • Delaune D
      • Delpeyroux F
      • Rousset D
      • Bessaud M.
      Coxsackievirus A24 variant associated with acute haemorrhagic conjunctivitis cases, French Guiana, 2017.
      ), North America (
      • Fragoso-Fonseca DE
      • Escobar-Escamilla N
      • Rodríguez-Maldonado AP
      • Barrera-Badillo G
      • Garcés-Ayala F
      • Mendieta-Condado E
      • et al.
      Complete genome sequence of a Coxsackievirus type A24 variant causing an outbreak of acute haemorrhagic conjunctivitis in southeastern Mexico in 2017.
      ), and South America (
      • Tavares FN
      • Campos Rde M
      • Burlandy FM
      • Fontella R
      • de Melo MM
      • da Costa EV
      • et al.
      Molecular characterization and phylogenetic study of coxsackievirus A24v causing outbreaks of acute hemorrhagic conjunctivitis (AHC) in Brazil.
      ), in the presence of numerous AHC outbreaks and several pandemics. Currently, phylogenetic analysis based on the VP1 coding sequence has been broadly used for typing CVA24v genotypes and shedding light on important topics in CVA24v epidemiology. Previous studies classified CVA24v into four genotypes (GⅠ-GⅣ) (
      • Chu PY
      • Ke GM
      • Chang CH
      • Lin JC
      • Sun CY
      • Huang WL
      • et al.
      Molecular epidemiology of coxsackie A type 24 variant in Taiwan, 2000–2007.
      ;
      • De W
      • Huanying Z
      • Hui L
      • Corina M
      • Xue G
      • Leng L
      • et al.
      Phylogenetic and molecular characterization of coxsackievirus A24 variant isolates from a 2010 acute hemorrhagic conjunctivitis outbreak in Guangdong, China.
      ;
      • Yen YC
      • Chu PH
      • Lu PL
      • Lin YC
      • Shi YY
      • Chou LC
      • et al.
      Phylodynamic characterization of an ocular-tropism coxsackievirus A24 variant.
      ). Nevertheless, GⅠ-GⅢ became silent over the past decades, and GⅣ emerged to be a predominant genotype and was characterized by extremely rapid spread worldwide after 2000 (
      • De W
      • Huanying Z
      • Hui L
      • Corina M
      • Xue G
      • Leng L
      • et al.
      Phylogenetic and molecular characterization of coxsackievirus A24 variant isolates from a 2010 acute hemorrhagic conjunctivitis outbreak in Guangdong, China.
      ;
      • Fonseca MC
      • Pupo-Meriño M
      • García-González LA
      • Resik S
      • Hung LH
      • Muné M
      • et al.
      Molecular evolution of coxsackievirus A24v in Cuba over 23-years, 1986–2009 [sci rep].
      ). Previous studies described some evolutionary aspects of CVA24v and the estimated temporal origin of the virus dating around the 1960s based on VP1 sequence, using simple and robust methods (
      • Chu PY
      • Ke GM
      • Chang CH
      • Lin JC
      • Sun CY
      • Huang WL
      • et al.
      Molecular epidemiology of coxsackie A type 24 variant in Taiwan, 2000–2007.
      ;
      • Fonseca MC
      • Pupo-Meriño M
      • García-González LA
      • Resik S
      • Hung LH
      • Muné M
      • et al.
      Molecular evolution of coxsackievirus A24v in Cuba over 23-years, 1986–2009 [sci rep].
      ;
      • Ishiko H
      • Takeda N
      • Miyamura K
      • Kato N
      • Tanimura M
      • Lin KH
      • et al.
      Phylogenetic analysis of a coxsackievirus A24 variant: the most recent worldwide pandemic was caused by progenies of a virus prevalent around 1981.
      ,
      • Ishiko H
      • Takeda N
      • Miyamura K
      • Tanimura M
      • Yamanaka T
      • Kasuga K
      • et al.
      Phylogenetically different strains of a variant of coxsackievirus A 24 were repeatedly introduced but discontinued circulating in Japan.
      ;
      • Nidaira M
      • Kuba Y
      • Saitoh M
      • Taira K
      • Maeshiro N
      • Mahoe Y
      • et al.
      Molecular evolution of VP3, VP1, 3C(pro) and 3D(pol) coding regions in Coxsackievirus group A type 24 variant isolates from acute hemorrhagic conjunctivitis in 2011 in Okinawa.
      ;
      • Yen YC
      • Chu PH
      • Lu PL
      • Lin YC
      • Shi YY
      • Chou LC
      • et al.
      Phylodynamic characterization of an ocular-tropism coxsackievirus A24 variant.
      ). On the other hand, CVA24v shedding in the stool of patients with AHC has been well-described since the first AHC epidemics as a classical EV; however, this aspect as others was underestimated, despite the fact that several authors had emphasized its importance (
      • Fonseca MC
      • Sarmiento L
      • Resik S
      • Pereda N
      • Rodríguez H
      • Kourí V
      • et al.
      Isolation of Coxsackievirus A24 variant from patients with hemorrhagic conjunctivitis in Cuba, 2008–2009.
      ,
      • Fonseca MC
      • Pupo-Meriño M
      • García-González LA
      • Muné M
      • Resik S
      • Norder H
      • Sarmiento L.
      Molecular characterization of coxsackievirus A24v from feces and conjunctiva reveals epidemiological links.
      ;
      • Lim KH
      • Yin-Murphy M.
      An epidemic of conjunctivitis in Singapore in 1970.
      ;
      • Yin-Murphy M
      • Lim KH
      • Ho YM.
      A Coxsackievirus type A24 epidemic of acute conjunctivitis.
      ). Moreover, identification with comparable prevalence both in ocular and gut (
      • Palacios G
      • Oberste MS.
      Enteroviruses as agents of emerging infectious diseases.
      ) questions the tissue tropism and pathogenicity of CVA24v or at least suggests the rare genetic exchanges with other gutty tropic EVs.
      In Shandong province, China, the acute flaccid paralysis (AFP) surveillance has been conducted for detecting EV circulation since 1988 (
      • Chen P
      • Liu Y
      • Wang H
      • Liu G
      • Lin X
      • Zhang W
      • et al.
      Environmental surveillance complements case-based surveillance of acute flaccid paralysis in polio endgame strategy 2019–2023.
      ). Precisely, our study supported many new findings with a robust analysis in terms of data, covering a long period and considering both CVA24v sequences from patients with AHC and AFP cases, sweeping the entire clinical scenario of this variant. Phylogenetic and phylogeographic analyses were combined to evaluate global CVA24v history, spatiotemporal dynamics, and spreading patterns. The current study shed a new light on selective pressure analysis on the CVA24v VP1 sequence, and sites under positive selection were exhibited in the capsid tertiary structure.

      Methods

      Virus isolation, reverse transcriptase-polymerase chain reaction, and sequencing

      Stool samples were collected from AFP cases at least 24 hours apart and within 14 days of paralysis onset. All samples were transported on ice to the Poliovirus Laboratory at Shandong Center for Disease Control and Prevention, China. After arrival, stool samples were processed according to the
      World Health Organization
      Isolation and identification of polioviruses, WHO polio laboratory manual.
      -recommended procedure, and supernatants were frozen at -20°C and inoculated onto human rhabdomyosarcoma and human epidermoid carcinoma cell lines for virus isolation. As described previously, virus strains isolated from previous AHC outbreaks were also examined (
      • Yang J
      • Lin Y
      • Wang HY
      • Tao ZX
      • Li Y
      • Chen P
      • et al.
      Identification and genetic characterization of coxsackievirus A24 isolated from patients with acute hemorrhagic conjunctivitis in Shandong Province.
      ). All virus isolates were stored at -80°C for further processing. Viral RNA was extracted using a QIAmp viral RNA mini kit according to the manufacturers’ recommended procedures (
      • Chen P
      • Wang H
      • Tao Z
      • Xu A
      • Lin X
      • Zhou N
      • et al.
      Multiple transmission chains of coxsackievirus A4 co-circulating in China and neighboring countries in recent years: phylogenetic and spatiotemporal analyses based on virological surveillance.
      ). The Access Reverse Transcriptase-Polymerase Chain Reaction System (Promega, USA) was used to perform reverse transcriptase-polymerase chain reaction to amplify the complete VP1 sequence (
      • Tao Z
      • Wang H
      • Liu Y
      • Li Y
      • Jiang P
      • Liu G
      • et al.
      Non-polio enteroviruses from acute flaccid paralysis surveillance in Shandong Province, China, 1988–2013.
      ). PCR amplicons were sequenced bidirectionally using the BigDye Terminator v3.0 Cycle Sequencing kit. Molecular typing based on the complete VP1 sequences was performed to identify CVA24v genotypes using online Enterovirus Genotyping Tool v0.1 (
      • Kroneman A
      • Vennema H
      • Deforche K
      • v d Avoort H
      • Peñaranda S
      • Oberste MS
      • et al.
      An automated genotyping tool for enteroviruses and noroviruses.
      ).

      Dataset preparation

      Dataset 1 was designed to include the prototype strain (EH24/70) and 22 complete VP1 sequences of CVA24v obtained in this study. For global CVA24v strains, VP1 sequences were downloaded using the nucleotide program on the National Center for Biotechnology Information website (accessed on April 9, 2022) and annotated with the strain names, GenBank numbers, geographical regions, and collection dates. Before the analysis, we followed the common strategy to reduce the computational load by choosing one closely related sequence that represented sequences with similarity >99% collected at the same location and time point. To deal with unbalanced sequence availability, which could bias the results, two representative VP1 sequence datasets were prepared to alternatively benefit from the higher sequence length or number. Dataset 2 comprised 155 complete VP1 sequences of CVA24v, including 20 Shandong isolates and 135 global strains obtained from 20 countries during 1970-2018. Dataset 3 was based on a partial VP1 region (285 bp), where the higher number of CVA24v sequences had full coverage; altogether, 14 Shandong isolates and 195 global strains obtained from 27 countries were collected. All currently available nucleotide sequences and deduced amino acid sequences of CVA24v were multiply aligned by BioEdit software v7.0.5.3 (
      • Hall TA.
      BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT.
      ).

      Phylogenetic and phylogeographic analyses

      The substitution saturation of all CVA24v sequences was tested with the saturation test of
      • Xia X.
      DAMBE7: New and improved tools for data analysis in molecular biology and evolution.
      , using the software DAMBE v7.0.1, and the root-to-tip regression analysis was conducted to test the temporal signal of the datasets (
      • Rambaut A
      • Lam TT
      • Max Carvalho L
      • Pybus OG
      Exploring the temporal structure of heterochronous sequences using TempEst (formerly Path-O-Gen).
      ). On each of the three datasets, phylogenetic analysis was performed using the maximum likelihood method implemented in MEGA 6 (
      • Tamura K
      • Stecher G
      • Peterson D
      • Filipski A
      • Kumar S.
      MEGA6: molecular evolutionary genetics analysis version 6.0.
      ) after the estimation of the best-fit substitution model with jModelTest v0.1.1 (
      • Posada D.
      JModelTest: phylogenetic model averaging.
      ). For datasets 2 and 3, the tip-dated serial coalescent analyses were performed using the Bayesian Markov chain Monte Carlo method implemented in BEAST 1.10.4 (
      • Suchard MA
      • Lemey P
      • Baele G
      • Ayres DL
      • Drummond AJ
      • Rambaut A.
      Bayesian phylogenetic and phylodynamic data integration using BEAST 1.10.
      ). Multiple combinations of molecular clock models, coalescent models, and discrete trait substitution models were explored and selected by path sampling and stepping-stone approach (
      • Baele G
      • Lemey P
      • Bedford T
      • Rambaut A
      • Suchard MA
      • Alekseyenko AV.
      Improving the accuracy of demographic and molecular clock model comparison while accommodating phylogenetic uncertainty.
      ). The model parameters were estimated with a 60 million generation Markov chain Monte Carlo chain. Run results were assessed by TRACER v1.5. All parameter estimations were summarized as median and 95% highest posterior density (HPD).
      For dataset 2, the continuous time Markov chain process was used to identify the spreading process with the Bayesian stochastic search variable selection model. SPREAD v1.0.6 was used to visualize the results, and a Bayes factor of >10.0 was considered well-supported (
      • Bielejec F
      • Rambaut A
      • Suchard MA
      • Lemey P.
      SPREAD: spatial phylogenetic reconstruction of evolutionary dynamics.
      ). In addition, the maximum clade credibility (MCC) tree was estimated using the TreeAnnotator tool and then visualized with FigTree v1.4.2. A Bayesian skyline plot of CVA24v was generated in TRACER 1.5 by measuring the dynamics of the complete VP1 sequence through time (
      • Drummond AJ
      • Rambaut A
      • Shapiro B
      • Pybus OG.
      Bayesian coalescent inference of past population dynamics from molecular sequences.
      ).

      Selective pressure analysis

      Dataset 2 was selected to estimate positive selection sites using the CODEML algorithm provided in EasyCodeML 1.21 (
      • Gao F
      • Chen C
      • Arab DA
      • Du Z
      • He Y
      • Ho SYW.
      EasyCodeML: a visual tool for analysis of selection using codeml.
      ). The site-specific model was performed to calculate nonsynonymous (dN) and synonymous (dS) substitution rates. Seven substitution models were analyzed and positively selected sites (ω > 1, where ω = dN/dS) were identified in comparisons of M0 (one-ratio) versus M3 (discrete), M1a (nearly neutral) versus M2a (positive selection), M7 (β) versus M8 (β+ωs>1) and M8a (β+ωs=1) versus M8. Among all comparisons, M7 versus M8 is the most stringent test (
      • Anisimova M
      • Bielawski JP
      • Yang Z.
      Accuracy and power of the likelihood ratio test in detecting adaptive molecular evolution.
      ), and representations of the results of the model comparisons were evaluated using likelihood ratio tests (LRTs), and the significance threshold was set to P <0.05. Furthermore, an empirical Bayes method was implemented to predict codon sites under positive selection, with the posterior probability providing good evidence of reliability (
      • Nielsen R
      • Yang Z.
      Likelihood models for detecting positively selected amino acid sites and applications to the HIV-1 envelope gene.
      ;
      • Yang Z
      • Nielsen R
      • Goldman N
      • Pedersen AM.
      Codon-substitution models for heterogeneous selection pressure at amino acid sites.
      ).

      Protein structure prediction

      The tertiary structures of CVA24v VP1 and P1 proteins were generated by a homology model using the Protein Homology/analogY Recognition Engine (Phyre 2.0) (
      • Kelley LA
      • Mezulis S
      • Yates CM
      • Wass MN
      • Sternberg MJ.
      The Phyre2 web portal for protein modeling, prediction and analysis.
      ). Specifically, the homology model of VP1 protein was obtained by submission of the amino acid sequences to the Phyre 2.0. A 100% confidence model covering 90% of the submitted sequence was obtained. The crystal structure of CVA24v (PDB 4Q4V) (
      • Zocher G
      • Mistry N
      • Frank M
      • Hähnlein-Schick I
      • Ekström JO
      • Arnberg N
      • et al.
      A sialic acid binding site in a human picornavirus.
      ) was used as a reference point to predict the P1 structure. All tertiary structures of proteins were visualized, and positively selected sites were mapped using PyMol 2.4.0.

      GenBank accession numbers

      The complete VP1 sequences of 12 CVA24v isolates from AFP cases in this study were deposited in the GenBank database under accession numbers OM386801-OM386812 and JQ28980-JQ28989.

      Results

      Virus isolation and new genotypes identification

      From 1988 to 2020, a total of 10,826 AFP cases were examined for EVs in Shandong province, China. In all, 961 nonpolio EVs were detected, and the overall percentage of nonpolio EV positivity in AFP cases was 8.9%, with the data varying from 15.6% (2003) to 0.4% (2020) in the 33-year surveillance. Of note, 12 isolates, 11 from AFP cases and one from contact, were typed as CVA24v with the molecular method (Table 1). The first isolation of CVA24v can be dated back to 1988, and the last strain was detected in 2010. Also, an outbreak of AHC was reported in Qingdao and Linyi cities in Shandong province in September 2010, and CVA24v was identified as the causative agent, with 10 isolates examined. This was consistent with the detection of CVA24v from AFP surveillance in 2010, suggesting that the case-based AFP surveillance in Shandong province exhibited sufficient sensitivity in detecting CVA24v circulation.
      Table 1Virus characterization and VP1 homology of coxsackievirus A24 variant isolated in Shandong province, China (1988-2020).
      IsolatesCollection dateGenotypeSourceGenBank accession no.VP1 homology with other genotype strains (%)
      880321988Stool/AFP caseOM38680188.594.6-95.092.7-93.788.7-92.190.8-91.492.4-93.0
      94290Sep 1994Stool/AFP caseOM38680987.492.7-93.292.4-93.187.7-91.490.8-91.490.9-91.2
      94405Oct 1994Stool/AFP caseOM38681087.291.9-92.392.0-92.687.6-90.990.2-90.990.7-91.0
      95338Aug 1995Stool/AFP caseOM38681187.592.6-93.192.5-93.287.8-91.590.9-91.591.0-91.3
      97092Apr 1997Stool/AFP caseOM38681286.991.1-91.591.6-92.187.2-90.689.9-90.689.7-90.0
      03353Sep 2003Stool/AFP caseOM38680286.390.1-90.394.5-95.191.8-99.695.3-96.189.3-89.9
      07023Jan 2007Stool/AFP caseOM38680385.688.6-88.891.9-92.292.6-99.792.4-93.188.3-88.8
      07038CFeb 2007Stool/AFP contactOM38680485.688.6-88.891.9-92.292.4-99.792.4-93.188.3-88.8
      10273Sep 2010Stool/AFP caseOM38680586.989.0-89.292.3-92.593.2-99.592.0-92.488.3-88.8
      10276Sep 2010Stool/AFP caseOM38680686.688.7-88.992.2-92.492.8-99.291.9-92.387.7-88.3
      10347Oct 2010Stool/AFP caseOM38680787.289.0-89.292.3-92.593.2-99.592.0-92.488.5-89.0
      10358Nov 2010Stool/AFP caseOM38680886.888.9-89.192.2-92.493.1-99.691.9-92.388.1-88.7
      LY001Sep 2010Conjunctival swab/AHC caseJQ72898086.989.0-89.292.3-92.593.2-99.792.0-92.488.3-88.8
      LY002Sep 2010Conjunctival swab/AHC caseJQ72898186.989.0-89.292.3-92.593.3-99.691.9-92.388.3-88.8
      LY003Sep 2010Conjunctival swab/AHC caseJQ72898287.389.1-89.392.2-92.493.1-99.691.9-92.388.4-88.9
      LY004Sep 2010Conjunctival swab/AHC caseJQ72898386.989.0-89.292.3-92.593.2-99.792.0-92.488.3-88.8
      LY007Sep 2010Conjunctival swab/AHC caseJQ72898486.989.0-89.292.3-92.593.3-99.691.9-92.388.3-88.8
      QD012Sep 2010Conjunctival swab/AHC caseJQ72898587.389.1-89.392.2-92.493.1-99.691.9-92.388.4-88.9
      QD015Sep 2010Conjunctival swab/AHC caseJQ72898686.989.0-89.292.3-92.593.2-99.792.0-92.488.3-88.8
      QD016Sep 2010Conjunctival swab/AHC caseJQ72898787.189.1-89.392.4-92.693.3-99.892.1-92.588.4-88.9
      QD017Sep 2010Conjunctival swab/AHC caseJQ72898886.989.0-89.292.3-92.593.2-99.792.0-92.488.3-88.8
      QD019Sep 2010Conjunctival swab/AHC caseJQ72898987.289.0-89.292.3-92.593.2-99.792.0-92.488.4-88.9
      AFP, Acute flaccid paralysis; AHC, Acute hemorrhagic conjunctivitis
      Of these 22 Shandong CVA24v isolates, the complete VP1 sequences (dataset 1) were used to analyze the genetic characteristics. Shandong CVA24v sequences had 85.6-88.5% (94.4-96.3% amino acid sequence identities) nucleotide acid similarities with the EH24/70 prototype strain (Table 1, Figure 1a). Of note, five CVA24v isolates, 88032, 94290, 94405, 95338, and 97092, showed a relatively remote relationship with the known genotypes. The basic local alignment search tool results indicated that these isolates exhibited the highest homology (97.7%-99.7%) with the Beijing strains (MZ171074-MZ171086) that were detected from patients with AHC in 1988 and 1994. According to the Enterovirus Genotyping Tool, these five isolates were also determined as two novel genotypes of CVA24v. As such, strain 88032 was provisionally assigned to GⅧ and the rest four strains to GⅦ (Table 1). To the best of our knowledge, this is the first report of the GⅦ and GⅧ of CVA24v in the world, and these results proved that not a single CVA24v genotype had existed previously in China and the AFP surveillance in Shandong province was advantageous in identifying virus circulation.
      Figure 1
      Figure 1Homology and phylogenetic relationship of CVA24v viruses isolated from Shandong province, China during 1988-2020. (a) Color-coded pairwise identity matrix generated from 23 CVA24v (22 Shandong CVA24v isolates and the prototype strain) VP1 sequences. Each colored cell represents a percentage identity score between two sequences (one indicated horizontally to the left and the other vertically at the bottom). A colored key indicates the correspondence between pairwise identities and the colors displayed in the matrix. (b) Phylogenetic tree of CVA24v based on the 23 complete VP1 sequences. Triangle indicates the prototype strain.
      Abbreviations in isolate names are as follows: AFP, Acute flaccid paralysis; AHC, Acute hemorrhagic conjunctivitis; SIN, Singapore.
      Abbreviation: CVA24v, coxsackievirus A24 variant.

      Phylogenetic analysis

      Phylogenetic analysis revealed that all Shandong isolates were segregated into three distinct partitions (GⅣ, GⅦ, and GⅧ), with at least 6.2% nucleotide diversity among these three genotypes (Figure 1b). Also, a significant temporal aggregation characteristic could be detected from the phylogenetic tree. To elucidate the genetic characteristics of global CVA24v, phylogenetic trees based on dataset 2 (Figure 2) and dataset 3 (Supplementary material Figure S1) were constructed. CVA24v strains in these two phylogenetic trees had similar topologies, and both temporal and regional aggregation characteristics could be detected. In our study, the global CVA24v strains were segregated into eight distinct partitions (GⅠ-GⅧ), and a limited temporal overlap was observed among these genotypes that circulated worldwide. The phylogenetic analysis showed that GⅠ included only the EH24/70 prototype strain, and GⅡ comprised strains from Brazil, Cuba, and Jamaica during 1986-1987. CVA24v strains in GⅢ were mainly from Cuba (in 1997), Dominican Republic (in 1993), France (in 1999), Philippines (in 2009), and the United States (in 1998). Most of the CVA24v strains were clustered in GⅣ, which included viruses detected since 2002. GⅤ comprised strains from Malaysia, Philippines, and Taiwan province, China, whereas CVA24v strains isolated from Taiwan province, China from 1985 to 1986 formed the separate GⅥ. For the five Shandong CVA24v strains isolated during 1988-1997, four strains were segregated into GⅦ, with strains from Beijing city, China isolated in 1994, and the other one strain isolated in 1988 belonged to GⅧ in the phylogenetic tree (Figure 2, Supplementary material Figure S1).
      Figure 2
      Figure 2(a) Phylogenetic analysis of CVA24v based on 155 complete VP1 sequences from 20 countries by using the maximum likelihood method. ●, Viruses isolated in conjunctival/eye/corneal swab. ○, Viruses isolated in other specimens (nasopharyngeal swab, stool or saliva sample), ■, Viral sources are acute hemorrhagic conjunctivitis cases. □, Viral sources are other populations (acute flaccid paralysis, Hopkins syndrome or normal person). Specimens and sources with blank indicate unknown. Diamonds indicate isolates detected in this study, and triangle indicates the prototype strain. Abbreviations of geographic locations are as follows: AUS, Australia; BRA, Brazil; DOM, Dominican Republic; ESP, Spain; GAB, Gabon; GUF, French Guiana; IND, India; JAM, Jamaica; JPN, Japan; KEN, Kenya; KOR, Korea; MAS, Malaysia; MEX, Mexico; REF, French Reunion Island; SIN, Singapore; THA, Thailand; UGA, Uganda; USA, United States of America; and abbreviations of GD-CHN, JS-CHN, TW-CHN, YN-CHN, ZJ-CHN represent the provinces of Guangdong, Jiangsu, Taiwan, Yunnan, Zhejiang in China. (b) Mutated amino acid residues detected in the 137 VP1 protein alignments according to the prototype sequence. Abbreviations of residues are as follows: A, Alanine; C, Cysteine; D, Aspartic acid; E, Glutamic acid; F, Phenylalanine; G, Glycine; H, Histidine; I, Isoleucine; K, Lysine; L, Leucine; M, Methionine; N, Asparagine; P, Proline; Q, Glutamine; R, Arginine; S, Serine; T, Threonine; V, Valine; Y, Tyrosine.
      Abbreviation: CVA24v, coxsackievirus A24 variant.
      By inspecting the VP1 protein alignment in dataset 2, a total of 81 amino acids were found to mutate according to the sequence of the prototype strain (Figure 2b). Of these, 15 highly variable amino acid residues possessed the specificity of the mutation in genotypes and thus were analyzed in detail. All CVA24v strains shared nine common mutations (S11T, V51A, I56V, M89I, Y151H, I196M, F250Y, I256T, and I297T) in VP1 protein from the prototype strain (GⅠ), and the strains belonging to GⅡ-GⅥ and GⅧ also had T146A mutation. In addition, amino acid residues at positions 25 and 32 were also specific to different genotypes. GⅡ, GⅢ, GⅦ, and GⅧ possessed L25S mutation, GⅥ and GⅤ strains possessed L25T and L25P mutation, respectively, whereas GⅣ had both the L25P and L25H mutations. Apart from GⅣ strains, which possessed S32L, the other CVA24v strains possessed S32P mutation in VP1. Furthermore, several specific amino acid mutations were also identified in the VP1 protein of some genotype strains. For example, the E100D mutation was found in GⅣ and GⅤ, whereas the K103R mutation was specific to GⅣ-GⅥ. Similar to amino acid mutation at position 25, the N301T mutation was detected in GⅥ, whereas GⅣ possessed both the N301S and N301D mutations based on the current data.

      Coxsackievirus A24 variant origin, evolutionary rate, and population dynamics

      The analysis performed on dataset 2 provided a time to the most recent common ancestor (tMRCA) of 53.27 years (95% HPD: 50.21-56.70), which corresponds to the original date of March 1965 (95% HPD: October 1961-March 1967) (Table 2). The molecular evolution rate of global CVA24v was inferred as 5.573 × 10−3 substitutions/site/year, with the 95% HPD of 4.858 × 10−3-6.341 × 10−3 substitutions/site/year (Table 2). Moreover, the substitution rate of GⅣ strains was estimated as 5.396 × 10−3 substitutions/site/year (95% HPD: 4.481 × 10−3-6.377 × 10−3 substitutions/site/year), which was very homogeneous with the global mean evolutionary rate. We also estimated the tMRCA and evolutionary rate of global CVA24v based on dataset 3 (Supplementary material Table S1). All estimated parameters were very homogeneous with dataset 2, and the 95% HPD interval of the estimates did not allow the rejection of the hypothesis that the parameters estimated on the VP1 sequence are the same between datasets 2 and 3.
      Table 2Estimates of the tMRCA and ancestral origin of coxsackievirus A24 variant on the complete VP1 gene.
      GenotypetMRCA (year)Ancestral origin
      Median value95% highest posterior densityEffective sample size valueLocationPosterior probability
      Global virusMar 1965Sep 1961, Mar 19681472.8Singapore1.0
      Dec 1986Jun 1986, Apr 19875256.7Brazil1.0
      Jul 1992Jul 1991, Feb 19932681.5Dominican Republic/ Taiwan province, China1.0/1.0
      Aug 2001Jan 2001, Jan 2002726.9Taiwan province, China1.0
      Feb 2000Sep 1999, Jun 200012651.8Taiwan province, China1.0
      Jan 1985Jul 1984, Apr 19854266.2Taiwan province, China1.0
      Apr 1993Sep 1992, Aug 19931082.6Shandong province, China0.9994
      Sep 1987Apr 1987, Feb 19883972.5Shandong province, China1.0
      tMRCA, the most recent common ancestor.
      The evolutionary epidemiology of global CVA24v was performed based on dataset 2 (Figure 3). The skyline plot clearly showed that the relative genetic diversity presented nearly constant before 2000; however, a smooth and steady rise was observed between 2001 and 2003, indicating the emergence of GⅣ on a global scale (Figure 3a and Table 2). A sharp but transient increase in relative genetic diversity was also observed for the CVA24v between 2010 and 2011, peaking approximately at the turn of the year. This trend coincided with the evolutionary epidemiological characteristics of GⅣ in Figure 3b. In addition, the global CVA24v data were summarized by the corresponding location and collection date (Figure 3c and 3d). Clearly, the genetic diversity trend showed an evident consistency with the global CVA24v data, suggesting that elevation in genetic diversity was consistent with the epidemiological data on the large scale.
      Figure 3
      Figure 3Demographic history of global CVA24v through time. Bayesian skyline plot estimates the relative genetic diversity of global CVA24v (a) and GⅣ strains (b) through time. The solid line stands for the median and shaded areas correspond to the 95% HPD intervals. The global reporting of CVA24v isolation from 1970 to 2020 was analyzed, compiled from the published GenBank database (c). The circles indicate CVA24v is detected in the corresponding country and year, and different colors are used to distinguish virus from different continents. We also depict the number of countries detected CVA24v per year in d. Abbreviations of geographic locations are the same as in Figure 2.
      Abbreviation: CVA24v, coxsackievirus A24 variant.

      Phylogeographic analysis

      The phylogenetic topology of MCC tree for global CVA24v was highly congruent with the maximum likelihood trees in genotype grouping (Figure 4). Ancestral CVA24v initially diverged into eight genotypes (GⅠ-GⅧ) during 1964-1997 (posterior probability ranging from 0.9994 to 1.0). The MCC tree showed that all these genotypes, except the GⅥ, were relatively conservative due to the limited sequences, whereas GⅥ strains were more active as imported over short or similar timescales and then co-circulated (Figure 4b). Except for the prototype strain, the most probable ancestral geographical location for all contemporary CVA24v strains was found to be Taiwan province, China, with posterior probability 1.0 (Supplementary material Figure S2). However, the prototype strain, which originated from Singapore, had been extinct since the last detection in 1970. The ancestral state for GⅢ-GⅥ was also found to be from Taiwan province, China (Table 2), inferring Taiwan province, China as the main source sink for the dispersion of the virus. The alignment of amino acid sequences, with the prototype strain EH24/70 as a reference strain, suggested that CVA24v strains harbored several amino acid mutations throughout the historical background (Figure 4b). Among these variable amino acid sites, 14 residues were predicted to substitute in about 38.7 years ago. In addition, amino acids at positions 25, 32, 103, 146, and 301 showed distinct features of dynamic evolution over time and were consistent with the temporal dynamics of CVA24v genotypes (Figure 4c).
      Figure 4
      Figure 4Phylogenetic and molecular clock analysis of global CVA24v. (a) Root-to-tip regression analysis estimated with a maximum likelihood tree. R2, coefficient of determination. (b) MCC tree of complete VP1 sequences of global CVA24v visualized in FigTree. The colors of the branches correspond to the variable amino acid residues, and colors of the strain names correspond to their geographic locations. The width of a branch reflects the evolution rate of individual sequences and their reconstructed ancestors. The node labels depict the median of the tMRCA estimates. Diamonds indicate isolates detected in this study, and triangle indicates the prototype strain. (c) An illustration of the temporal amino acid evolution of the CVA24v VP1 gene. Amino acid changes are inferred using the EH24/70 as the parental strain. Observed amino acid substitutions are plotted with their estimated times of occurrence, as obtained from the MCC tree in panel B. Abbreviations of amino acid residues and geographic locations are the same as in Figure 2.
      Abbreviation: CVA24v, coxsackievirus A24 variant; MCC, maximum clade credibility.
      Phylogeographic reconstructions of the CVA24v spread over time were estimated in terms of well-supported (Bayes factor >16 and posterior probability >0.5) contacts among locations (Figure 5). The Bayesian stochastic search variable selection analysis supported the presence of four interconnected nuclei of viral transmission, corresponding to Africa, Asia, North America, and South America, with Africa and Asia acting as a central role in the dispersal of CVA24v. These regions showed significant transmission links with each other, irrespective of their geographical proximity (Figure 5a and 5b). Within these main areas, three significant geographic migrations were observed from Gabon to Brazil, from Brazil to Taiwan province, China, and from Japan to Uganda. Interestingly, Gabon was at the core of the CVA24v diffusion process as being the most likely responsible for viral spread to other locations (Figure 5c).
      Figure 5
      Figure 5Geographic distribution and inferred spatial dynamics of CVA24v. (a) Map showing the diffusion process inferred using SPREAD software and visualized in Google™ Earth (http://earth.google.com). Arrows between locations indicate the inferred migration routes of CVA24v, and colors are corresponding to both the BF values in the legend and the posterior probability in panel B. (b) The posterior probability calculated for the migration route pairs. (c) Histogram of the total number of location state transitions.
      Abbreviations: BF, Bayes Factor; CVA24v, coxsackievirus A24 variant The geographic abbreviations are the same as in Figure 2.

      Selective pressure analysis

      Selective pressure analysis on the codon sites of the VP1 gene of CVA24v was investigated, and the log likelihood values and their LRTs are shown in Table 3. Among these models, M0 estimated an average ω across all codon sites of 0.08673, which indicated evolution by strong purifying selection. For the other two models, which also permitted an ω>1 (M2 and M8), there was a less degree of overlap in the positively selected amino acid sites in CVA24v. Moreover, LRTs indicated a strong evidence of positive selection on M7 versus M8 model comparison, and one site (146T) was found to potentially be under positive selection, with a posterior probability greater than 0.95.
      Table 3Parameter estimates and log likelihood values for coxsackievirus A24 variant under different models of variable ω ratios among VP1 sites.
      ModelParameter estimatesLikelihood scoresModel comparedLikelihood ratio test P- valuePositively selected sites
      Sites with a posterior probability >50% of having ω>1. Positively selected sites with posterior probability ≥90 are underlined, 0.7-0.8 in bold, 0.5-0.7 in plain text. Model detected posterior probabilities based on Bayes Empirical Bayes analysis. Letters refer to the amino acids residues found in the prototype strain EH24/70 (D90457).
      M3: Discreteω0=0.04216, ω1=0.64265, ω2=34.69808, (p0=0.91849, p1=0.08151, p2<0.00001)-6043.748402M0 vs M3< 0.000000001-
      M0: One-ratioω0=0.08673-6096.656498Not Allowed
      M2a: Positive selectionω0=0.05062, ω1=1.00000, ω2=22.60055, (p0=0.94446, p1=0.05554, p2<0.00001)-6046.645442M1a vs M2a0.999175340146T, 0.779
      M1a: Nearly neutralω0=0.05062, ω1=1.00000, (p0=0.94446, p1=0.05554)-6046.644617Not Allowed
      M8: β+ωs>1p0=0.96256, p1= 0.03744, p=0.02919, q=0.19408, ω= 1.00000-6046.948267M7 vs M80.00003788925L, 0.716; 146T, 0.976; 269R, 0.518; 301N, 0.591
      M7: βp=0.23839, q=1.76052-6057.129118Not Allowed
      M8a: β+ωs=1p0=0.96237, p1= 0.03763, p=0.58217, q=8.92196, ω= 1.00000-6043.004368M8a vs M80.004976942Not Allowed
      M0: one-ratio ω value is the average for all codon sites, whereas M2a, M3 and M8 models ω values are the estimated values for the positively selected codon sites under respective models.
      a Sites with a posterior probability >50% of having ω>1. Positively selected sites with posterior probability ≥90 are underlined, 0.7-0.8 in bold, 0.5-0.7 in plain text. Model detected posterior probabilities based on Bayes Empirical Bayes analysis. Letters refer to the amino acids residues found in the prototype strain EH24/70 (D90457).

      Protein structure prediction

      The three-dimensional structures of both VP1 (Figure 6a and 6b) and P1 (Figure 6c and 6d) genes for CVA24v were successfully predicted with high confidence, and 277 (91%) and 858 (97%) residues modeled at >90% confidence for VP1 and P1, respectively. From the predicted structures, the number of the two common secondary structures for α-helix and β-sheet existed in the VP1 gene were six (11.8% residues) and nine (32.1% residues), and in the P1 gene were 21 (15.8% residues) and 38 (26.8% residues), both of which covered one positively selected site (at position 269), as seen in Table 3. Evidently, all four positively selected residues were located on the capsid protein surface (Figure 6b and 6d). We finally predicted the residue preference and mutation effect on the four residues under positive selection (Supplementary material Figure S3). Compared with the amino acid residues at positions 25, 146, and 301, the residue at position 269 showed a relatively high mutational effect, especially as the site was associated with either a D, G, P, or V mutation.
      Figure 6
      Figure 6Tertiary structure of the VP1 protein (a and b) and P1 protein (c and d) estimated with homology model. The α-helix and β-sheet existed in the protein structure are highlighted in blue and yellow, respectively. When located on the protein surfaces, sites under positive selection are highlighted in red. Letters and numbers refer to the amino acids residues and positions in the prototype strain EH24/70 (D90457).

      Discussion

      CVA24v was demonstrated to circulate in four different periods, and four genotypes had been detected according to previous studies (
      • Oh MD
      • Park S
      • Choi Y
      • Kim H
      • Lee K
      • Park W
      • et al.
      Acute hemorrhagic conjunctivitis caused by coxsackievirus A24 variant, South Korea, 2002.
      ;
      • Wu B
      • Qi X
      • Xu K
      • Ji H
      • Zhu Y
      • Tang F
      • et al.
      Genetic characteristics of the coxsackievirus A24 variant causing outbreaks of acute hemorrhagic conjunctivitis in Jiangsu, China, 2010.
      ). In this study, however, our results provided a more comprehensive picture of CVA24v, and at least eight genotypes were detected with several genotypes co-circulating simultaneously at times in the past. Compared with most previous phylogenetic analyses on the 3Cpro region (
      • De W
      • Huanying Z
      • Hui L
      • Corina M
      • Xue G
      • Leng L
      • et al.
      Phylogenetic and molecular characterization of coxsackievirus A24 variant isolates from a 2010 acute hemorrhagic conjunctivitis outbreak in Guangdong, China.
      ;
      • Fonseca MC
      • Pupo-Meriño M
      • García-González LA
      • Resik S
      • Hung LH
      • Muné M
      • et al.
      Molecular evolution of coxsackievirus A24v in Cuba over 23-years, 1986–2009 [sci rep].
      ;
      • Yen YC
      • Chu PH
      • Lu PL
      • Lin YC
      • Shi YY
      • Chou LC
      • et al.
      Phylodynamic characterization of an ocular-tropism coxsackievirus A24 variant.
      ), our study initially depicted GⅦ and GⅧ of CVA24v, which included strains from Beijing city (1988 and 1994) and Shandong province (1988–1997), China. To the best of our knowledge, this is the first report of GⅦ and GⅧ of CVA24v in the world, and our study reveals that at least three genotypes had existed previously in mainland China. This cyclical circulation of the pandemic genotypes could be one distinctive characteristic of CVA24v, as described previously (
      • Chu PY
      • Ke GM
      • Chang CH
      • Lin JC
      • Sun CY
      • Huang WL
      • et al.
      Molecular epidemiology of coxsackie A type 24 variant in Taiwan, 2000–2007.
      ;
      • De W
      • Huanying Z
      • Hui L
      • Corina M
      • Xue G
      • Leng L
      • et al.
      Phylogenetic and molecular characterization of coxsackievirus A24 variant isolates from a 2010 acute hemorrhagic conjunctivitis outbreak in Guangdong, China.
      ;
      • Fonseca MC
      • Pupo-Meriño M
      • García-González LA
      • Resik S
      • Hung LH
      • Muné M
      • et al.
      Molecular evolution of coxsackievirus A24v in Cuba over 23-years, 1986–2009 [sci rep].
      ). Notably, the detection of GⅦ and GⅧ strains from AFP cases revealed that this variant possessed a potential neuroinvasion capability besides causing AHC, and the case-based AFP surveillance in Shandong province exhibited sufficient sensitivity in detecting CVA24v circulation. Although at least 75% VP1 nucleotide identity was recommended as the prevailing standard to classify EV serotypes (
      • Oberste MS
      • Schnurr D
      • Maher K
      • Al-Busaidy S
      • Pallansch MA.
      Molecular identification of new picornaviruses and characterization of a proposed enterovirus 73 serotype.
      ), no nucleotide identity standard was established to define the genotypes of CVA24v. Generally, CVA24v genotypes were defined based on the phylogenetic relationships obtained in the dendrograms, which occasionally varied depending on the selected sequences. Consequently, we recommended revising the CVA24v genotypes proposed in accordance with the previous phylogenies. In our study, the VP1 nucleotide divergences between different CVA24v genotypes were not higher than 9%. Hence, further studies are needed, and a consensus is necessary to obtain a better classification of CVA24v; moreover, new genotypes should be proposed according to chronological order.
      As is well known, CVA24v is an ocular tropism pathogen and is associated with AHC. Hypothetically, if a virus shows the characteristics of the preferential targeting of specific tissue, a special growth environment may have a great impact on the variation of this virus. Previously, several studies have meritoriously attempted to investigate the CVA24v origin and dynamic evolution, and the preliminary results suggested a recent CVA24v emergence, which was similar to the isolation year of the prototype strain (
      • Fonseca MC
      • Pupo-Meriño M
      • García-González LA
      • Resik S
      • Hung LH
      • Muné M
      • et al.
      Molecular evolution of coxsackievirus A24v in Cuba over 23-years, 1986–2009 [sci rep].
      ;
      • Nidaira M
      • Kuba Y
      • Saitoh M
      • Taira K
      • Maeshiro N
      • Mahoe Y
      • et al.
      Molecular evolution of VP3, VP1, 3C(pro) and 3D(pol) coding regions in Coxsackievirus group A type 24 variant isolates from acute hemorrhagic conjunctivitis in 2011 in Okinawa.
      ;
      • Yen YC
      • Chu PH
      • Lu PL
      • Lin YC
      • Shi YY
      • Chou LC
      • et al.
      Phylodynamic characterization of an ocular-tropism coxsackievirus A24 variant.
      ). Based on these premises, the current study reconstructed the evolutionary phylogeography of CVA24v based on a high quality and updated sequence dataset, spanning a sampling time of nearly 50 years (1970-2018). Our study supported a more ancient origin but a lower substitution rate of CVA24v. Our analysis of sequence divergence in the VP1 region showed a slightly lower rate of sequence change than what was published for CVB1 (7.0 × 10−3 substitutions/site/year), echovirus 14 (9.2 × 10−3 substitutions/site/year), echovirus 6 (7.0 × 10−3 substitutions/site/year), and EV-A71 (7.2 × 10−3 substitutions/site/year) (
      • Abdelkhalek I
      • Seghier M
      • Yahia AB
      • Touzi H
      • Meddeb Z
      • Triki H
      • et al.
      Molecular epidemiology of Coxsackievirus type B1.
      ;
      • Chen P
      • Li Y
      • Tao Z
      • Wang H
      • Lin X
      • Liu Y
      • et al.
      Evolutionary phylogeography and transmission pattern of echovirus 14: an exploration of spatiotemporal dynamics based on the 26-year acute flaccid paralysis surveillance in Shandong, China.
      ;
      • Kyriakopoulou Z
      • Bletsa M
      • Tsakogiannis D
      • Dimitriou TG
      • Amoutzias GD
      • Gartzonika C
      • et al.
      Molecular epidemiology and evolutionary dynamics of echovirus 3 serotype.
      ;
      • McWilliam Leitch EC
      • Cabrerizo M
      • Cardosa J
      • Harvala H
      • Ivanova OE
      • Koike S
      • et al.
      The association of recombination events in the founding and emergence of subgenogroup evolutionary lineages of human enterovirus 71.
      ;
      • Tao Z
      • Wang H
      • Li Y
      • Xu A
      • Zhang Y
      • Song L
      • et al.
      Cocirculation of two transmission lineages of echovirus 6 in jinan, china, as revealed by environmental surveillance and sequence analysis.
      ), and all results turned out to be 20% (about 1.0 × 10−3 substitutions/site/year) higher than CVA24v, suggesting special growth environment may limit the variation of this virus. Several studies also pointed out the occurrence of a “time-dependent rate phenomenon” in the virus evolutionary process, as the evolutionary rates appeared to vary over time and continuously decreased with the timescale of rate measurement (
      • Aiewsakun P
      • Katzourakis A.
      Time-dependent rate phenomenon in viruses.
      ;
      • Holmes EC.
      What does virus evolution tell us about virus origins?.
      ).
      A particular focus of this study was to determine the phylogeographic relationship of global CVA24v. This uncontrolled viral circulation could be explained by the consistent CVA24v identification and by its temporal aggregation characteristics (
      • Enfissi A
      • Joffret ML
      • Delaune D
      • Delpeyroux F
      • Rousset D
      • Bessaud M.
      Coxsackievirus A24 variant associated with acute haemorrhagic conjunctivitis cases, French Guiana, 2017.
      ;
      • Fragoso-Fonseca DE
      • Escobar-Escamilla N
      • Rodríguez-Maldonado AP
      • Barrera-Badillo G
      • Garcés-Ayala F
      • Mendieta-Condado E
      • et al.
      Complete genome sequence of a Coxsackievirus type A24 variant causing an outbreak of acute haemorrhagic conjunctivitis in southeastern Mexico in 2017.
      ;
      • Sousa Jr, IP
      • Burlandy FM
      • Ferreira JL
      • Alves JCS
      • Sousa-Júnior EC
      • Tavares FN
      • et al.
      Re-emergence of a coxsackievirus A24 variant causing acute hemorrhagic conjunctivitis in Brazil from 2017 to 2018.
      ). In particular, CVA24v from four countries seemed to have played a central role in spreading worldwide, with Africa and Asia acting important roles in the dispersal of CVA24v. The similar spatiotemporal dynamics of CVA24v could also be detected in the earlier studies with Asia (especially China) playing a key role in spreading the virus worldwide (
      • Fonseca MC
      • Pupo-Meriño M
      • García-González LA
      • Resik S
      • Hung LH
      • Muné M
      • et al.
      Molecular evolution of coxsackievirus A24v in Cuba over 23-years, 1986–2009 [sci rep].
      ;
      • Yen YC
      • Chu PH
      • Lu PL
      • Lin YC
      • Shi YY
      • Chou LC
      • et al.
      Phylodynamic characterization of an ocular-tropism coxsackievirus A24 variant.
      ). We believe that large-scale patterns in human mobility between these locations will provide a more precise explanation on CVA24v emergence and introduction.
      Despite at least eight genotypes having been reported in the last 50-year CVA24v circulation, two increases in relative genetic diversity related to GⅣ were observed using a skyline plot. As mentioned previously, the ancient origin of CVA24v also implies a lower evolutionary rate than other EV pathogens. All this evidence could suggest a lower intensity of diversifying selective pressure shaping the evolution of the virus. The research also found that biodiversity loss of amino acids was tightly linked to the emergence of the AHC-causing CVA24v, probably because of viral tropism or a population bottleneck (
      • Baggen J
      • Hurdiss DL
      • Zocher G
      • Mistry N
      • Roberts RW
      • Slager JJ
      • et al.
      Role of enhanced receptor engagement in the evolution of a pandemic acute hemorrhagic conjunctivitis virus.
      ;
      • Smura T
      • Blomqvist S
      • Vuorinen T
      • Ivanova O
      • Samoilovich E
      • Al-Hello H
      • et al.
      The evolution of Vp1 gene in enterovirus C species sub-group that contains types CVA-21, CVA-24, EV-C95, EV-C96 and EV-C99.
      ). Thus, the identification of positive selection residues in CVA24v VP1 may provide a reference for future studies to investigate whether changes in the capsid coding region have contributed to the ocular tropism manifestation of AHC. These results, coupled with the location on the capsid protein surface of sites under positive selection, suggested at least a certain role of the host immune response in shaping CVA24v evolution. Although some studies had attempted to map the CVA24v receptor in in vitro models and investigate the virus immunopathogenesis and interaction within the host (
      • Baggen J
      • Hurdiss DL
      • Zocher G
      • Mistry N
      • Roberts RW
      • Slager JJ
      • et al.
      Role of enhanced receptor engagement in the evolution of a pandemic acute hemorrhagic conjunctivitis virus.
      ), the field of CVA24v immunology was still in its infancy, and further confirmation should be needed to determine what had facilitated the initial adaptation of CVA24v to the eye, which was likely the most relevant determinants of its evolution.

      Conclusion

      The global CVA24v was divided into eight genotypes, and five genotypes circulated previously in China. The current study reported two novel CVA24v genotypes (GⅦ and GⅧ) first and brought new insights into the virus epidemiological understanding and phylogenetic relationship. The phylogeographic analysis provided an updated representation of CVA24v origin, spatiotemporal dynamics, and evolution, suggesting that special ocular tropism might impact the variation of CVA24v. A complex and nondirectional viral flow network was evidenced, indicating a long-lasting transmission pattern of CVA24v and with China being the main source sink for the global dispersion. Combined with the selective pressure analysis, all these results could contribute to the evaluation of the actual CVA24v relevance for AHC outbreaks, and possibly, to the determinants of viral tropism and pathogenicity.

      Declarations of competing interest

      The authors have no competing interests to declare.

      Funding

      This work was supported by grants from the National Natural Science Foundation of China (grant numbers 82003510, 81302481), the Shandong provincial Natural Science Foundation, China (grant number ZR2019BH085), the Taishan Scholar Program for Young Experts (tsqn202103187), and the Medicine and Health Science Technology Development Program of Shandong province (202012061240). None of the funders had any influence over the study design, analysis, or decision to submit for publication.

      Ethical approval

      Ethical approval was given by the ethics review committee of Shandong Center for Disease Control and Prevention, and the methods were carried out in accordance with the principles of the Declaration of Helsinki. Written informed consents for the use of clinical samples of acute flaccid paralysis cases were obtained from all subjects (the legal guardians of the patients and contacts).

      Authors’ contributions

      AX conceived the study design and critically reviewed the manuscript. PC contributed to the study idea, virological investigation, data analysis and writing of the paper. ZT contributed to the study design and provided management support of the epidemiological data. XL, SW and FJ contributed to laboratory investigation and participated in data interpretation. YL (Yan Li) and YL (Yao Liu) revised the manuscript. All authors read and approved the final manuscript.

      References

        • Abdelkhalek I
        • Seghier M
        • Yahia AB
        • Touzi H
        • Meddeb Z
        • Triki H
        • et al.
        Molecular epidemiology of Coxsackievirus type B1.
        Arch Virol. 2015; 160: 2815-2821
        • Aiewsakun P
        • Katzourakis A.
        Time-dependent rate phenomenon in viruses.
        J Virol. 2016; 90: 7184-7195
        • Akçay E
        • Çarhan A
        • Hondur G
        • Tufan ZK
        • Duru N
        • Kılıç S
        • et al.
        Molecular identification of viral agents associated with acute conjunctivitis: a prospective controlled study.
        Braz J Infect Dis. 2017; 21: 391-395
        • Anisimova M
        • Bielawski JP
        • Yang Z.
        Accuracy and power of the likelihood ratio test in detecting adaptive molecular evolution.
        Mol Biol Evol. 2001; 18: 1585-1592
        • Baele G
        • Lemey P
        • Bedford T
        • Rambaut A
        • Suchard MA
        • Alekseyenko AV.
        Improving the accuracy of demographic and molecular clock model comparison while accommodating phylogenetic uncertainty.
        Mol Biol Evol. 2012; 29: 2157-2167
        • Baggen J
        • Hurdiss DL
        • Zocher G
        • Mistry N
        • Roberts RW
        • Slager JJ
        • et al.
        Role of enhanced receptor engagement in the evolution of a pandemic acute hemorrhagic conjunctivitis virus.
        Proc Natl Acad Sci U S A. 2018; 115: 397-402
        • Bielejec F
        • Rambaut A
        • Suchard MA
        • Lemey P.
        SPREAD: spatial phylogenetic reconstruction of evolutionary dynamics.
        Bioinformatics. 2011; 27: 2910-2912
        • Burr SE
        • Sillah A
        • Joof H
        • Bailey RL
        • Holland MJ.
        An outbreak of acute haemorrhagic conjunctivitis associated with coxsackievirus A24 variant in the Gambia, West Africa.
        BMC Res Notes. 2017; 10: 692
        • Chen P
        • Li Y
        • Tao Z
        • Wang H
        • Lin X
        • Liu Y
        • et al.
        Evolutionary phylogeography and transmission pattern of echovirus 14: an exploration of spatiotemporal dynamics based on the 26-year acute flaccid paralysis surveillance in Shandong, China.
        BMC Genomics. 2017; 18: 48
        • Chen P
        • Liu Y
        • Wang H
        • Liu G
        • Lin X
        • Zhang W
        • et al.
        Environmental surveillance complements case-based surveillance of acute flaccid paralysis in polio endgame strategy 2019–2023.
        Appl Environ Microbiol. 2020; 86: e00702-e00720
        • Chen P
        • Wang H
        • Tao Z
        • Xu A
        • Lin X
        • Zhou N
        • et al.
        Multiple transmission chains of coxsackievirus A4 co-circulating in China and neighboring countries in recent years: phylogenetic and spatiotemporal analyses based on virological surveillance.
        Mol Phylogenet Evol. 2018; 118: 23-31
        • Chu PY
        • Ke GM
        • Chang CH
        • Lin JC
        • Sun CY
        • Huang WL
        • et al.
        Molecular epidemiology of coxsackie A type 24 variant in Taiwan, 2000–2007.
        J Clin Virol. 2009; 45: 285-291
        • De W
        • Huanying Z
        • Hui L
        • Corina M
        • Xue G
        • Leng L
        • et al.
        Phylogenetic and molecular characterization of coxsackievirus A24 variant isolates from a 2010 acute hemorrhagic conjunctivitis outbreak in Guangdong, China.
        Virol J. 2012; 9: 41
        • Drummond AJ
        • Rambaut A
        • Shapiro B
        • Pybus OG.
        Bayesian coalescent inference of past population dynamics from molecular sequences.
        Mol Biol Evol. 2005; 22: 1185-1192
        • Enfissi A
        • Joffret ML
        • Delaune D
        • Delpeyroux F
        • Rousset D
        • Bessaud M.
        Coxsackievirus A24 variant associated with acute haemorrhagic conjunctivitis cases, French Guiana, 2017.
        Intervirology. 2017; 60: 271-275
        • Fonseca MC
        • Sarmiento L
        • Resik S
        • Pereda N
        • Rodríguez H
        • Kourí V
        • et al.
        Isolation of Coxsackievirus A24 variant from patients with hemorrhagic conjunctivitis in Cuba, 2008–2009.
        J Clin Virol. 2012; 53: 77-81
        • Fonseca MC
        • Pupo-Meriño M
        • García-González LA
        • Resik S
        • Hung LH
        • Muné M
        • et al.
        Molecular evolution of coxsackievirus A24v in Cuba over 23-years, 1986–2009 [sci rep].
        Sci Rep. 2020; 10: 13761
        • Fonseca MC
        • Pupo-Meriño M
        • García-González LA
        • Muné M
        • Resik S
        • Norder H
        • Sarmiento L.
        Molecular characterization of coxsackievirus A24v from feces and conjunctiva reveals epidemiological links.
        Microorganisms. 2021; 9: 531
        • Fragoso-Fonseca DE
        • Escobar-Escamilla N
        • Rodríguez-Maldonado AP
        • Barrera-Badillo G
        • Garcés-Ayala F
        • Mendieta-Condado E
        • et al.
        Complete genome sequence of a Coxsackievirus type A24 variant causing an outbreak of acute haemorrhagic conjunctivitis in southeastern Mexico in 2017.
        Arch Virol. 2020; 165: 1015-1018
        • Gao F
        • Chen C
        • Arab DA
        • Du Z
        • He Y
        • Ho SYW.
        EasyCodeML: a visual tool for analysis of selection using codeml.
        Ecol Evol. 2019; 9: 3891-3898
        • Hall TA.
        BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT.
        Nucl Acids Symp Ser. 1999; 41: 95-98
        • Holmes EC.
        What does virus evolution tell us about virus origins?.
        J Virol. 2011; 85: 5247-5251
        • Ishiko H
        • Takeda N
        • Miyamura K
        • Kato N
        • Tanimura M
        • Lin KH
        • et al.
        Phylogenetic analysis of a coxsackievirus A24 variant: the most recent worldwide pandemic was caused by progenies of a virus prevalent around 1981.
        Virology. 1992; 187: 748-759
        • Ishiko H
        • Takeda N
        • Miyamura K
        • Tanimura M
        • Yamanaka T
        • Kasuga K
        • et al.
        Phylogenetically different strains of a variant of coxsackievirus A 24 were repeatedly introduced but discontinued circulating in Japan.
        Arch Virol. 1992; 126: 179-193
        • Kelley LA
        • Mezulis S
        • Yates CM
        • Wass MN
        • Sternberg MJ.
        The Phyre2 web portal for protein modeling, prediction and analysis.
        Nat Protoc. 2015; 10: 845-858
        • Kroneman A
        • Vennema H
        • Deforche K
        • v d Avoort H
        • Peñaranda S
        • Oberste MS
        • et al.
        An automated genotyping tool for enteroviruses and noroviruses.
        J Clin Virol. 2011; 51: 121-125
        • Kyriakopoulou Z
        • Bletsa M
        • Tsakogiannis D
        • Dimitriou TG
        • Amoutzias GD
        • Gartzonika C
        • et al.
        Molecular epidemiology and evolutionary dynamics of echovirus 3 serotype.
        Infect Genet Evol. 2015; 32: 305-312
        • Lim KH
        • Yin-Murphy M.
        An epidemic of conjunctivitis in Singapore in 1970.
        Singapore Med J. 1971; 12: 247-249
        • McWilliam Leitch EC
        • Cabrerizo M
        • Cardosa J
        • Harvala H
        • Ivanova OE
        • Koike S
        • et al.
        The association of recombination events in the founding and emergence of subgenogroup evolutionary lineages of human enterovirus 71.
        J Virol. 2012; 86: 2676-2685
        • Mirkovic RR
        • Schmidt NJ
        • Yin-Murphy M
        • Melnick JL.
        Enterovirus etiology of the 1970 Singapore epidemic of acute conjunctivitis.
        Intervirology. 1974; 4: 119-127
        • Nidaira M
        • Kuba Y
        • Saitoh M
        • Taira K
        • Maeshiro N
        • Mahoe Y
        • et al.
        Molecular evolution of VP3, VP1, 3C(pro) and 3D(pol) coding regions in Coxsackievirus group A type 24 variant isolates from acute hemorrhagic conjunctivitis in 2011 in Okinawa.
        Japan. Microbiol Immunol. 2014; 58: 227-238
        • Nielsen R
        • Yang Z.
        Likelihood models for detecting positively selected amino acid sites and applications to the HIV-1 envelope gene.
        Genetics. 1998; 148: 929-936
        • Oberste MS
        • Schnurr D
        • Maher K
        • Al-Busaidy S
        • Pallansch MA.
        Molecular identification of new picornaviruses and characterization of a proposed enterovirus 73 serotype.
        J Gen Virol. 2001; 82: 409-416
        • Oh MD
        • Park S
        • Choi Y
        • Kim H
        • Lee K
        • Park W
        • et al.
        Acute hemorrhagic conjunctivitis caused by coxsackievirus A24 variant, South Korea, 2002.
        Emerg Infect Dis. 2003; 9: 1010-1012
        • Palacios G
        • Oberste MS.
        Enteroviruses as agents of emerging infectious diseases.
        J Neurovirol. 2005; 11: 424-433
        • Posada D.
        JModelTest: phylogenetic model averaging.
        Mol Biol Evol. 2008; 25: 1253-1256
        • Rambaut A
        • Lam TT
        • Max Carvalho L
        • Pybus OG
        Exploring the temporal structure of heterochronous sequences using TempEst (formerly Path-O-Gen).
        Virus Evol. 2016; 2 (vew007)
        • Smura T
        • Blomqvist S
        • Vuorinen T
        • Ivanova O
        • Samoilovich E
        • Al-Hello H
        • et al.
        The evolution of Vp1 gene in enterovirus C species sub-group that contains types CVA-21, CVA-24, EV-C95, EV-C96 and EV-C99.
        PLOS ONE. 2014; 9: e93737
        • Sousa Jr, IP
        • Burlandy FM
        • Ferreira JL
        • Alves JCS
        • Sousa-Júnior EC
        • Tavares FN
        • et al.
        Re-emergence of a coxsackievirus A24 variant causing acute hemorrhagic conjunctivitis in Brazil from 2017 to 2018.
        Arch Virol. 2019; 164: 1181-1185
        • Suchard MA
        • Lemey P
        • Baele G
        • Ayres DL
        • Drummond AJ
        • Rambaut A.
        Bayesian phylogenetic and phylodynamic data integration using BEAST 1.10.
        Virus Evol. 2018; 4 (vey016)
        • Tamura K
        • Stecher G
        • Peterson D
        • Filipski A
        • Kumar S.
        MEGA6: molecular evolutionary genetics analysis version 6.0.
        Mol Biol Evol. 2013; 30: 2725-2729
        • Tao Z
        • Wang H
        • Li Y
        • Xu A
        • Zhang Y
        • Song L
        • et al.
        Cocirculation of two transmission lineages of echovirus 6 in jinan, china, as revealed by environmental surveillance and sequence analysis.
        Appl Environ Microbiol. 2011; 77: 3786-3792
        • Tao Z
        • Wang H
        • Liu Y
        • Li Y
        • Jiang P
        • Liu G
        • et al.
        Non-polio enteroviruses from acute flaccid paralysis surveillance in Shandong Province, China, 1988–2013.
        Sci Rep. 2014; 4: 6167
        • Tavares FN
        • Campos Rde M
        • Burlandy FM
        • Fontella R
        • de Melo MM
        • da Costa EV
        • et al.
        Molecular characterization and phylogenetic study of coxsackievirus A24v causing outbreaks of acute hemorrhagic conjunctivitis (AHC) in Brazil.
        PLoS One. 2011; 6: e23206
        • World Health Organization
        Isolation and identification of polioviruses, WHO polio laboratory manual.
        4th ed. World Health Organization, Geneva2004
        • Wu B
        • Qi X
        • Xu K
        • Ji H
        • Zhu Y
        • Tang F
        • et al.
        Genetic characteristics of the coxsackievirus A24 variant causing outbreaks of acute hemorrhagic conjunctivitis in Jiangsu, China, 2010.
        PLoS One. 2014; 9: e86883
        • Xia X.
        DAMBE7: New and improved tools for data analysis in molecular biology and evolution.
        Mol Biol Evol. 2018; 35: 1550-1552
        • Yang J
        • Lin Y
        • Wang HY
        • Tao ZX
        • Li Y
        • Chen P
        • et al.
        Identification and genetic characterization of coxsackievirus A24 isolated from patients with acute hemorrhagic conjunctivitis in Shandong Province.
        Bing Du Xue Bao. 2012; 28: 663-669
        • Yang Z
        • Nielsen R
        • Goldman N
        • Pedersen AM.
        Codon-substitution models for heterogeneous selection pressure at amino acid sites.
        Genetics. 2000; 155: 431-449
        • Yen YC
        • Chu PH
        • Lu PL
        • Lin YC
        • Shi YY
        • Chou LC
        • et al.
        Phylodynamic characterization of an ocular-tropism coxsackievirus A24 variant.
        PLoS One. 2016; 11e0160672
        • Yin-Murphy M
        • Lim KH
        • Ho YM.
        A Coxsackievirus type A24 epidemic of acute conjunctivitis.
        Southeast Asian J Trop Med Public Health. 1976; 1: 1-5
        • Zhang L
        • Zhao N
        • Huang X
        • Jin X
        • Geng X
        • Chan TC
        • et al.
        Molecular epidemiology of acute hemorrhagic conjunctivitis caused by coxsackie A type 24 variant in China, 2004–2014.
        Sci Rep. 2017; 7: 45202
        • Zocher G
        • Mistry N
        • Frank M
        • Hähnlein-Schick I
        • Ekström JO
        • Arnberg N
        • et al.
        A sialic acid binding site in a human picornavirus.
        PLoS Pathog. 2014; 10e1004401