Isolation and evolutionary analyses of porcine epidemic diarrhea virus in Asia

Virus detection and isolation

In 2018, an outbreak of diarrhea occurred in a CV777-vaccinated pig farm (numbers of sows ≥ 100) in Hubei Province in China. Many pigs in the farm suffered from severe watery diarrhea, and some of them died. Samples of intestinal tissues were collected from dead pigs and sent to the Veterinary Diagnostic Laboratory of Hubei Academy of Agricultural Sciences in Wuhan, China, for diagnosis. Tissues were immersed with Dulbecco’s modified Eagle medium (DMEM; Gibco, Grand Island, NY, USA), and were then homogenized using a QIAGEN TissueLyser II (QIAGEN, Dusseldorf, Nordrhein-Westfalen, Germany). The sample homogenates were then frozen at −80 °C and thawed for three times. After that, the supernatants were filtered through a 0.22-µm membrane and were harvested for RNA and virus isolation. Total RNAs were extracted using TRIzol (Thermo, Waltham, MA, USA) and were reverse transcribed to cDNA using a Thermo Scientific First Strand cDNA Synthesis kit (Thermo, Waltham, MA, USA). Viral nucleic acids were detected by RT-PCR assays using the cDNA as templates and the primers specific for PEDV (F: 5′-TTCGGTTCTATTCCCGTTGATG-3′, R: 5′-CCCATGAAGCACTTTCTCACTATC-3′), TGEV (transmissible gastroenteritis virus) (F: 5′-TTACAAACTCGCTATCGCATGG-3′, R: 5′-TCTTGTCACATCACCTTTACCTGC-3′) and PoRV (porcine rotavirus) (F: 5′-CCCCGGTATTGAATATACCACAGT-3′, R: 5′-TTTCTGTTGGCCACCCTTTAGT-3′), respectively.

Vero cells (Purchased from ATCC, Manassas, VA, USA) were used for virus isolation. In brief, homogenate supernatants and trypsin (5 µg/ml) were inoculated into monolayers of Vero cells, which were then incubated in a 37 °C incubator supplemented with 5% CO₂. Cells with obvious cytopathic effects (CPEs) were harvested, thawed, and refrozen multiple times. The harvested virus suspension was then inoculated into newly prepared Vero cells for passages, and the propagation was continuously performed for 20 passages (F20). Virus RNA was extracted every five passage for RT-PCR detection of the virus nucleic acids.

Virus titration and serum neutralization

Virus titers were measured on 96-well plates using 10-fold serial dilutions of culture supernatant in triplicate per dilution to determine the quantity of viruses required to produce CPEs in 50% of cells. After incubating for enough time, no more CPEs appeared, and TCID50 was calculated using the Reed-Muench method (Reed & Muench, 1983). The virus titer was also determined by plaque assay using Vero cells and expressed as plaque-forming units (PFU) per mL. The serum neutralization (SN) test was performed in 96-well plates with inactivated serum collected from the guinea pigs infected with the vaccine strain CV777. The virus was diluted in serum-free DMEM to make 200 TCID in a 50 µL volume, and mixed with 50 µL of 2-fold serial dilution serum. The mixture was added to cells cultured in 96-well plates and incubated at 37 °C for 1 h. After removing the mixture and thoroughly washing three times with PBS, the cells were incubated at 37 °C with 5% CO₂ for 2 days. Neutralization titers were calculated as the reciprocal of the highest dilution of serum that inhibits CPEs.

Genome sequencing and annotation

Genomic RNA was extracted using the TAKARA RNA extraction kit (Takara, Kusatsu, Shiga, Japan) following the manufacture instruction. The quantity and quality of the extracted RNA were measured by using a Nanodrop spectrophotometer (Thermo, Waltham, MA, USA). The RNA was then subjected to reverse transcription for cDNA using a cDNA synthesis kit (Thermo, Waltham, MA, USA). Genome sequencing was performed with a paired-end library constructed by using a NEB-Next^® DNA Library Prep Master Mix Set for Illumina (NEB, Ipswich, MA, USA) and subsequently sequenced on an Illumina NextSeq 500 with 2 × 150 paired end sequencing chemistry. After filtering, the clean reads were assembled using SPAdes v3.10.1 (Bankevich et al., 2012) and assembled sequences were mapped to the reference genome. The prediction of the genes and proteins were conducted with Prokka v1.12 and RAST Serve (http://rast.nmpdr.org) (Aziz et al., 2008). The complete genome sequence as well as its annotations were deposited into NCBI GenBank under the accession number MT166307.

Comparative genomics and bioinformatical analysis

The NCBI data was search for “porcine epidemic diarrhea virus” and a total of 207 complete genome sequences were publicly available for PEDV isolates representing different parts of Asia (See Table S1). All of these 207 sequences were downloaded for further analysis. The average nucleotide sequence identity between the genomes of HB2018 and CV777 was calculated by ANI calculator (Goris et al., 2007). Sequence alignments were performed using MAFFT v7.4.02 (Katoh & Standley, 2013). Nucleotide sequence similarity and the putative recombination sites was assessed by SimPlot v.3.5.1 (Lole et al., 1999), with a sliding window size of 500 bp, step size of 100 nucleotides, and 1,000 bootstrap replicates, using gap-stripped alignments and the F84 (ML) distance model. Phylogenetic trees based on complete genome sequences were generated by using MEGA X software with 1,000 bootstrapping (Kumar et al., 2018).The evolutionary history was inferred by using the Maximum Likelihood method and Tamura-Nei model (Tamura & Nei, 1993). Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using the Maximum Composite Likelihood (MCL) approach, and then selecting the topology with superior log likelihood value. The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. A maximum likelihood tree was also generated using the BEAST 2 package (version 2.6.3) (Bouckaert et al., 2019). Gamma correction for site heterogeneity and the GTR model (Gatto, Catanzaro & Milinkovitch, 2007) were selected for the tree generation. Both of the trees were annotated and visualized by using the iTOL v.4 online tool (Interactive Tree of Life, http://itol.embl.de/) (Letunic & Bork, 2019). Single nucleotide polymorphisms (SNPs) between two genome sequences were determined by the MAUVE package (version 2.4.0) (Darling et al., 2004), and the coding effect of these SNPs were analyzed using a previously reported local Perl command (Peng et al., 2016). Protein structure was generated using SWISS-MODEL (https://swissmodel.expasy.org). Protein N-glycosylation sites were predicted using online software (http://www.cbs.dtu.dk/services/NetNGlyc/). Threshold values of greater than 0.5 and Jury agreement 9/9 were used for the high-specificity N-glycosylation sites determination (Sagesser et al., 1997).

Isolation of PEDV HB2018 and its genomic characteristics

RT-PCR detection of the viral nucleic acids revealed that the intestinal samples from pigs suffered and died from severe watery diarrhea were positive for PEDV but negative for TGEV and PoRV (Fig. S1). Through virus isolation and purification using Vero cells and determination of PEDV nucleic acids using RT-PCR, a PEDV strain was finally recovered and designated HB2018. The TCID50/0.1 mL value of HB2018 was 10^5.3. The complete genome sequence of PEDV strain HB2018 was 28,138 bp in length. This 2.8-kb genome contained seven open reading frames (ORFs): ORF1a (nucleotide positions 281 to 12,634), ORF1b (positions 12,664 to 20,625), S gene (positions 24,782 to 25,456), ORF3 (positions 25,675 to 25,667), E gene (positions 25,437 to 25,667), M gene (positions 25,675 to 26,355), and N gene (positions 26,367 to 27,692).

Phylogenetic analysis based on the complete genome sequence showed that HB2018 was phylogenetically distinct from the vaccine strain CV777 (Fig. 1A). According to the genotyping system based on a full-length genomic sequence analysis (Guo et al., 2019; Wang, Fang & Xiao, 2016a), HB2018 and CV777 belonged to two different genotype: HB2018 was assigned as a type GII strain while CV777 was a GI strain (Fig. 1A). The average nucleotide identity between the genomes of HB2018 and CV777 (GenBank accession no. AF353511) was 96.06% (Fig. S2). The ORF1, ORF3, E, M, and N genes of HB2018 as well as their encoding proteins were highly homologous to those of CV777 (nucleotide identity ≥ 95% for genes; amino acid similarity ≥ 95% for proteins) (Fig. 1B; Table 1). However, the identity of the S genes and proteins between the two strains was relatively low: the homology for nucleotide and amino acid sequences between HB2018 and CV777 were 93.76% and 93.44%, respectively (Figs. 1B & 1C; Table 1). SNP analysis determined a total of 946 SNPs in the genome sequence of HB2018 when compared to the genome sequence of the reference strain CV777. Among these SNPs, 925 SNPs including 262 non-synonymous substitutions and 663 synonymous substitutions were located with the ORF regions, with an overall ratio of nonsynonymous to synonymous substitutions (dN/dS) of 0.39 (Table 2). The dN/dS ratios in each of the ORFs encoded by the HB2018 genome ranged from 0.15 to 0.61, with the S protein had the highest dN/dS ratio (Table 2).

Compared to the S protein of CV777, the S protein of HB2018 had changes, deletions, and/or insertions of amino acids at multiple sites (Table S2 and Fig. S3). Notably, most of these mutations occurred in the N-terminal domain (NTD, 19-233aa) of the S protein (Fig. 1C; Fig. S3). Interestingly, some of these mutations were located within the neutralizing epitopes of PEDV (COE (499–638), SS2 (748–755), SS6 (764–771) and 2C10 (1368–1374)). In addition, these mutations led to a structural change at some parts of the HB2018 S protein compared to the CV777 S protein (Figs. 1D & 1E).

Phylogenetic analysis of Aisan PEDV isolates

To explore the phylogenic relationships of the PEDVs currently circulating in Asia, we generated two maximum likelihood trees based on the whole genome sequences, either by using the MEGA X software with the Tamura-Nei model (Fig. 2A) or by using the BEAST 2 package with the GTR model (Fig. 2B). Both of the results revealed that the 208 PEDV strains in Asia, representing the 207 genome sequences publicly available in NCBI and the HB2018 sequence was divided into two genogroups: GI (classical) and GII (variant). Interestingly, isolates in China before 2010 and the vaccine strain CV777 were included within the GI genogroup. However, most of the PEDV isolates from China as well as the other Asian countries after 2010 belonged to GII genogroup (Figs. 2A & 2B). The phylogenetic trees also showed that the two genogroups consisted of several subgroups: the genogroup GI was divided into two subgroups, GI-a and GI-b, while the genogroups GII was divided into three subgroups, GII-a, GII-b, and GII-c (Figs. 2A & 2B). The GI-a and GI-b subgroups included isolates from China before 2010 and several Chinese isolates between 2010 and 2015 (Figs. 2A & 2B). Most of the GII isolates from China and South Korea and all GII isolates from Japan were included within the GII-a subgroup, while less proportion of the Chinese GII isolates and most of the GII-a isolates from Southeast Asia (Vietnam and Thailand) were included within GII-b subgroup (Figs. 2A & 2B). Interestingly, the GII-c subgroup only consisted of isolates from China (Figs. 2A & 2B).

By analyzing isolation years and genogroups of PEDVs, the history of PEDV and the evolution in China are speculated. Between 1986 and 2008, only five PEDV strains were sequenced in China, and all of them belonged to G1 (Figs. 2A & 2B). However, the number of PEDV sequences increased significantly after 2010 (Fig. 2C). While several PEDV sequences belonged to genogroup GI after 2010, most sequences from China were GII strains (Figs. 2A–2C).

Analysis on the S protein

Compared to the S proteins of the Chinese GI-a strains, amino acid changes, deletions, and/or insertions were observed at multiple sites within the S proteins of the Chinese GI-b strains (Table S3 and Txt S1). Compared to the S proteins of the Chinese GI strains, the S proteins of the Chinese GII strains commonly had amino acid changes, deletions, and/or insertions at several sites (Table S4 and Txt S1). Most these mutations occurred in S-NTD (19-233aa) of the S protein (Fig. 3A; Txt S1).

Compared to S proteins of the Chinese GII-a strains, S proteins of most isolates from Japan, South Korea, and Vietnam did not contain characteristic amino acid mutations, with the exception of S proteins of two Japanese strains (NIG-2/JPN/2014 and KMM-1/JPN/2014) which had a continuous deletion of 194 amino acids at sites 23–216 (Fig. 3B; Txt S2). In addition, four isolates from Japan (GenBank accession numbers LC063844 –LC063847) and three isolates from South Korea (KNU-141112-S DEL5, KNU-141112-S DEL5ORF3, KNU-1406-1) had a continuous deletion of 5 amino acids (GENQG) at sites 56-60 in their S proteins compared to the S proteins of the Chinese GII-a strains (Fig. 3B; Txt S2). In addition to the GII-a strains, several isolates from China, South Korea, Thailand, and Vietnam were GII-b strains (Fig. 1). Compared to S proteins of most of the Chinese GII-b strains, S proteins of the GII-b isolates from Thailand and Vietnam had amino acid changes at sites 130-131 (SI →DN), 182 (Y →H), 287 (I →M), 324 (N →D), 327 (S →A), 358 (A →T), 367 (I →T), 433 (D →G), 558-559 (TN →PT), 1287 (E →K), and 1317 (L →F) (Txt S3).

The N-glycosylation sites in the S proteins of the Asian strains studied were investigated herein. Bioinformatical analysis revealed that most of the 208 Asian isolates contained 7∼9 high-specificity N-glycosylation sites in their S proteins. When combining all high-specificity N-glycosylation sites determined and deleting the duplicates, eleven sites appeared in most isolates, including 57NSTW60, 112NATA115, 127NKTL130, 212NVTS215, 320NDTS323, 347NSSD350, 510NITV513, 552NVTN555, 777NISI780, 1245NKTL1248, and 1257NRTG1260 (Table 3). Among these sites, 212NVTS215, 777NISI780, and 1245NKTL1248 were conserved in almost all of the 208 Asian strains. However, the N-glycosylation sites at sites 57–60 of some strains were “NSSS” rather than “NSTW”. This is because the S proteins of these strains had the 59QGVN62 deletion compared to the S proteins of the other strains. Due to the amino acids changes, the N-glycosylation sites at sites 347–350 in some strains also changed from “NSSD” to “NSSN” or “NSTN”. The 510NITV513 and 552NVTN555 N-glycosylation sites were missing in S proteins of the isolates from South Korea and Vietnam. Due to the large deletion of amino acids at sites 23–216, the 57NSTW60, 112NATA115, 127NKTL130, and 212NVTS215 N-glycosylation sites were missing in S proteins of two Japanese strains NIG-2/JPN/2014 and KMM-1/JPN/2014, and seven N-glycosylation sites, including 130NDTS133, 157NSSN160, 494NVTS497, 549NCTE552, 587NISI590, 1055NKTL1058, and 1067NRTG1070 were retained in these two strains.

Analysis of the ORF3-E-M-N proteins

Unlike the S protein, ORF3 is a conserved protein among the PEDV isolates (Wang et al., 2016b). However, several PEDVs were found to have characteristic amino acid mutations in the ORF3 (Fig. 4). Compared to ORF3 proteins of many G1 strains and all GII strains, nine GI-a PEDVs (ZJUG12013, 85-7, 85-7-mutant 1, 85-7-mutant 2, 85-7-mutant 3, 85-7-mutant 4, 85-7-mutant 5, 85-7-A40, and 85-7-C40) had a continuous deletion of 70 amino acids at their N-terminal (positions 1–70) of ORF3; two GI-a strains (85-7-mutant 2 and 85-7-mutant 4) had a continuous deletion of 47 amino acids at their C-terminal (positions 178–224) of ORF3; while nine GI-b strains (JS2008, AH-M, SD-M, SQ2014, SC1402, HLJBY, PEDV-SX, JSLS-12015 and JS-22015) had a continuous deletion of 133 amino acids at their C-terminal (positions 92-224) of ORF3 (Fig. 4). Compared to ORF3 proteins of many GI and GII strains, ORF3 protein of a GII-a strain (CHSXYL2016) had a continuous deletion of 14 amino acids at positions 211-224; ORF3 protein of another GII-a strain (NW17) had a continuous deletion of 6 amino acids (DLYLAI) at positions 168-173; while ORF3 proteins of six GII-b strains (YN15, YN30, YN60, YN90, YN144, YN200) had a continuous deletion of 79 amino acids at their C-terminal (positions 146-224) (Fig. 4). Compared to ORF3 proteins of the GII-a isolates, more than half of the GII-b isolates had amino acid changes at positions 25 (L →S), 70 (I →V), 80 (V →F), 107 (C →F), 168 (D →N), and 182 (Q →H) (Fig. 4). Interestingly, these amino acid changes were also found in the ORF3 proteins of many GII-c isolates after 2016.

Sequence comparisons revealed that there were no common INDELs or mutations in E proteins of one subgroup of GII strains compared to E proteins of other subgroups of GII strains (Fig. 5A). M proteins of most GII-a strains had a glutamine (Q) at site 13; however, all GII-b strains isolated between 2011 and 2012 had a glutamic acid (E) at the same position in their M proteins, and this amino acid change (Q →E) occurred frequently in M proteins of GII-b since 2013 (Fig. 5B). A similar phenomenon was also observed in the M proteins of the GII-c strains, as most of the GII-c strains isolated before 2016 had a glutamine (Q) at position 13 in their M proteins, but a Q →E change at position 13 was seen in the M proteins of more frequently in strains isolated after 2016. In addition, amino acid changes at positions 192 (G →S) and 214 (S →A) appeared simultaneously in M proteins of some GII-b and GII-c strains. Similarly, in N proteins, amino acid changes at positions 216 (M →V) and 241 (R →K) appeared simultaneously in many GII strains (Fig. 5C).