Genetic characteristics of coxsackievirus A16 associated with hand, foot, and mouth disease in Nanjing, China

Introduction: Coxsackievirus A16 (CVA16) is a main pathogen in hand, foot, and mouth disease (HFMD) worldwide. This study intended to clarify the genetic characteristics of CVA16 associated with HFMD in a defined area in Nanjing, China. Methodology: A total of 175 CVA16 strains isolated from throat swabs between 2011 and 2013 were obtained through sentinel hospitals in Nanjing. Multiplex polymerase chain reaction (PCR) was used to amplify the VP1 sequence of local CVA16 strains, and their genetic relationship with 138 CVA16 strains isolated in China and other countries of the world was compared. Results: Phylogenetic analysis based on complete VP1 sequences revealed that subgenotype B1a and B1b were predominantly circulating in Nanjing and B1b strains were spread more widely. The evolution of CVA16 strains is very conservative, with a mean distance of less than 9%. Moreover, six reported conservative regions in VP1 protein were examined, and three of them exhibited high conservation in all CVA16 genotypes except the G-10 prototype and may serve for further vaccine research. Conclusions: The CVA16 strains circulating in Nanjing, China, in 2011 to 2013 belonged to different genotypes and evolved in a conservative way. To provide further evidence for epidemiological linkage and evolutionary recombination events in CVA16, persistent surveillance of HFMD-associated pathogens is required.


Introduction
Coxsackievirus A16 (CVA16) is the first identified human enterovirus (EV) as a causative agent of hand, foot, and mouth disease (HFMD) [1]. HFMD is a common and highly contagious infection that affects infants and children and is caused by Enterovirus of the Picornaviridae family. It is typically characterized by a mild fever followed by small red spots and bumps that may blister on the skin of palms, soles, in the oral cavity, and occasionally on the buttocks or/and genitalia [2]. For decades, multiple outbreaks of HFMD have occurred in many different countries, and CVA16 has played a critical role in the epidemic [2][3][4][5]. In recent years, the morbidity of HFMD has increased significantly in Asian countries including China, Japan, Malaysia, Thailand, and Korea, causing a huge burden for human health and economy [3,4,6].
Since the first outbreak of HFMD in Anhui Province of China in 2008 [7], millions of cases have been reported, and CVA16 was responsible for nearly 40% of all the laboratory-confirmed HFMD [8,9]. Compared with EV71, another major causative agent of HFMD, the patients infected with CVA16 usually showed relatively mild symptoms followed by fever, maculopapular rashes on the skin of the hands, feet, and oral cavity, with a certain self-limiting trend [5,10]. However, severe complications occasionally occurred in individuals infected with CVA16, such as encephalitis [11,12], myocarditis and intractable shock [13], fatal pneumonitis [14], severe rhabdomyolysis, renal failure [15], and even death [13][14][15]. Apparently, CVA16 infection is not always a benign one, but so far, the molecular evolution and pathogenesis of CVA16 have not been fully described.
From the year 2009, we conducted and gradually strengthened virological surveillance on EVs in Nanjing area. The surveillance program is an effective tool in assessing the circulation of EVs, identifying EV outbreaks (especially HFMD epidemics), monitoring the emergence of novel lineages genetically distinct from currently recognized strains, and providing possible public health measures. However, until now, information regarding the pathogenic role of CVA16, their antigenic characteristics and molecular evolution in Nanjing area was limited. Here, we focused on the genetic variability and molecular evolution of CVA16 in a local area. Epidemiological profiles of 175 different CVA16 strains circulating in Nanjing spanning three years were evaluated based on VP1 sequence analysis.

Sample collection and ethical approval
Throat swab samples from patients with clinical diagnoses of HFMD or suspected EV infection were obtained from sentinel hospitals associated with the Nanjing Center for Disease Control and Prevention for EV surveillance. A standardized submission form was used for each specimen, which included basic information such as name, sex, age, residency, symptoms, date of onset, date of sample collection, type of case (hospitalized or outpatient), and clinical diagnosis. The guardians of all the participating children provided written consent for the study.

Extraction of RNA and EV detection
Throat swabs were mixed thoroughly in 4 mL phosphate-buffered saline (pH 7.4) to yield homogeneous suspensions. Viral RNA was extracted from 140 μL suspension using a QIAamp viral RNA Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. Extracted RNA was stored at -70°C for later use. Virus nucleotides were screened by a fluorogenic polymerase chain reaction (PCR) kit screened for pan-enterovirus, EV71, and CVA16 (bioPerfectus technologies, Taizhou, China).

Reverse transcription PCR and nucleotide sequencing
Synthesis of cDNA was carried out as described previously [16]. An amplicon of the full-length VP1 region (891 bp) was then generated from each viral cDNA by PCR using Applied Biosystems Veriti Thermal Cycler (Life Technology, New York, USA). The specific primers were as follows: sense, 5'-CCTATTGCAGACATGATTGACCAG-3'; antisense, 5'-TGTTGTTATCTTGTCTCTACTAGTG-3'.
All products were bidirectionally sequenced via sequencing service supplied by the Shanghai branch of Life Technology (Life Technology, Shanghai, China).

Phylogenetic analysis
The entire VP1 nucleotide sequences of 175 CVA16-positive samples were determined and used for phylogenetic analysis. The nucleotide sequence data were trimmed and assembled with the LaserGene package (DNAstar version 7.1) (DNASTAR, Madison, USA). Another 138 complete VP1 encoding sequences of CVA16 were obtained from NCBI GenBank (data not published). Phylogenetic dendrograms were constructed with the neighbor-joining method, and the mean evolutionary distances were calculated using the Kimura-2 model with MEGA program (version 5.22, Tamura, Peterson, Stecher, Nei, and Kumar 2011). Bootstrap testing with 1,000 replicates was used to estimate the reliability of the phylogenetic analysis. Bootstrap values below 80% were hidden.

Alignment of antigenic positions of VP1 protein
Ninety randomly selected local CVA16 strains and representative strains of each genotype were evaluated. Six reported antigenic epitope sequences were aligned using MEGA program (version 5.22).

Nucleotide sequence accession numbers
The 175 VP1 sequences of CVA16 strains reported here have been submitted to the GenBank database with the accession numbers KP751450 to KP751624.

Results
HFMD occurred throughout the whole period between 2011 and 2013 in Nanjing area. No significant outbreaks were observed in this period. All of the cases were sporadic, and the morbidity showed the same tendency each year, rising from April and with a peak between May and July. From the local surveillance data, a total of 4,440 clinically diagnosed HFMD specimens were collected within three days after the suspected symptoms appeared; 2,203 (49.61%) samples were determined to be EV-positive in the laboratory. CVA16-positive samples represented 25.15% (554/2,203) of these laboratory-confirmed EV-positive HFMD cases, and a total of 175/2,203 (7.94%) CVA16positive samples were randomly selected to be evaluated in the present study. Of these cases, all the patients were children, ranging in age from 3 months to 17 years, and more than 90% were under 5 years of age. A male predominance was noticed, with a male-tofemale ratio of 1.51:1.

Phylogenetic analysis of CVA16 strains
Complete sequences of the VP1 gene (891 nucleotides, positions 2,446-3,336, relative to strain CVA16/G-10) from the 175 different CVA16 specimens were then examined. A phylogenetic tree constructed based on the VP1 sequences demonstrated the existence of at least three genetically distinct groups ( Figure 1). The nucleotide divergences between the three groups were 6.4% (Gp1/Gp2), 9.9% (Gp1/Gp3), and 9.7% (Gp2/Gp3). CVA16 strains from 2011 to 2013 were scattered in these three groups.   However, only one strain from 2012 and none of the strains from 2013 belonged to Group 2. It was suggested that the three genetically distinct clusters cocirculated in the local area in a period spanning three years, with Group 2 no longer dominant since the year 2012.
To investigate the genetic characteristic of the CVA16 strains circulating in Nanjing, another 90 CVA16 strains isolated in other areas of China together with 48 international CVA16 sequences were retrieved from GenBank database and analyzed (Table 1). It is known that CVA16 can be grouped into genotype A and B. The prototype G-10 is the only strain of genotype A, and genotype B consists of two subgenotypes, B1 and B2. The B1 subgenotype of CVA16 can be further divided into three clusters: B1a, B1b, and B1c [17]. Here, the phylogenetic analysis showed that all 175 local CVA16 strains were grouped into the B1 genotype ( Figure 2). No B2 strain was detected. Within the B1 group, 29/175 (16.57%, Gp3) CVA16 strains segregated into a B1a cluster, and 146/175 (83.43%, Gp1+Gp2) strains segregated into a B1b cluster, indicating that the endemic prevalent subgenotypes of CVA16 strains in Nanjing are B1a and B1b, consistent with the subgenotype distribution in China. With respect to the geographic distribution of these cases, the B1b subgenotype came from all the 12 counties of Nanjing city (including Gulou, Xuanwu, Qinhuai, Jianye, Pukou, Qixia, Dachang, Yuhua, Jiangning, Luhe, Lishui, Gaochun), while the B1a subgenotype only covered 9 counties (except Qinhuai, Pukou, Gaochun). According to the amount and geographic distribution, the B1b subgenotype is much more widespread than B1a in Nanjing area.
Compared with CVA16 strains isolated worldwide, the Nanjing local CVA16 strains showed a closer relationship to the strains isolated in China, with a 94.0% identity, but only displayed a 90.8% identity with those strains abroad (Figure 2). In the B1a subgenotype, the 29 local strains had a higher homology with those isolated in Henan, Zhejiang, Beijing, Anhui, Guangzhou, and Yunnan from 2010 to 2013, but a little far from the ones isolated in Shenzhen from 2001 to 2003 and Shandong, Gansu from 2007 to 2008. Within the B1b subgenotype, Group 1 and Group 2 showed different similarities to separate geographic origins. Group 1 was closer to the ones isolated in Guangzhou, Henan, Yunnan, Wuhan and other areas from 2009 to 2013, and Group 2 shared a higher similarity to part of the Wuhan strains and those isolated in Anhui from 2010 to 2011. Other B1b strains isolated from 2007 to 2008 in mainland China were appreciably apart from Nanjing local strains, forming a relatively independent cluster. It is obvious that the relationship between homology and isolation time is closer and more concentrated than the geographic origins; thus, 2008 to 2009 could have been a pivotal period in CVA16 evolution in which apparent variation may occurred, leading strains isolated before or after this period to separate clusters in the phylogenetic tree.
However, although different clustering manners have been observed in different years and geographic locations, the total CVA16 strains manifested extreme similarity to each other (Table 1). Overall, all the CVA16 strains involved in this study shared a high nucleotide identity (93.8%). The similarity within the B1a and B1b clusters was 95.2% and 96.1%, respectively, and with an 8.8% divergence between the two clusters. As for the 175 Nanjing local strains, identity within the B1a and B1b genotype was 96.5% and 96.4%, respectively, and the divergence between local B1a and B1b was 9.9%, both a little higher than the ones calculated on the total strains. Comparatively, the divergence of amino acids in the VP1 region was much smaller, all less than 1%, indicating that evolution of CVA16 is highly conserved both locally and globally.

Amino acid alignment in antigenic regions of the VP1 protein
The protein coded by VP1 sequence is a critical component of CVA16 viral capsid, which is the most external and immunodominant of all the four capsid proteins (VP1, VP2, VP3, and VP4) [18]. Recently, six conserved regions in the VP1 protein have been  B1b, B1c, and B2) were randomly selected, and the prototype G-10 strain was used as a reference. The EV71 prototype strain EV71/BrCr and one representative strain of the local predominant EV71 subgenotype C4a, NJ01/JS/CHN/2010, were used as outgroup. As indicated in Figure 3, antigenic Region 4 was the most conserved sequence in all the CVA16 genotypes with no amino acid substitution.
Remarkably, the sequences in the same region of EV71 prototype and C4a subgenotype were exactly consistent with CVA16 strains, strongly suggesting a possible cross-reaction between EV71 and CVA16 in this region. Region 3 was also highly conserved, and only one amino acid change was observed in local strain PK11-040/NJ/CHN/2011. For Region 2 and Region 5, compared with the prototype G-10 strain, five mutations were found in all the other genotypes at positions 119 (M→L), 213 (A→E), 215 (P→L), 217 (S→A), and 220 (A→L). Noticeably, sporadic substitutions occurred more frequently in Region 1 and Region 6, especially among local strains in 2011-2012, suggesting that Region 1 and Region 6 require more study if they are to be considered as candidates for CVA16 vaccine development.

Discussion
CVA16 is a main causative agent of HFMD and is associated with a significantly high transmission among infected children. It belongs to Human enterovirus A, whose coding region is divided into three sub-regions: P1, P2, and P3. The P1 region encodes four structural proteins VP1, VP2, VP3 and VP4, whilst the P2 and P3 regions encode the nonstructural proteins including 2A, 2B, 2C, 3A, 3B, 3C, and 3D. The VP1 protein is located at the external side of the capsid and is constantly exposed to host immune pressure [20]. It has been demonstrated that the VP1 region is closely related to EV serotypes, and has been well accepted and extensively used as a target gene for EV genetic profile description and phylogenetic analysis [18,21,22]. Here, we used this method to investigate 175 CVA16 strains associated with HFMD occurred in Nanjing city, China, from 2011 and 2013, and established the molecular epidemiologic characteristics of the local circulating viruses. Phylogenetic analysis based on complete VP1 sequences has preliminary clustered these local 175 CVA16 strains into three groups (1, 2, and 3), suggesting that the local CVA16 strains may evolve through three distinct evolutionary chains.
It has been indicated that all the CVA16 strains so far can be clustered into genotypes A and B, and further into three subgenotypes (A, B1, and B2), based on the phylogenetic analysis of complete VP1 sequences [23]. The profiles of phylogenetic strain segregation in this study were generally concordant with this criterion for CVA16 genotyping, which revealed that the total Nanjing local CVA16 strains belonged to the B1 cluster. Separately, Group 1 and 2 classed into genotype B1b and Group 3 segregated into genotype B1a, while the other member of B1 cluster (B1c) has not been detected in the defined area, indicating that B1a and B1b are predominant CVA16 subgenotypes circulating in Nanjing area. It has been reported that the subgenotype B1a and B1b strains were first identified in 1995 and 1998, respectively, in Japan, and have played a predominant role in HFMD outbreaks worldwide since 2000. As for China, despite the existence of several CVA16 genotypes circulating globally, B1a and B1b are the only two subgenotypes of CVA16 present at home so far [23]. However, different subgenotypes have been observed in neighboring countries such as B2 in Japan [24] and B1c, B2 in Malaysia [25], suggesting a risk of import transmission of these adventive genotypes. Long-term and sustainable surveillance on HFMD helps to monitor the variation of existing pathogens and presence of new ones, and needs to be carried out continually and effectively.
It has been reported that the mean evolutionary rate (synonymous substitutions per site per year) determined for CVA16 B1a and B2 was slightly lower than the value derived from EV71 [21,26]. In this study, Nanjing CVA16 strains also exhibited high nucleotide sequence identities in complete VP1 sequences, both within or between the subgenotype B1a and B1b strains, with differences of less than 4% and 9%, respectively. The divergence of amino acids in the VP1 region is even as low as less than 1% regardless of within or between B1a and B1b groups, exhibiting high conservation in the evolutionary process. In addition, the CVA16 VP1 sequences in this study did not show successive accumulation of non-synonymous mutations through time, which was in agreement with the fact that CVA16 did not cause outbreaks in local communities for a period of time in Nanjing.
For decades, an autoimmune mechanism has been considered to be associated with EV infection [27,28]. Finding appropriate epitopes that elicit only neutralizing and not deleterious antibodies as vaccine candidates turned to be a tough task in the vaccine research field. It has been determined that the neutralization epitopes of EV are mainly located on capsid proteins, especially protein VP1 [23,29]. Several neutralizing epitopes on the VP1 structural protein of EV71 have been reported [30,31], but those capable of eliciting CVA16-neutralizing antibodies remained to be elucidated for a long time. Recently, a report indicated the identification of six neutralizing linear epitopes of CVA16 on the VP1 protein [23]. According to our alignment, the sequence of aa 187-201 was the most conserved among both CVA16 and EV71 strains, which may cause antigenic cross-reaction. The regions of aa 109-123, aa 163-177, and aa 211-225 would become good candidates when considered in CVA16 vaccine development because these sequences are highly conserved among all the CVA16 strains except CVA16 A genotype, which has not shown up for quite a long time. The sequences of aa 94-108 and aa 271-285 seemed not so suitable for CVA16 vaccine research since the nucleotide substitutions occurred frequently in these regions, and large amounts of CVA16 strains need to be evaluated to verify the possibility of being candidate for vaccine development.
Among the 175 strains investigated in this study, 23 (13.14%) were from inpatients with strain names beginning with NJ, and 152 (86.86%) were from outpatients. The inpatients presented with more severe symptoms compared to the outpatients; however, no distinctive genetic features of the VP1 region was observed in isolates obtained from inpatients.

Conclusions
Genotypes B1a and B1b have been steadily circulating with low divergence in Nanjing from 2011 to 2013. The information in this study will serve as the starting point for the molecular surveillance of CVA16 in Nanjing area. Although the number of CVA16 isolates was limited, the observation has nevertheless showed the genetic background of CVA16 circulating in Nanjing and would assist with both individual patient diagnosis and public health actions. To provide further evidence of epidemiological linkage and evolutionary recombination events in CVA16, substantial efforts to improve the EV surveillance are necessary.