Molecular Characterization and Phylogenetic Analysis of Porcine Epidemic Diarrhea Viruses Associated with Outbreaks of Severe Diarrhea in Piglets in Jiangxi, China 2013

Porcine epidemic diarrhea (PED), caused by porcine epidemic diarrhea virus (PEDV), is a highly contagious, acute enteric viral disease of swine characterized by vomiting, watery diarrhea, dehydration and death. To identify and characterize the field PEDVs associated with the outbreaks of severe diarrhea in piglets in Jiangxi, 2013, the complete genome sequences of two representative strains of PEDV, designated CH/JX-1/2013 and CH/JX-2/2013, were determined and analyzed. The genome sequences of both emergent Jiangxi PEDV strains, CH/JX-1/2013 and CH/JX-2/2013, were 28,038 nucleotides in length excluding 3’ poly (A) tail. Compared to the PEDV CV777 strain, CH/JX-1/2013 and CH/JX-2/2013 had some unique genetic characteristics in the proximal region of the 5´-UTRs. Phylogenetic analysis of the complete genomes and the structural proteins revealed that CH/JX-1/2013 and CH/JX-2/2013 had a close relationship with post-2010 Chinese PEDV strains and US strains identified in 2013. The nucleotide identity between the two Jiangxi strains (CH/JX-1/2013 and CH/JX-2/2013) and 30 strains of PEDV identified ante-2010 and post-2010 ranged from 96.3–97.0% and 97.3–99.7%, respectively. Multiple nucleotide and deduced amino acid mutations were observed in the ORF1a/b, S, ORF3, E, M and N genes among the current field PEDV strains when compared to the CV777 strain. Some of the mutations altered the amino acid charge and hydrophilicity, and notably, there was an amino acid substitution in the middle of one neutralizing epitope (L1371I) of the S gene of both CH/JX-1/2013 and CH/JX-2/2013. Taken together, the accumulated genetic variations of the current field PEDV strains might have led to antigenic changes of the viruses, which might confer the less effectiveness or failure of the CV777-based vaccines currently being widely used in Jiangxi, China.


Ethics statement
The ethics committee of Jiangxi Agricultural University approved the animal protocol for this study (protocol number P-2013-03). All the procedures involving animals in this study were carried out in accordance with The Care and Use Guidelines of Experimental Animals established by the Ministry of Agriculture of China. Intestinal and fecal samples were collected according to the approved procedures. The intestinal samples used in this study were obtained from the dead piglets and the fecal samples were non-invasively collected immediately after excretion from healthy and diarrheal pigs from premises with PED outbreaks in Jiangxi, China (S1 Table).

Virus identification
Fecal and intestinal samples (N = 125) were collected from < 10-day-old suckling piglets with severe watery diarrhea, vomiting and dehydration from different premises in Jiangxi Province, China (25-29°N and 114-118°E) in 2013. All the samples were examined by a RT-PCR established in our laboratory. Briefly, the total RNA was extracted from the samples using RNAplus Reagent (TaKaRa, Japan) following the manufacturer's instructions. The one-step RT-PCR for determination of PEDV-positive samples was carried out using 50-200 ng of extracted RNA, forward primer (5´-GTATTGGTGGTGAGCGGAAT-3´), reverse primer (5´-CCTGTTCC GCCATTCTATCA-3´) and OneStep RT-PCR kit (Qiagen, Valencia, CA, USA) according to the manufacturer's protocol. To eliminate possible co-infections caused by porcine transmissible gastroenteritis virus (TGEV), porcine rotavirus (PoRV), and other common pathogenic intestinal pathogens, differentiation assays were performed using standard protocols. Ninety-six out of 125 samples were tested positive for PEDV, and all samples had been confirmed to be free of TGEV, PoRV and other common enteric pathogens (data not shown). To address the genetic/antigenic variations, and phylogenetic characteristics of PEDV strains associated with 2013 Jiangxi PED outbreaks, two representative PEDV strains, designated CH/JX-1/2013 and CH/JX-2/2013 (the GenBank accession numbers are KF760557 and KJ526096, respectively), were used for sequencing the full-length genome.

Full-length genome sequencing
Total RNAs were extracted from the feces and small intestinal homogenates by RNAplus Reagent (TaKaRa, Japan) according to the manufacturer's instructions. The concentrations of the extracted RNAs were measured by NanoDrop 2000 spectrophotometer (Thermo Scientific, USA) and then stored at −80°C until use. The first-strand cDNA synthesis was performed at 42°C for 50 min and then 95°C for 5 min to inactivate the M-MLV reverse transcriptase (TaKaRa, Japan) and followed by 4°C for 5 min. The entire PEDV genome was amplified by 33 pairs of primer designed with Primer 3 software (http://primer3.ut.ee/) based on the conserved regions determined by a multiple alignment analysis of the reference strains and CV777 (S2 Table). Fragments were amplified on the conditions of a denaturation at 94°C for 4 min, 35 cycles (94°C x 45 sec, 53°C x 45 sec, 72°C x 1.5 min), and then with a final extension at 72°C for 10 min. PCR products obtained were subjected to gel purification using a gel extraction kit (TaKaRa, Japan), and afterwards cloned into pMD 18-T vectors (TaKaRa, Japan) following the manufacturer's protocols. Five positive clones of each amplicon were submitted to a commercial sequencing company (Sangon Biotech, Shanghai, China) for sequencing at both directions by Sanger sequencing methodology. The 5'-and 3'-RACE for the determination of the terminal sequences of both CH/JX-1/2013 and CH/JX-2/2013 were performed by using 5'/3' SMARTer RACE kit (Clontech, Beijing, China) following the manufacturer's instructions.

Sequence and phylogenetic analysis
The raw sequence fragments were imported to SeqMan in DNAStar Lasergene V 7.10 (DNAStar, Inc., Madison, WI) for assembly and annotation. Nucleotide and deduced amino acid (aa) sequences of both CH/JX-1/2013 and CH/JX-2/2013 and 30 reference PEDV sequences retrieved from GenBank were comparatively analyzed. A summary of the background information of PEDVs used in this study is shown in Table 1. The complete genome sequences of CH/JX-1/2013 and CH/JX-2/2013 were deposited into GenBank. Phylogenetic trees based on the entire genomes, and deduced aa sequences of S, ORF3, E, M, and N genes were constructed using the neighbor-joining method of MEGA 5.2.2 (http://www.megasoftware.net/) with a bootstrap of 1,000 replicate datasets. Primary sequences of 5´-proximal region of 5´-UTRs (nt 42 to 133) were pairwise compared between the two Jiangxi strains (CH/JX-1/2013 and CH/JX-2/2013) and the reference strains. Antigenicity and hydrophilicity analyses based on the aa sequence from 1 to 350 at N-terminal of the S proteins were carried out by Protean software of DNAStar Lasergene V7.10 (DNAStar, Inc., Madison, WI).   Table). A total of 26 unique nucleotide substitutions were identified between two Jiangxi PEDV strains (CH/JX-1/2013 and CH/JX-2/2013) and 30 reference PEDV strains, which resulted in 14-aa changes, and most of them were located in ORF1a, ORF1b and S gene. Two of the nucleotide substitutions led to aa changes along with the charge variations ( Table 2). The phylogenetic tree based upon the full-length genome sequence of CH/JX-1/2013, CH/JX-2/2013 and 30 reference PEDVs indicated that the PEDV strains could be divided into two groups, designated for group1 (G1) and group 2 (G2): G1 could be further divided into three subgroups, i.e., 1a, 1b and R, a tentative cluster consisting of virulent DR13 isolated in South Korea in 1999 and CH/S isolated in China in 1986; and G2 was split into two subgroups, 2a and 2b (Fig. 1A). Notably, all the strains identified from 2011 to 2013 were fallen into G2; while the prototype of CV777, together with three Korean strains (SM98, virulent DR13 and attenuated DR13), and four earlier Chinese strains (JS2008, JS2008new, LZC and CH/S) were fallen into another group, i.e., group 1. Genetic characteristics were observed between the two groups: 1) compared to genome sequences of the members in G1, four insertions, 20803G, 20810CAGGGTGTCAA20820, 20830G, 21042AAT21044 and two deletions, 20842A, 21097CGTGAT21102, existed in the N-terminal domain (NTD) of the S protein in G 2 members; 2) the three field PEDV strains of JS2008, JS2008new and SD-M together with two attenuated PEDV strains, DR13 and vaccine_KC189944, were clustered into subgroup 1b. All of the five PEDV strains aforementioned had 24-nt deletions in nsp3 of ORF1a and 49-nt deletions in the C terminus of ORF3; 3) CH/JX-1/2013 and CH/JX-2/2013 along with three Chinese strains (GD-B, JS-HZ2012 and BJ-2011-1) were clustered into an independent clade, and CH/JX-1/2013 and CH/JX-2/2013 had 99.7%, 99.7% and 99.6% nucleotide identity with the three Chinese strains, respectively. These strains shared six additional unique nucleotide substitutions (T227A, C2342T, G2346T, T2724C, T6331C, C7233T) with the rest of reference PEDV strains. Of which, one led to an aa change (T682M in ORF1a, from hydrophilic polarity T to hydrophobic polarity M).

Comparative analysis of structural genes
The full-length of S gene of CH/JX-1/2013 and CH/JX-2/2013 was 4,161 nt in size, which was 9nt longer than that of the prototype of PEDV CV777 strain. The results from nucleotide sequence comparisons of the S gene of 34 strains of PEDV, including CH/JX-1/2013, CH/JX-2/2013 and other 32 reference PEDV strains from China, United States, UK, South Korea, Belgium, France and Japan showed that the two Jiangxi strains had a 99.8% nucleotide identity with each other, and 96.7% nucleotide identity with CV777. CH/JX-1/2013 and CH/JX-2/2013 shared 97.3-99.6% nt identity and 93.1-99.1% aa identity with the post-2010 PEDV strains, respectively. By contrast, the two Jiangxi strains showed only 93.5-95.0% nt identity, and 92.5-94.8% aa identity with the ante-2010 PEDV strains and attenuated vaccine strains, respectively (S4 Table).
The phylogenetic tree based on the aa sequences of the S protein of 34 strains of PEDV showed that they could be classified into two major groups, and each group contained two subgroups (Fig. 1B). Both CH/JX-1/2013 and CH/JX-2/2013 belonged to subgroup 2b, which also included five newly identified US strains (13-019349, Co/13, IA1, ISU13-22038-IA-homogenate and ISU13-22038-IA-P9) and 13 Chinese strains identified at the same period or later  than 2010. Compared to G1, there were two insertions (61VNST64 and 136N) and one deletion (155DG156) in the G2 members. Amino acid differences were also present between G1 and G2, and most of them were located in the area of aa 1 to 350 at N-terminal of the S protein ( Table 3). The analysis based on aa from positioned 1 to 350 of the S protein suggested that the antigenicity and hydrophilicity of S protein of CH/JX-1/2013, CH/JX-2/2013 might have already changed due to the mutated amino acids with changes of charges and polarities (Fig. 2).
Notably, an amino acid substitution was found in the middle of one neutralizing epitope Phylogenetic trees based on the complete genome, aa sequences of structural proteins and ORF3 of PEDV strains. The trees were constructed by the distance-based neighbor-joining algorithm using MEGA 5.2.2 software. Bootstrap was set in 1,000 replicates with a value >70% to assess the significance of the tree topology. A bar of 0.002/0.005 indicates nucleotide or amino acid substitutions per site. "•" indicates the strains identified in this study, "" indicates the strains from China, "◆" indicates the strains from Belgium, "■"indicates the strains from USA, "▲"indicates the strains from South Korea, "▼"indicates the strains from France, "◇"indicates the strains from Japan. 1A: Phylogenetic tree generated on the basis of nucleotide sequences of the complete genome of 33 PEDVs. 1B to 1F: Phylogenetic trees based on deduced amino acid sequences of S glycoprotein genes, ORF3, envelope, membrane, and nucleocapsid genes, respectively.  Table 4). The phylogenetic tree based on the deduced aa sequences of the ORF3 of the earlier European strains, CV777, Korean strains, US strains and Chinese strains identified ante or post-2010 revealed that those strains were grouped into two different groups (Fig. 1C). CV777 and two early strains SM98 and LZC, and two cell-adapted strains SD-M and attenuated DR13 were classified into group 1; while the 24 strains, including CH/JX-1/2013, CH/JX-2/2013, 15 Chinese PEDV strains, a Korean strain (virulent DR13), five newly determined US strains, and      (from negative charge/hydrophilic polarity E to neutral charge/ hydrophobic polarity Q) of the N terminus of the M protein, which might have an impact on its antigenicity/immunogenicity. The entire N gene of CH/JX-1/2013 and CH/JX-2/2013 was 1,326 nt long, encoding a protein of 441-aa. Sequence analyses revealed that the aa homology between CH/JX-1/2013 and CH/JX-2/2013 and other PEDV strains used in this study varied from 95-98.9%. The nucleotide and amino acid sequences of N protein were highly conserved and neither insertion nor deletion was found in all these PEDV strains analyzed in the study except a few sporadic mutations. Multi-alignment results of the aa sequences of the PEDV N proteins indicated that all the strains could be divided into two groups, i.e., Group 1 and Group 2. Each group contained two subgroups (Fig. 1F). The phylogenetic topology of N protein was similar to that of ORF3. Three specific aa changes were present between group 1 and group 2, at the residue 142 (A142T), 242 (H242L, from hydrophilic polarity I to hydrophobic polarity T) and 397 (Q397L). Five aa substitutions at the residue 84 (G84A), 205 (N205K, from neutral charge N to positive K), 381(L381P), 395 (L395Q, from hydrophobic L to hydrophilic Q), and 398 (H398N, from hydrophilic H to hydrophobic N) existed in subgroup 2a when compared with other three subgroups, i.e., subgroup 2b, subgroup1a and subgroup 1b.

Nucleotide sequence comparison and structure prediction of 5´-RNA terminus
The 5´UTRs of CH/JX-1/2013 and CH/JX-2/2013 were 292 nt in length, which was 4-nt shorter than that of CV777, resulting from 5-nt deletion and one nt insertion in the proximal region of 5´-UTRs of the two Jiangxi strains (Fig. 3). The core sequence (CUAAAC) of the PEDV leader transcription-regulating sequence (TRS) was extremely conserved with no nucleotide substitutions in all the PEDV strains. An 'A' deletion was observed immediately following the core sequence of these PEDV strains except the strains of CV777 and LZC. In comparison with CV777 and early isolated strains of SM98 and LZC, a 'U' insertion in loop 2 was found in both CH  in other four PEDV strains isolated earlier (CV777, Attenuated vaccine, LZC, and SM98). Compared to CV777, a 4-nt deletion (UUCC) existed in the stem region of SL4 of CH/JX-1-/2013, CH/JX-2/2013 and other six PEDV strains. However this deletion would not impact the secondary structure of the SL4 as demonstrated by an analysis (data not shown). The 3'UTRs of both CH/JX-1/2013 and CH/JX-2/2013 were 334 nt in size, and no insertion and deletion was observed (Table 2).

Discussion
PED affected massive pig farms in Jiangxi, China in 2013, causing substantial economic losses in the pig industry of Jiangxi. To elucidate the molecular characterization and phylogeny of the PEDV field strains associated with 2013 Jiangxi PED outbreaks, the entire genome sequences of two representative PEDV strains confirmed in Jiangxi were determined and analyzed. Phylogenetic analyses demonstrated that CH/JX-1/2013 and CH/JX-1/2013 defined in this study along with the strains determined post-2010 were clustered into group 2, whereas the strains detected ante-2010 were clustered into group 1. The results of phylogenetic analyses of those PEDV strains examined in this study were similar to that of described elsewhere [22]. CH/JX-1/2013 and CH/JX-1/2013 showed the highest nucleotide identity (>99%) with six newly confirmed US strains (IA1, Co/13, USA/Indiana, ISU13-22038-IA-homogenate, ISU13-22038-IA-P9, and 13-019349) and five Chinese strains (GD-B, AH2012, BJ-2011-1, CH/ZMD/ZY, and CH-HZ/2012) identified post-2010, suggesting these PEDV strains might have evolved from the same origin although the mechanisms of the evolution of the viruses are roughly unknown yet. Three field PEDV strains, JS2008, JS2008new and SD-M together with two attenuated PEDV strains, DR13 and vaccine_KC189944, were clustered into an independent cluster and showed 24-nt deletions in nsp 3 of ORF1a and 49-nt deletions in the C terminus of ORF3. These unique characteristics suggested that the three field PEDV strains might derive from the same source, and a recombination event might have occurred between these three viruses and the two vaccine strains in the same subgroup. Interestingly, a comparison analysis revealed that CH/JX-1/2013 and CH/JX-1/2013 had 99.7% nucleotide identity to the strain of GD-B (GenBank accession NO.JX088695) isolated from Guangdong, an adjacent province of Jiangxi, in which the severe PED outbreaks emerged in October 2010 [11,23]. Pigs are frequently traded between Jiangxi and Guangdong, which might be one of important factors causing the cross dissemination of PEDVs in this region. However, further study needs to be performed.
The findings from this study demonstrated that the 5´-UTR and ORF3 protein of CH/JX-1/2013 and CH/JX-2/2013 and all other PEDV strains analyzed were highly conserved, which were consistent with the studies reported previously [5,9]. In general, strict conservation of these regions is essential for the PEDV life cycle. It is known that the 5´-UTRs of coronaviruses form conserved RNA structural elements which are critical for viral replication, sgmRNA transcription, and translation [24]. Studies on the betacoronavirus have indicated that the stemloop 2 (SL2) in 5´UTR is extremely conserved in all coronaviruses, and plays a cis-acting reaction on transcription and replication, and more importantly the replication of the virus requires the conservative sequence and a firm number of nucleotides with specific properties in SL2 in 5´UTR [25]. Although the 'U' insertion in SL2 does not alter the RNA secondary structure, it might slightly affect the efficiency of virus replication, which needs to be addressed in the future studies. The core sequence (5´-CUAAAC-3´) in TRS of PEDV, which was also present in CH/JX-1/2013 and CH/JX-2/2013, has reported to be a determinant factor in transcriptional regulation in coronavirus because the synthesis of sgmRNA requires the appropriate tertiary structure of core sequence [26]. Both PEDV and TGEV belong to the alphacoronavirus genus within the Coronaviridae family and possess the conserved core sequence structure.
Studies based on the infectious genomic TGEV cDNAs have proved that the nucleotides immediately flanking the TRS sequence could apparently affect the expression of the sgmRNA. And the nucleotides adjacent to core sequences of the leader TRS (CS-L) by the 3´region are more decisive for mRNA synthesis than nucleotides in the 5´region. The motif of 5´-GAAA-3´within 3´CS-L (5´-CUAAACGAAA-3´) shows a higher expression than that of other motifs [27]. In respect of the PEDV, an 'A' deletion immediately followed the CS-L sequence in the 5´UTR of CH/JX-1/2013 and CH/JX-2/2013 as well as other post-2010 PEDV strains analyzed in the study and thus, this deletion formed a four base oligonucleotide (5´-GAAA-3) adjacent to the CS-L by the 3´region, which might up-regulate the sgmRNA expression in PEDV [24].
The S protein of PEDV is known to play pivotal roles in viral entry and inducing the neutralizing antibodies in natural hosts, and thus makes it become a primary target for the development of effective vaccines against PEDV [28][29][30][31]. It was demonstrated that there were significant genetic variations in the S gene between the newly determined PEDV field strains and early isolates [20,32]. The results in this study agreed with previous documentations. Additionally, CH/JX-1/2013 and CH/JX-2/2013 showed extremely high nucleotide and aa identity with each other but less to CV777 and attenuated strains (attenuated vaccine_KC18944 and attenuated DR13). There were significant aa differences between the G1 and G2 members, and most of them were located in NTD of the S protein of PEDVs. When compared to the group 1 strains of PEDV, CH/JX-1/2013 and CH/JX-2/2013, members of group 2, showed significant antigenic and hydrophobic differences, and some of which were located in the neutralizing epitope regions. Those genetic variations made the CH/JX-1/2013 and CH/JX-2/2013 and post-2010 field PEDV strains different from the ante-2010 strains, especially the CV777. It might explain why the recent PED outbreaks were significantly different from the previous sporadic outbreaks. Present PED became a devastating enteric viral disease causing substantial economic losses in the pig industry in major pig-raising countries in the world, and recently made it become a reportable disease by USDA [33]. Moreover, the present PEDV field strains displayed considerable genetic variations from CV777-based vaccine strain, a strain being widely used for production of PEDV vaccines in China for many years. Notably, Leucine, a highly hydrophobic aa residue in the middle of one neutralizing epitope (L1371I) of the S gene of CV777 was substituted by an Isoleucine, a higher hydrophobic index aa residue in both CH/JX-1/2013 and CH/JX-2/2013 strains. This variation might provide a possible mechanism for a poor protection on swine vaccinated with the CV777-based vaccines.
Park et al. [29] have reported that a reduced cell-culture-adapted PEDV strain (attenuated DR13, passage 100) has a truncated ORF3 of a 51-nt deletion, suggesting that this gene may be involved in cell tropism and is essential for the virulence of PEDVs. By contrast, CH/JX-1/2013 and CH/JX-2/2013, together with the post-2010 variant PEDV strains had an intact ORF3 of 675 nucleotides encoding a protein of 224 aa. The phylogenic analysis demonstrated that the ORF3-based tree had a similar topology with the one generated from the entire genome sequences. And thus the ORF3 of PEDV may serve as a useful target gene for phylogenetic analysis of newly emerging field PEDV strains since its small size. As recognized previously [34,35], there is a large deletion region existing in PEDV isolates attenuated DR13 (51-nt at the position of 245-295), SM98 (210 nt at the position of 1-210) and SD-M (51-nt at the position of 245-295), a Chinese strain propagated only four passages on cell lines.
In summary, the findings obtained in this study provide some insight into the genetic-/phylogenetic variations and molecular characterizations of the Jiangxi field PEDV strains associated with the outbreaks in piglets in Jiangxi, China 2013. The comparisons of the fulllength genome and the structural protein genes and deduced aa sequences revealed that two Jiangxi strains CH/JX-1/2013 and CH/JX-2/2013 defined in this study had a close relationship with the recent prevailing field PEDV strains in China and the United States. Differences at the level of the nucleotides and deduced aa between the present field PEDV strains and CV777, especially the aa substitutions in the neutralizing epitope of PEDV field strains, might have conferred the less effectiveness of the vaccines currently widely being used in China. It might be urgently needed to develop improved efficacious and safe vaccines against field PEDVs being currently circulating in China.