Molecular Epidemiology And Genetic Characteristics Of Coxsackievirus A16 Circulating In Children In Beijing,China,2010 To 2019

Coxsackievirus A16（CVA16）is a member of species Enterovirus A in genus Enterovirus（EV）,family Picornaviridae.The prototype strain G-10 was first isolated in South Africa in 1951.Up to now,the global CVA16 strains can be divided into three genotypes,A,B and D,based on the phylogenetic tree and genetic diversity of VP1gene.Genotype B can be further divided into three sub-genotypes,B1,B2 and B3.Sub-genotype B2 was the predominant type before 2000,then replaced by B1 strains.Sub-genotype B1 contains three clades,B1a,B1b and B1c.CVA16,EV-A71 and CVA6 are the main pathogens of hand,foot,and mouth disease（HFMD）,and mainly infected children under 5 years old.In recent decades,CVA16 shows a high prevalence in Asian-Pacific region,and was responsible for several HFMD outbreaks.HFMD caused by CVA16 infection is generally mild and self-limiting compared with that by EV-A71.However,CVA16 can occasionally cause severe and fatal cases.EV contains a non-segmented positive-strand RNA genome,coding a RNA-dependent RNA polymerase（Rd Rp）which lacks a proof-reading mechanism.Recombination and mutation caused by the viral error-prone polymerase are the main mechanism for EV evolution.Mutations in the key sites and recombination can contribute to the appearance of new lineages with altered properties,including the change of antigenicity and pathogenicity,which may result in immune escape and resistance to antiviral,and increase the transmission of EV.Objective:To investigate the molecular epidemiological features and evolution characteristics of CVA16 circulating in children with HFMD and suspected HFMD in Beijing from 2010 to 2019.To identify the predominant genotype in Beijing and to explore the genetic diversity,evolution rate,amino acid mutation and recombination of the local CVA16 strains.Methods:The screening and typing results of EV were analyzed to investigate retrospectively the epidemic features of CVA16.Patients involved in this study visited the Infectious Department of Children’s Hospital of Capital Institute of Pediatrics during the period from March,2010 to October,2019.Throat swabs were collected from patients with clinical diagnoses of HFMD,herpangina and rash and fever illness.CVA16 positive samples were selected randomly and proportionally from each year to amplify the complete VP1 gene by RT-PCR with specific primers,and the amplified products were sequenced by Sanger sequencing.Sequences were edited using DNAStar v.5.01 software.To investigate the molecular epidemiology of CVA16 and the diversity of VP1 gene,homology and phylogenetic analyses were conducted by MEGA v.6.05 software.BEAST v.1.10.4 software were used to reconstructed maximum clade credibility tree to confirm the genotype and molecular evolution features of VP1 gene of Beijing CVA16 strains.The evolution rate,the most recent common ancestor（TMRCA）and Bayesian skyline plot of VP1 gene were estimated.The natural selected sites were inferred by Datamonkey.To understand the genetic characteristics of CVA16 more completely,CVA16 strains were selected randomly and proportionally from each year to acquire the whole-genome sequences by RT-PCR and Sanger sequencing.The genotype of CVA16 strains was analyzed by phylogenetic analysis.The diversity of nucleotide and amino acid sequences of CVA16genomes and the mutation of neutralizing epitopes of capsid protein were also analyzed.To investigate the evolution characteristics of coding genes,the evolution rate,TMRCA and natural selected sites were estimated.Recombination was detected by RDP 4 and Sim Plot v.3.5.1 software and verified by the phylogenetic tree.The recombination forms of Beijing CVA16 strains were identified by phylogenetic trees based on VP1and 3D genes.Results:1)In the retrospective analysis,a total of 4 709 throat swabs were screened.EV were detected in 3 180 specimens,and 814 were positive for CVA16.Among CVA16-positive children,81.6%were younger than 5 years,and most of them were in the age of 1-3 years old.From 2010 to 2019,the prevalence of CVA16 in Beijing showed obvious periodic fluctuations.CVA6 replaced CVA16 and EV-A71 as the dominant pathogen in 2013.Although CVA16 and EV-A71 became the main pathogens again in 2014,CVA6 predominated again from 2015 to 2019.However,the prevalence and effective population size of CVA16 increased steadily during 2018～2019.According to the 172 complete VP1 sequences acquired,the estimated mean evolution rate was 4.49×10^-3 substitution/site/year(95%HPD:3.98×10^-3～4.98×10^-3).All the Beijing CVA16 strains belonged to sub-genotype B1,and B1b was the most predominant strain circulating continuously during the period of 2010-2019.One B1c strain was detected in Beijing for the first time in 2016.The Methionine fixed at site-23 of VP1 since 2012,which may be related to the persistent spread of clade B1b.VP1gene of Beijing CVA16 strains underwent strong purifying selection.Only two sites were detected under episodic positive selection,one of which（site-223）located in neutralizing linear epitope PEP71.2)In total,81 whole-genome sequences of CVA16strains were obtained.All of the strains belonged to sub-genotype B1,including 8 B1a and 73 B1b strains.The similarity of nucleotide and amino acid sequences of P1 region were higher than those of P2 and P3.As for P1,the similarities of nucleotide and amino acid sequences of VP1,VP2 and VP3 genes were higher than that of VP4.Different encoding gene had a different evolution rate and shared similar emerging time of TMRCA.The diversity over time of different gene presented a different trend,but the trend of VP1 gene was almost consistent with that of the whole genome.Sporadic amino acid mutations occurred in the neutralizing epitopes of capsid protein.Most genes were under purifying selection.Many mutations located at the positive selected sites of 2A and 3D protein,which might be associated with the prevalent of CVA16.Two recombination strains were detected by bioinformatics analysis.One recombination occurred between clade B1a and B1b of CVA16,the other occurred between CVA16 and CVA2.The recombination breakpoints located in 3B and 5’UTR gene,respectively.Beijing CVA16 strains involved into two recombination forms,RF-B and RF-D,and the most strains belonged to RF-D.Conclusions:From 2010 to 2019,the circulating Beijing CVA16 strains belonged to B1 genotype,clade B1b was the predominant strain.Variation in the genome of Beijing CVA16 strains was extensive,and mutation was the significant evolutionary mode.An increased epidemic trend and effective population size of CVA16 was observed from 2018 to 2019,which might be caused by the mutations in key sites and intertypic recombination event.Continuous surveillance of the molecular epidemiology and genetic background of CVA16 is necessary.