Coxsackievirus B4: a Undervalued Pathogen that Associated with Hand, Foot and Mouth Disease Outbreak

Objectives: This study was conducted to discover the causes of Coxsackievirus B4 (CVB4) -induced hand, foot, and mouth disease (HFMD) outbreaks and its evolutionary characteristics. Methods: In this study, we sequenced isolates obtained during the outbreak for comparative analyses with previously sequenced strains. Phylogenetic analysis and evolutionary dynamics were performed to illustrate the genetic characteristics of CVB4 in China and worldwide. Results: The nucleotide sequence of CVB4 isolated during the outbreak in 2011 was more similar to that of CVB4 isolated in Shandong Province, China in 2010 (95.7–99.4%) than to other CVB4 isolated in China (90.9–98.8%). A phylogenetic analysis showed that CVB4 originated from a common ancestor in Shandong. CVB4 strains isolated worldwide could be classied into genotypes A–E according to the VP1 region. All CVB4 strains in China belonged to genotype E. The global population diversity of CVB4 uctuated substantially over time, and CVB4 isolated from China accounted for a signicant increase in diversity of CVB4. The average nucleotide substitution rate in VP1 of Chinese CVB4 (5.20 × 10 -3 substitutions/site/year) was slightly higher than that of global CVB4 (4.82 × 10 -3 substitutions/site/year). Conclusions: These ndings explain both the 2011 outbreak and a global increase in CVB4 diversity. In addition to improving our understanding of a major outbreak, these ndings provide a basis for the development of surveillance strategies.


Introduction
Coxsackieviruses belong to genus Enterovirus within the family Picornaviridae, which are small non-enveloped RNA viruses with a single-stranded positive sense genome of approximately 7500 nucleotides (nt). They are associated with various diseases, such as paralytic diseases, hand, foot, and mouth disease (HFMD), meningitis-encephalitis, myocarditis, febrile diseases, and respiratory infections [1][2][3].
According to their pathogenicity in suckling mice, coxsackieviruses can be divided into two groups: A and B [4].
CVB has caused many disease outbreaks worldwide [13][14][15]. Coxsackievirus B3 (CVB3) and Coxsackievirus B5 (CVB5) are considered the main CVB serotypes causing HFMD and aseptic meningitis [16][17][18]. Surveillance data in China have shown that the incidence of HFMD in Shandong Province is high, including an outbreak in 2011. Coxsackievirus B4 (CVB4) infections were found in 21 patients with HFMD. There are few reports of HFMD outbreaks caused by CVB4 in China and elsewhere, which makes the study meaningful.
CVB3 is a frequent cause of disease outbreaks; accordingly, studies of the evolution and clinical characteristics of CVB4 are needed [19].
In this study, the entire VP1 region of 35 CVB4 strains (21 strains isolated from patients and 14 strains from healthy individuals) were obtained. Fourteen CVB4 strains were isolated from HFMD surveillance specimens in Shandong Province in 2010, and 21 CVB4 strains were isolated from HFMD outbreak specimens in Shandong Province in 2011. All available CVB4 genome sequences were analyzed by a bioinformatics approach and analyzed with respect to clinical characteristics. Combined with HFMD surveillance results obtained before the outbreak, an epidemiological study of the HFMD outbreak caused by CVB4 in Shandong Province was performed.

Ethical considerations
This study did not involve any human experimentation. The main purpose of the study was to contribute to public health and the prevention of HFMD outbreaks caused by enteroviruses. Samples were obtained from the national laboratory network established by the National Laboratory for Poliomyelitis in the National Institute for Viral Disease Control and Prevention (IVDC) of the Chinese Center for Disease Control and Prevention (China CDC). Written informed consent for the use of clinical samples was obtained from the parents of all patients. Parents agreed to the use of specimens for research purposes and to the collection of samples by trained professional staff. This study was approved by the second session of the Ethics Review Committee of IVDC. All experimental protocols were approved by IVDC, and the methods were carried out in accordance with the approved guidelines.

Sample collection and virus isolation
In total, 197 stool samples from healthy children and 178 stool samples from patients with HFMD in an outbreak were collected in Shandong Province in 2010 and 2011, respectively. The diagnostic criteria were children with rashes and herpes on the hands, feet, and mouth accompanied by fever symptoms. Most patients had mild symptoms. A small number of patients with severe infections may develop aseptic meningitis, viral myocarditis, encephalitis, acute accid paralysis, and other symptoms, and further development may lead to cardiopulmonary insu ciency and eventual death [20,21]. The 375 specimens were processed according to the Polio Laboratory Manual (World Health Organization, Geneva, Switzerland). The processed samples were inoculated into human rhabdomyosarcoma provided by the WHO Global Poliovirus Specialty Laboratory in the United States. Continuous monitoring was performed and harvested cell cultures after the cells exhibited a complete EV-like cytopathic effect (CPE).

Molecular typing and VP1 sequencing
Viral RNA was extracted from the harvested cell culture using the QIAamp Viral RNA Mini Kit (Qiagen, Hilden, Germany).
The One-step RT-PCR Kit (Perfect Real Time, TaKaRa, Dalian, China) was used for enterovirus-speci c real-time polymerase chain reaction (real-time PCR) [22]. Enterovirus-positive specimens determined by real-time PCR were subjected to reverse transcription polymerase chain reaction (RT-PCR) using the PrimeScript One Step RT-PCR Kit Ver. 2 (TaKaRa) and primers E486/ E488, E490/E492, and E494/E496 were used to amplify the partial VP1 region [23]. The PCR product was puri ed using a QIAquick PCR Puri cation Kit (Qiagen), followed by sequencing using an ABI 3130 Genetic Analyzer (Applied Biosystems, Foster City, CA, USA). Sequences were analyzed using the EV Genotyping Tool and the BLAST server [24]. To obtain the entire VP1 sequence of CVB4, speci c primers were designed (upstream primer, CVB4-F: GTTATAGTTCCAGCCGAAGCG; downstream primer, CVB4-R: CCACACACAACTCTGCCAATC). The PCR product was puri ed and sequenced.

Bioinformatics Analysis
A total of 209 CVB4 VP1 sequences were obtained from the GenBank database and Muscle (v3.8.31_i86linux32) was used to generate alignments [25,26]. After screening and removing low-quality sequences (incomplete VP1 region or sequencing error), they were combined with 35 entire VP1 sequences of CVB4 strains isolated from Shandong Province in 2010 and 2011 for a phylogenetic analysis. RAxML (v8.2.12) was used to construct maximum likelihood trees based on 215 complete VP1 coding regions of CVB4 [27]. Nucleotide substitution models were GTR+G, as calculated using ModelGenerator 0.85 [28]. Support was estimated from 1000 bootstrap replicates, and the results were visualized using FigTree (v1.4.4). Maximum likelihood trees based on sequence set were constructed using MEGA (v7.0) [29] with 1000 bootstrap replicates to verify the tree topology obtained using RAxML. According to phylogenetic trees, genotypes were de ned using a 15-25% divergence threshold for the VP1 coding region [30,31] [32]. BioEdit (v7.0.9.0) was used to obtain nucleotide and amino acid sequence similarities [33].
To evaluate the global epidemic dynamics of CVB4 over time, 200 entire VP1 sequences of CVB4 containing time and geographic information were analyzed. Sequences were analyzed using BEAST (v1.8.4) [34][35][36]. The nucleotide substitution model was GTR+G, as calculated using ModelGenerator 0.85. The optimal combination of molecular clock model and tree prior was determined. A Bayesian Markov chain Monte Carlo (MCMC) analysis was run for 2 × 10 8 generations to ensure the convergence of each parameter. The sampling frequency was set to 2 × 10 4 generations. Tracer (v1.6) was used to check convergence and optimize the molecular clock and tree prior and to determine the estimated sample sizes (ESS) [37]. The ESS value should be at least 200. To analyze the dynamics of population diversity, a Bayesian skyline plot based on 200 CVB4 entire VP1 sequences was constructed using Tracer [38]. TreeAnnotator (v1.8.4) was used to construct a Bayesian maximum clade credibility (MCC) tree. The rst 10% of sampled trees were removed with the burn-in option. The phylogenetic tree was visualized using FigTree.
In total, 70 VP1 sequences of CVB4 strains isolated from China were selected for an analysis of population dynamics using BEAST. The strict clock model and the Bayesian skyline tree prior were used. A Bayesian MCMC was run for 3 × 10 8 generations to ensure that each parameter can converge. The sampling frequency was set to 3 × 10 4 generations. Tracer was used to evaluate parameters. An MCC tree was constructed using TreeAnnotator. The burn-in option was used to remove the rst 10% of sampled trees. FigTree was used to visualize the phylogenetic tree. Nucleotide and amino acid sequence similarities for comparisons between VP1 sequences of CVB4 strains isolated from Shandong Province in 2011 and other CVB4 VP1 sequences were estimated using BioEdit.  (Table 1). A phylogenetic tree based on all entire VP1 sequences showed that CVB4 can be divided into genotypes A, B, C, D, and E (Fig. 2). Genotype A included the prototype strain (X05690/JVB/New York/US/1951) and the other two CVB4 strains (GenBank accessions S39291 and DQ480420). Genotype B-D contain strains from several countries other than China. Genotype E consisted of CVB4 strains isolated from China (N = 70) and Australia (N = 3). CVB4 strains isolated from China were all classi ed as genotype E, including the new sequences obtained in this study. CVB4 in China formed essentially an independent lineage in the phylogeny. Genotype E was the most prevalent CVB4 genotype in China. Basic data for each genotype are summarized in Table 1. With an uncorrelated relaxed molecular clock model (exponential) and exponential growth tree prior, MCC trees based on sequences of 200 globally distributed CVB4 were obtained using BEAST (Fig. 3). According to a Bayesian skyline plot of the entire CVB4 VP1 sequence (Fig. 3 the population of CVB4 showed a slow decline, followed by a rapid decline. Isolates mainly belonged to genotype D during the rst half of this stage and to genotype E, mainly in China, in the second half of the stage. In the fth stage, beginning in 2011, the diversity of CVB4 increased rapidly and then remained steady. During this period, there were many epidemics and outbreaks of CVB4, which may be related to the rapid increase in diversity. Gradual stabilization since then may be due to a decrease in the quantity of CVB4 found in disease surveillance.

Phylogenetic analysis of CVB4 in China
A total of 70 CVB4 strains isolated from China were analyzed using BEAST with strict molecular clock model and Bayesian skyline tree prior (Fig. 4). The MCC tree showed that CVB4 strains isolated from the HFMD outbreak in Shandong Province in 2011 and a portion of CVB4 strains isolated from healthy individuals in 2010 diverged from a common ancestor (with a posterior probability of up to 0.89). The nucleotide sequence similarity between VP1 sequences of CVB4 in Shandong in 2011 was 99.6-100% and the amino acid sequence similarity was 100%. The nucleotide sequence similarity in comparison with CVB4 strains isolated in Shandong Province in 2010 was 95.7-99.4% and amino acid sequence similarity was 97.8-98.9%. The nucleotide sequence similarity with other CVB4 isolates from China was 90.9-98.8% and amino acid sequence similarity was 97.5-99.2%.

Discussion
Frequent outbreaks of HFMD worldwide, disproportionately affecting children under 5 years of age, are a major public health issue. If patients are not treated in time, infant death may occur [40][41][42][43]. The main causal pathogens are EV-A71, CVA16, CVA6, and CVA10 [44][45][46][47]. CVB3 and CVB5 are related to HFMD and aseptic meningitis [16,18] We further found that the global population diversity experienced three periods of growth. In the rst period , CVB4 gradually increased and was prevalent in Europe, North America, and Oceania. In the second period, 1990-1991, diversity showed a temporary upward trend. The global prevalence of genotype B and D increased, explaining the overall increase in diversity. In the third period, beginning in 2011, a sharp and steep rise in diversity was detected. During this period, there were many outbreaks of CVB4, including the outbreak in Shandong Province evaluated in this study, which may explain the general increase in diversity. More than half a century of evolution has allowed CVB4 to spread and cause an epidemic, dominated by genotypes D and E. Decreases in other genotypes may be attributed to a lack of clinical symptoms or inadequate monitoring. These genotypes have been reported in Europe, North America, Africa, Australia, and other regions, and monitoring should be strengthened, particularly since extant genotypes have the potential ability to cause disease outbreaks.
The average nucleotide substitution rate of the CVB4 VP1 coding region was higher than those of genotypes B and C of EV-A71 VP1 (4.6 × 10 − 3 substitutions/site/year). Additionally, the substitution rate of China CVB4 was slightly higher than that of the global CVB4 VP1 coding region. These results indicate that in China, CVB4 is active and has a high mutation rate; this may also explain the outbreak of HFMD in Shandong Province. Furthermore, current strains of CVB4 in China are genotype E. The strain responsible for the outbreak of CVB4 in Shandong Province is derived from local CVB4 and not from abroad. Some CVB4 strains isolated from healthy children in 2010 were closely related to CVB4 strains isolated from HFMD samples from the outbreak in 2011. They shared a common ancestral sequence and exhibited higher nucleotide sequence similarity with each other than with other CVB4 strains isolated in China at other times. These results indicate that the outbreak of HFMD caused by CVB4 in Shandong Province in 2011 could be attributed to a pathogen that existing in the population as early as 2010. Eventually, an outbreak occurred in Shandong in 2011, proving that CVB4 has the ability to cause an outbreak but usually causes a subclinical or asymptomatic infection in the healthy population.

Conclusions
Our analysis of 35 entire VP1 sequences of CVB4 strains isolated in this study combined with CVB4 data available in GenBank provided insight into the global phylogenetic and epidemic characteristics of CVB4. After more than half a century, co-circulation of genotypes D and E resulted in the reported epidemics and outbreaks [48,39]. The outbreak of HFMD caused by CVB4 in Shandong Province is an example. These results broaden the pathogenic spectrum of HFMD.
In addition, CVB4 isolated from the healthy population in 2010 was closely related to CVB4 isolated from patients with HFMD in 2011, indicating that the strain spread in the healthy population before the outbreak. When rare pathogens exhibit a sudden increase in the healthy population, careful monitoring is needed owing to the risk of HFMD or other diseases.

Declarations
Ethics approval and consent to participate This study was approved by the Ethics Review Committee of the National Institute for Viral Diseases Control and Prevention, Chinese Center for Disease Control and Prevention. Written informed consent for the use of clinical samples was obtained from all individuals, including written informed consent from guardians of children.

Consent to participate
Written informed consent was obtained from the parents.

Consent for publication
Not applicable.

Availability of data and materials
The 35 CVB4 VP1 sequences determined in this study have been deposited in GenBank under accession numbers MT109013-MT109047.

Competing interests
The authors declare that no competing interests exist.