Molecular Characterization of Coxsackievirus B1 Strains Isolated from Patients with Hand Foot and Mouth Disease in Yunnan, Southwest China

Coxsackievirus (CV) B belongs to the species Enterovirus B, genus Enterovirus of the family Picornaviridae. Enterovirus B (EV-B) includes 63 serotypes: CVB1–6; CVA9; echoviruses E1–7, 9, 11–21, 24–27, and 29–33; EV-B69, EV-B 73–75, EV-B 77–88, EV-B 93, EV-B 97–101, EV-B 106–107, and EV-B 110–114 (www.picornaviridae.com/ensavirinae/enterovirus/ev-b/ev-b.htm). Like other EV-Bs, CVBs are single-stranded, small non-enveloped and positive-sense RNA viruses. The CVB genome is ≥ 7,400 nucleotides long and includes a large open reading frame (ORF), encoding a single pre-polyprotein that is further processed into P1, P2, and P3. P1 is cleaved to VP4, VP2, VP3, and VP1 structural proteins and P2, and P3 is cleaved into seven nonstructural proteins (2A, 2B, 2C, 3A, 3B, 3C, and 3D)^[1]. Mature structural proteins can be used for viral assembly, while nonstructural proteins can be used for replication, induction of apoptosis, suppression of innate immunity, and shutting down host cell translation^[2].

CVB1 is commonly associated with mild infections [fever, irritation, and hand, foot and moth disease (HFMD)]; however, some patients may develop serious illness, such as aseptic meningitis, myocarditis, acute flaccid myelitis, and necrotizing hepatitis. CVB1 is also associated with type 1 diabetes^[3]. Since 2008, CVB1 has frequently been detected in HFMD cases and associated with severe disease^[4]. In 2019, CVB1 was one of the most frequent causative pathogens of HFMD in Yunnan, China^[5]. Eleven complete genomes of CVB1 are currently accessible in GenBank, and three strains, CVB1/XM0108 (MG780414), MSH/KM9/2009 (JN596588), and CVB1SD2011CHN (JX976769), were isolated in China between 2009 and 2011. Hence, in the study, the six full-length genome sequences of CVB1 strains were determined and compared with those of other EV-Bs.

During the HFMD surveillance in 2019, 242 stool samples were collected from children with HFMD under 5 years old in Yunnan, Southwest China. The green monkey kidney (Vero), the human rhabdomyosarcoma (RD) and human embryonic lung diploid fibroblast (KMB17) cell lines were used to isolate the samples, respectively. Viral RNA was extracted from positive culture supernatants that produced typical EV-like cytopathic effect (CPE) using the QIAamp Viral RNA Mini Kit (Qiagen, CA, USA). Reverse transcription polymerase chain reaction (RT-PCR) was performed using the PrimeScript™ One Step RT-PCR Kit Ver.2 (TaKaRa, Dalian, China). The primers BOAS and BOS were used, and the positive RT-PCR products were sequenced on an ABI 3130 Genetic Analyzer (Applied Biosystems, USA) at Kunming Qingke Biotechnology Co. (Yunnan, China). Six CVB1 strains (7V3/YN/CHN/2019, 46V3/YN/CHN/2019, 53V3/YN/CHN/2019, 55V3/YN/CHN/2019, 62V3/YN/CHN/2019, and 64V3/YN/CHN/2019) were identified as CVB1 by enterovirus typing online software (https://www.rivm.nl/mpf/typingtool/enterovirus/). However, the CVB1 isolates were only isolated from Vero cells, whereas none were from RD and KMB17 cells. This indicated CVB1 is sensitive to Vero cells. Then, two long-fragment PCR amplifications of viral RNA were performed using PrimeScript^TM One Step RT-PCR Kit Ver.2 (TaKaRa) and the following primers: E201F (5’ TTAAAACAGCCTGTGGGGTTG 3’) and B12R (5’ TTAGTGGGCCGTGAACTAAT 3’) and B13F (5’ AACCCAAGCACGTGAAGGCG 3’) and EV8R (5’ CACCGAATGCGGAGAATTTA 3’). The sequencing primers for the complete genome were designed according to the primer walking strategy^[6]. All primers used in this study are shown in Supplementary Table S1 (available in www.besjournal.com). Nucleotide and amino acid sequence comparisons were carried out using Geneious Basic 5.4.1 Beta software. Phylogenetic analyses were performed using MEGA 7.037. Recombination Detection Program (RDP) 4.101 software in default mode was used to determine recombinant sequences using seven algorithms (RDP, GENECONV, BootScan, Maxchi, Chimaera, SiScan, and 3Seq). The complete genome sequences of the six CVB1 strains isolated in this study were submitted to GenBank under accession numbers OQ259520–OQ259525.

The whole VP1 sequences (834 nucleotides) of the six isolates showed the highest identity (97.7%–97.8%) with CVB1 strain CV-B1/SWG69/SD/CHN/2018 (MN541037), isolated from sewage in Shandong, China, 2018. The six isolates shared 79.3%–80.2% nucleotide identity and 93.9%–94.2% amino acid identity with the whole VP1 sequence of the CVB1 prototype strain Conn-5 (M16560) and 80.1%–97.7% nucleotide and 93.9%–100% amino acid identity with other CVB1 strains. The complete VP1 nucleotide and amino acid identities among the six isolates were 94.4%–100% and 98.2%–100%, respectively.

The 131 whole VP1 sequences available in GenBank were contained in this study (Supplementary Figure S1, available in www.besjournal.com). The VP1 sequence divergences between EV genotypes was 15% at the nucleotide level^[7]. All CVB1 strains were divided into genotypes A–E. D and E genotypes were further divided into D1 and D2, and E1 and E2 subgenotypes, respectively (Supplementary Table S2, available in www.besjournal.com). The Chinese CVB1 strains belonged to D1 (2003–2018), E1 (1989–1994), and E2 (1999–2019) subgenotypes. Along with most Chinese strains, the six isolates in the study belonged to E2 subgenotype. In addition, the average differences between the VP1 nucleotide and amino acid sequences of these six isolates and all full VP1 gene sequences of other Chinese strains in GenBank were 10.65% (2.20%–19.10%) and 3.25% (0.4%–6.1%). Compared with the three Chinese CVB1 strains in GenBank, CVB1/XM0108 (MG780414), MSH/KM9/2009 (JN596588), and CVB1SD2011CHN (JX976769), the differences between their whole genome nucleotide and amino acid sequences were 8.90% (4.70%–13.10%) and 6.15% (1.50%–10.80%), respectively. This indicates that CVB1 strains have evolved in China.

Figure S1.  Phylogenetic tree based on complete VP1 sequences (834 bp) of all CVB1 isolates available in GenBank. λStrains isolated in this study. ▲Prototype strain.

The whole genomes of the six CVB1 strains (7V3/YN/CHN/2019, 46V3/YN/CHN/2019, 53V3/YN/CHN/2019, 55V3/YN/CHN/2019, 62V3/YN/CHN/2019 and 64V3/YN/CHN/2019) isolated in Yunnan, Southwest China, 2019 were obtained. The length of their genomes was 7,402–7,412 nucleotides, containing a single ORF of 6,549 nucleotides, and encoding a polyprotein of 2,182 amino acids. The ORF was flanked by a 5′-untranslated region (UTR) of 750–757 nucleotides and 3′-UTR of 105–113 nucleotides. The whole genome nucleotide and amino acid identities of the six strains were 99.6%–99.8% and 99.1%–99.8%, respectively. The overall base compositions were 28.1%–28.2% A, 22.9%–23.0% G, 22.6%–24.7% C, and 24.2%–24.3% U. The whole genome nucleotide sequences and the deduced amino acid sequences of the six strains had > 99.6% mutual identity; therefore, the 46V3/YN/CHN/2019 strain was selected as a representative strain for further analysis. The 46V3/YN/CHN/20193 strain showed 79.0% and 79.0%–95.2% whole genome nucleotide identities and 95.0% and 95.0%–98.8% amino acid sequence identities with the prototype strain Conn-5 and other CVB1 strains, respectively (Table 1).

Phylogenetic analyses were performed on the P1, P2, and P3 regions of all available CVB1 strains and EV-B prototype strains in GenBank, as well as the six Yunnan CVB1 isolates and six EV-B epidemic strains (CVA9/VEN/2018-23122/MK652143, E6/P735/2013/China/KP289439, CVB4/B405/MW015040, CVB4/GX/10/JX308222, CVB4/B404/MW015039, and CVB5/B504/MW015048) (Figure 1). In the P1 coding region, the six isolates clustered with all CVB1 strains including the CVB1 prototype Conn-5 strain, which confirmed the initial VP1 typing result. The CVB1 strains clustered with different EV-B strains to form different clusters in the P2 and P3 coding regions (Figure 1). In the P2 coding region, the six isolates clustered with three CVB1 isolates (CVB1SD2011CHN/JX976769, XM0108/MG780414, and MSH/KM9/2009/JN596588) and CVA9/MK652143/VEN/2018-23122 strain. CVB1 strain RO-98-1-74/LS451286 clustered with E11 strain ROU-9191/AJ577594. CVB1 Nm/EU147493, Conn-5/M16560, and CVB1 strain NC 001472 clustered with E6 prototype lytic/U16283. E3 prototype Morrisey/AY302553 clustered with CVA9 prototype strains Griggs/D00627, pMP1.24/AY186745, pMP1.23/AY186746, pAMP2.17/AY186747, and pAMP2.22/AY186748. The five CVB1 strains, CVB1-Chi07/KJ849619, Groningen/2011/LN854562, 1167438 pmMC/JN797615, B1L071615/MG845887, and 2104/08/KP260537, exhibited clustering. In the P3 region, the six isolates only clustered with one CVB1 strain, namely CVB1 SD2011CHN/JX976769, along with five EV-B strains (E6/P735/2013/China/KP289439, CVB5 /B504/MW015048, CVB4/GX/10/JX308222, CVB4 /B404 /MW015039, and CVB4//B405 /MW015040). The prototype Conn-5 strain formed a single lineage with the E4 prototype strain Pesacek/AY302557. The findings suggest that multiple instances of recombination have taken place among the CVB1 strains (including this study) and other EV-B strains in the nonstructural coding regions. Furthermore, it is plausible that these CVB1 isolates represent distinct recombinant strains. The online BLAST tool was utilized to screen for sequences that exhibited the top three identities with strain 46V3/YN/CHN/2019 (Supplementary Table S3, available in www.besjournal.com). The 3A to 3D regions of strain 46V3/YN/CHN/2019 shared the highest identity with that of CVB1SD2011CHN, and followed with other EVBs (CVB5 strain B504 and E6 strain E6/E6SD11CHN). In the 3′-UTR, strain 46V3/YN/CHN/2019 shared the greatest identity (93.58%) with E6 strain E6/P735/2013/China.

Figure 1.  Phylogenetic trees based on the P1, P2, and P3 regions of all CVB1 strains and EV-B prototype strains available in GenBank, the six Yunnan isolates, and five EV-B strains: one CVB5 strain B504, three CVB4 strains (GX/10, B405, and B404), and one E6 strain E6/P735/2013/China). Nucleotide sequences were analyzed using the neighbor-joining algorithm (1,000 bootstrap repeats) in MEGA 7.0.26. The scale bars indicate the genetic distance. Numbers at the nodes indicate bootstrap support. Only high boot values (> 75%) are displayed. νStrains isolated in this study. λOther CVB1 strains. ▲CVB1 prototype strain.

Bayesian Markov Chain Monte Carlo (MCMC) analysis of the worldwide CVB1 VP1 sequences was carried out with the combination of the strict clock model, GTR + Gamma nucleotide substitution model and a Bayesian skyline tree prior. The MCMC chain length was set to 200 million states, and sample was taken every 20000 states. After produced through TreeAnnotator v1.10.4, the Maximum Clade Credibility (MCC) tree was visualized by FigTree v1.4.2 (Figure 2). The clustering structure of the MCC tree was highly similar to that of the phylogenetic analysis (Supplementary Figure S1). The mean base substitution rate of CVB1 VP1 was estimated to be 2.0 × 10^–3 substitutions/site/year, with the 95% highest posterior density (HPD) was 1.5 × 10^–3–2.5 × 10^–3. The most recent common ancestor (t_MRCA) of CVB1 worldwide could date back to 1810.1 (95% HPD: 1755.2–1852.9).

Figure 2.  Maximum clade credibility (MCC) tree based on the complete VP1 sequences of global CVB1 strains available in GenBank. The different colours of the tip labels corresponded to the countries where the strains were collected, and the colour panel is placed at the top left part. The time scale is given at the bottom of Figure 2.

Phylogeographic analysis was carried out using BEAST v1.10.4, employing Bayesian stochastic search variable selection (BSSVS) and the Asymmetric Substitution Model. The migration pathways and their Bayes factors (BF) were assessed utilizing SpreaD3 v0.9.6, and subsequently visualized on a map. (Supplementary Figure S2, available inwww.besjournal.com). The results showed that there were 17 major CVB1 migration routes (BF > 5) worldwide from 1948 to 2019. Notably, transmission originating from the United States, either directly or via Africa, to Europe, as well as inter-country transmission within Europe, predominates. Additionally, bidirectional transmission of CVB1 was observed between the United States and Australia, as well as between Mainland China and Europe.

Figure S2.  The spatial transmission of CVB1 worldwide. The green arrows indicate transmission paths with BF > 5.

RDP4 analysis revealed that this isolate was involved in a recombination event with CVB4/GX/10/JX308222. The recombination event occurred between 5,000 and 7,500 nucleotides (Supplementary Figure S3, available in www.besjournal.com). Recombination is one of the major mechanisms of EV-B evolution^[8]. Sequence analysis revealed that the P2 and P3 regions of the 46V3/YN/CHN/2019 strain had the highest homology with the CVB1SD2011CHN strain, which was isolated from a patient with HFMD in Shandong, China in 2010. The CVB1SD11 CHN shares sequences with E6SD11 CHN, CVA9-Alberta-2010, COXB5/Henan/2010, CVB4/GX/10, and E1-HN serotypes in the P2 and P3 regions, which indicates that multiple recombination events have occurred in the nonstructural region of CVB1SD2011CHN. In this study, we found several putative recombination events in nonstructural regions between strain CVB1 46V3/YN/CHN/2019 and CVB5 B504/MW015048, CVB4 strain (GX/10, B405, and B404) and EV-B E6/P735/2013/China. Based on the substantial sequence identity (≥ 90%), it can be inferred that these viruses serve as potential recombination parents. Among them, E6/P735/2013/China and CVB4 strain GX/10 were isolated from patients with HFMD in 2013 and 2010 in China, respectively, CVB4 strains B405 and B404 were isolated from raw sewage in Japan in 2010. CVB infections are often associated with HFMD^[4]. The numbers of HFMD cases infected with CVB4 and CVB5 are increasing in China. These viruses also often co-circulate to cause epidemics or outbreaks of HFMD in China. Yunnan, in southwestern China, is a major tourist province, with hundreds of people from all over the world (including other provinces of China) visiting every year. These characteristics may provide ample opportunities for recombination of these viruses. However, the recombination can make new variants with altered pathogenicity, and new variants have serotype-specific structural protein sequences but share nonstructural protein sequences from different sources (different viruses or identical viruses from different sources)^[9]. In addition, frequent recombination events have been found between CVB4 and EV-B multiple serotypes (CVB5, E9, CVA9, and E20) in the 3Dpol region. This recombination resulted in differences in the 3Dpol, which affected the fidelity of replication and ultimately the virus's environmental adaptability and survival^[10]. Thus, the new variant might cause outbreaks and epidemics. Attention should be paid to cases caused by CVB1 and its derivatives.

Figure S3.  RDP4 analysis of strain 46V3/YN/CHN/2019 with closely related strains. Recombination analysis was performed by seven algorithms: RDP, GENECONV, BootScan, Maxchi, Chimaera, SiScan, and 3Seq.

In conclusion, although the six CVB1 isolates in this study were not isolated from samples collected from throat swabs, blood, and other lesions of HFMD patients, the current prevalence of CVB1 suggests that CVB1 is closely related to HFMD. Recombination and mutation have driven the evolution of CVB1. The six CVB1 strains in this study were isolated only from patients with HFMD who did not present with aseptic meningitis, myocarditis, and acute flaccid myelitis, but systematic epidemiological surveillance is needed to evaluate the association between CVB1 and these diseases. The findings contribute to the expansion of the complete genome of CVB1 in the GenBank database. Nevertheless, it would be beneficial to elaborate on the potential implications of these findings for the understanding of CVB1 epidemiology and the development of preventive measures.