Genomic epidemiology of coxsackievirus A16 in mainland of China, 2000-18

Abstract

Hand, foot, and mouth disease (HFMD), which is a frequently reported and concerning disease worldwide, is a severe burden on societies globally, especially in the countries of East and Southeast Asia. Coxsackievirus A16 (CV-A16) is one of the most important causes of HFMD and a severe threat to human health, especially in children under 5 years of age. To investigate the epidemiological characteristics, spread dynamics, recombinant forms (RFs), and other features of CV-A16, we leveraged the continuous surveillance data of CV-A16-related HFMD cases collected over an 18-year period. With the advent of the EV-A71 vaccine since 2016, which targeted the EV-A71-related HFMD cases, EV-A71-related HFMD cases decreased dramatically, whereas the CV-A16-related HFMD cases showed an upward trend from 2017 to October 2019. The CV-A16 strains observed in this study were genetically related and widely distributed in the mainland of China. Our results show that three clusters (B1a-B1c) existed in the mainland of China and that the cluster of B1b dominates the diffusion of CV-A16 in China. We found that eastern China played a decisive role in seeding the diffusion of CV-A16 in China, with a more complex and variant transmission trend. Although EV-A71 vaccine was launched in China in 2016, it did not affect the genetic diversity of CV-A16, and its genetic diversity did not decline, which confirmed the epidemiological surveillance trend of CV-A16. Two discontinuous clusters (2000-13 and 2014-18) were observed in the full-length genome and arranged along the time gradient, which revealed the reason why the relative genetic diversity of CV-A16 increased and experienced more complex fluctuation model after 2014. In addition, the switch from RFs B (RF-B) and RF-C co-circulation to RF-D contributes to the prevalence of B1b cluster in China after 2008. The correlation between genotype and RFs partially explained the current prevalence of B1b. This study provides unprecedented full-length genomic sequences of CV-A16 in China, with a wider geographic distribution and a long-term time scale. The study presents valuable information about CV-A16, aimed at developing effective control strategies, as well as a call for a more robust surveillance system, especially in the Asia-Pacific region.

Keywords: coxsackievirus A16; enterovirus; epidemiology; phylodynamics; phylogeny.

Figure 1.

The numbers of probable cases, laboratory-confirmed cases, severe cases, and two major pathogens reported. The dashed line chart shows the sentinel surveillance data of probable HFMD cases reported from 2009 to October 2019, corresponding to the right vertical coordinates. The solid line chart represents the laboratory-confirmed cases, severe cases, EV-A71, and CV-A16, respectively, corresponding to the left vertical coordinates.

Figure 2.

The genetic diversity and statistic numbers of CV-A16 in this study. (A) The neighbor-joining phylogenetic tree based on the entire VP1 genome of seventy-one representative strains for genotyping. The numbers at each node show the bootstrap values constructed using the neighbor-joining method with 1,000 bootstrap replicates. (B) The numbers and geographic distributions of the CV-A16 sub-genogroups in this study (n = 166). (C) The bar chart showing the CV-A16 sample distributions of all representative full-length genomes across five geographic areas in China (n = 271). The corresponding genotypes, which were detected in China, are labeled in the bar plot. (D) The genetic diversity of the open reading frame (ORF) genome of CV-A16 in this study, which was calculated using a sliding window of 300 nucleotides with a step size of twenty-five nucleotides. The red lines show the partitive genomic region of P1, P2, and P3, respectively. The colored lines represent the corresponding areas, as shown in Fig. 2C.

Figure 3.

The temporal phylogenies and epidemic characteristics of CV-A16 estimated with the P1 genomic sequences. (A) The relative genetic diversity of the CV-A16 sequences in China. The x-axis represents the units of year, and the y-axis shows the measure of genetic diversity (logarithmic scale of Neτ, where Ne is the effective population size and τ is the generation time). The black line shows the median estimates of the CV-A16 population size and the green shading represents the 95 per cent CI. (B) The histogram of the average number of state transitions based on five geographic locations. (C) The MCC phylogenetic tree based on the entire P1 coding region in China and colored according to the different areas.

Figure 4.

Spatial diffusion of CV-A16 in China. (A) Spatial transmission pathways of CV-A16 inferred using the Bayesian Stochastic Search Variable Selection method. The solid black arrow shows decisive support for the diffusion pathway identified by the P1 and VP1 coding regions, with BF > 1,000 and indicator > 0.5. The dashed black arrow represents strong support for the diffusion pathway identified by one of the P1 and VP1 coding regions, with BF > 1,000 and indicator > 0.5. (B) The median inferred immigration events identified by the Marginal Structured Coalescent Approximation (MASCOT) method. The width of the arrows indicates the possibility of transmission. The dot sizes are proportional to the median inferred effective population sizes. The five regions of China are defined and colored.

Figure 5.

(A) and (C) The scatterplots show the first two principal components of the DAPC of the CV-A16 sequences in mainland China using years of sampling and locations of samples as prior clusters. Eigenvalues of the analysis are displayed in the inset. Groups are shown by different colors and dots representing individual isolates. (B) and (D) Assignment probabilities of CV-A16 individuals from the five sampling locations and years. The panels of CV-A16 individuals are shown according to the years of sampling and locations of samples. The scatterplot was analyzed using the CV-A16 ORF coding region.

Figure 6.

Similarity analysis of the CV-A16 strains with potential parents. The dataset of seventy-two ORF genomic sequences, including the representative genomic sequences available from GenBank and this study, was used. The red lines show the partitive genomic region of P1, P2, and P3, respectively. (A) Group 1, including three high identity sequences, was used as a query sequence. (B) Group 2 was used as a query sequence. (C) Group 3 was used as a query sequence.