Introduction

Porcine epidemic diarrhea virus (PEDV) has been recognized in veterinary medicine since the early 1970s in Europe and has subsequently been detected in many other swine-breeding areas throughout the world [27]. The virus is the causative agent of porcine epidemic diarrhea (PED), an acute and highly contagious enteric disease characterized by vomiting, watery diarrhea, dehydration, and high mortality in suckling piglets [8, 24]. Swine is the only host capable of a productive infection and serves as a reservoir of the virus [33]. Outbreaks of PED have been reported in many swine-raising countries in Europe and Asia, despite the use of vaccines [5, 6, 22]. Since the presence of PEDV in China was first confirmed in 1984, several PEDV strains have been isolated from some provinces of China. The use of bivalent inactivated or attenuated vaccines for transmissible gastroenteritis (TGE) and PED in China has brought a substantial reduction in prevalence of the disease [9]. However, a notable increase in PED cases has been observed on many swine farms in China since December 2010. The morbidity rate caused by PED ranged from 90% to 100%, with 70-100% mortality among neonatal piglets on affected swine farms, making it one of the most devastating enteric diseases of swine, which has resulted in huge economic losses to the pig-farming industry [12, 17].

PEDV, a member of the genus Alphacoronavirus (family Coronaviridae, order Nidovirales), is an enveloped, single-stranded, positive-sense RNA virus with a genome of approximately 28 kb in length [1, 25]. The PEDV genome contains seven open reading frames (ORFs) in the following order: ORF1a, ORF1b, the spike (S) glycoprotein, ORF3 (a hypothetical protein gene), envelope (E), membrane (M) and nucleocapsid (N) [21, 24]. Of the viral proteins, the S glycoprotein, containing 1383 amino acids (aa), with a molecular weight of about 150 kDa, is considered the most antigenic protein. Previous studies have shown that the S glycoprotein plays an important role in promoting the fusion of viral and cellular membranes, virulence, induction of neutralizing antibodies, and attachment of viral particles to receptors on host cells [16, 30]. This region is frequently used to determine the genetic relationships among different PEDV strains and the epidemiological status of PEDV strains circulating on swine farms. ORF3 is the only accessory gene in the PEDV genome. It has been demonstrated that the virulence of PEDV can be reduced by altering the ORF3 gene through cell culture adaptation, and it has been suggested to be a determinant of the virulence of the virus [26]. Analysis of the ORF3 gene can be used to discriminate between vaccine and wild-type PEDVs, providing a valuable tool for molecular epidemiological surveys of PEDV infections [19].

Although there have been several reports on the prevalence and genetic variation of PEDV in central China before 2016 [7, 15, 18, 23, 28, 32, 35], its prevalence in recent years remains unclear. For further molecular epidemiology and phylogenetic analysis of recent PEDV strains, intestinal tissues and fecal samples were obtained from diseased piglets during outbreaks of diarrhea between 2015 and 2019 on farms in Henan and Shanxi provinces of China where swine had been immunized with attenuated PEDV (CV777), and PEDV was screened by reverse transcription polymerase chain reaction (RT-PCR). The nucleotide and amino acid sequences of the S and ORF3 genes of 32 PEDV-positive samples were analyzed and characterized by comparing them with PEDV reference strains from Asia, the United States, South Korea, and Europe [10, 14, 15, 34, 36]. This study provides up-to-date information on PEDVs circulating in Henan and Shanxi provinces of China between 2015 and 2019. This information will be helpful for developing new strategies for prevention and control of PEDV.

Materials and methods

From January 2015 to December 2019, a total of 135 clinical samples (intestinal tissues and feces) were collected from diseased piglets during outbreaks of diarrhea on farms in Henan and Shanxi provinces of China where swine were immunized with attenuated PEDV (CV777). The homogenates of these clinical samples were diluted 1:5 in phosphate-buffered saline (PBS) and centrifuged at 8000 g for 5 min. Viral RNA was extracted from samples using TRIzol Reagent (Takara, Dalian, China) according to the manufacturer’s instructions, and nucleic acids were eluted in 30 μL of RNase-free water. Reverse transcription was carried out according to the manufacturer’s instructions (Vazyme, Nanjing, China).

Complete genome sequences of PEDV strains were downloaded from the GenBank database and aligned using the MegAlign program of DNAStar software (version 7.1, DNASTAR Inc., Madison, WI. USA). A pair of primers (M-F, 5’-CCTTATGGCTTGCATCACTCT-3’; M-R, 5’-CCCAAGCACTTTCTCACTATC-3’) was designed based on a highly conserved region within the M gene using Primer Premier software (version 5.0) (Primer 5.0) to detect PEDV with an amplicon of 419 bp. The reaction mixture consisted of 2 μL of cDNA, 12.5 μL of 2 × Ftaq PCR MasterMix (ZOMANBIO, Beijing, China), 1 μL of primer M-F (25 μM), 1 μL of primer M-R (25 μM), and 8.5 μL of RNase-free water in a total volume of 25 μL. The amplification parameters were as follows: 94 °C for 5 min, followed by 35 cycles of 94 °C for 30 s, 56 °C for 30 s, and 72 °C for 45 s, and a final elongation step for 10 min at 72 °C.

Two pairs of primers S1-F/S1-R (F1, 5’-GAAGGTAAGTTGCTAGTGCGTAA-3’; R1, 5’-AGGTAGCCAATACTGCCAGATTT-3’) and S2-F/S2-R (F2, 5’-GTGGC CTGTGTTGGTGTATAG-3’; R2, 5’-GGTGCCTCAAAGAAGACGCTT-3’) were designed to amplify two overlapping cDNA fragments spanning the entire S gene, and the full S gene was amplified by PCR from the PEDV-positive samples using the primer sets for the S gene. The full ORF3 gene was also amplified using the primers ORF3-F/ORF3-R (F3, 5’-GGCGTCCTAGACTTCAACCTT-3’; R3, 5’-GGACTGC GCTATTACACAACC-3’). PCR products were purified using a QIAquick Gel Extraction Kit (QIAGEN) according to the manufacturer’s instructions and cloned into pMD18-T vector (Takara) at 16°C overnight. The resulting plasmids were introduced into Escherichia coli DH-5α cells (Takara) by transformation according to the manufacturer’s instructions, and positive clones were visualized by β-galactosidase screening and isolated. The positive plasmids were sent to Sangon Biotech Shanghai Co., Ltd. for sequencing. All sequencing reactions were performed in duplicate.

The full S and ORF3 genes of PEDV strains were aligned using the MegAlign program of DNAStar software. Phylogenetic trees of the S and ORF3 genes were constructed by the neighbor-joining method in Molecular Evolutionary Genetics Analysis (MEGA) software (version 7.0) with a p-distance model, using 1,000 bootstrap replicates. Additionally, possible recombination events were identified using the RDP4.0 recombination detection program. The RDP, BOOTSCAN, MaxChi, Chimaera, SiScan and 3Seq programs embedded in the RDP4.0 software package were used for the analysis [3], and only those recombination events supported by more than four programs (P < 0.01) were considered [15].

Results and discussion

On the basis of the RT-PCR amplification results, 86 out of 135 clinical samples were found to be positive for PEDV, giving a PEDV-positive rate of 63.7% (86/135), which is consistent with a recent report [13]. This suggests that PEDV remains a major pathogenic agent in China. Thirty-two positive samples representing different cities or sampling time points in Henan and Shanxi provinces were selected to sequence the S and ORF3 genes of PEDV, and the sequences were subsequently submitted to the GenBank database. The origin, date, and GenBank accession numbers of these PEDV strains are summarized in Table 1.

Table 1 Information about PEDV strains used for sequence alignment and phylogenetic analysis

The S gene is one of the most variable genes in the PEDV genome and is considered a useful marker for analyzing genetic variation of PEDV strains [2]. In this study, of the 32 full S sequences examined, 26 were 4,161 nucleotides in length, encoding a protein of 1,387 aa, and the other six (SX/TY1/2017, SX/TY2/2017, HN/KF/2017, HN/HX/2019, HN/XZ/2017 and HN/JY/2018) were 4,158 nt long (1,386 aa). A multiple sequence alignment showed that the pairwise sequence identity of the 32 PEDV strains ranged from 96.2% to 99.9%, and these strains shared 92.1% to 99.7% nucleotide sequence identity with 31 reference strains. Phylogenetic analysis of the S gene demonstrated that the PEDV strains fell into two groups (Fig. 1). Group 1 (prototype PEDV) included a German strain (Br1/87), South Korean strains (DR13 Virulent and SM98), early domestic strains (LZC and CH-S), and the vaccine strain (CV777). All 32 PEDV strains in this study were clustered in group 2 (subgroups G2-a and G2-b), together with 17 Chinese strains and six foreign strains, including American strains (PC22A-P140.BI, USA/Colorado/2013 and USA/OK10240-6/2017), South Korean strains (KNU-1305 and KNU-1807) and a Mexican strain (PEDV/MEX/QRO/02/2017). These results suggest that the 32 PEDV strains might have evolved from a Chinese highly virulent PEDV strain after 2010. Phylogenetically, 27 of the 32 PEDV strains in this study and 16 reference strains were distributed within the G2-a genotype cluster and had 97.0%-99.9% nucleotide sequence identity. However, the 27 PEDV strains were in different clusters from the six strains from other countries. Notably, CH/ZMDZY/11, a virulent PEDV strain isolated from Henan province in 2011 [32], was clustered with the six foreign circulating strains mentioned above and formed a relatively independent sub-branch, sharing 96.6% to 98.9% nucleotide sequence identity with the 32 PEDV strains. These results suggested that the 32 PEDV strains have undergone rapid variation and evolution between 2015 and 2019. For subgroup G2-b, HN/XZ/2017, HN/KF/2017, HN/JY/2018, SX/TY2/2017 and SX/TY1/2017 were grouped in one cluster with seven Chinese reference strains (XY2013, CH/YNKM-8/2013, ZJCZ4, CHSD2014, AJ1102, CH/GDGZ/2012 and PEDV JS-A), with nucleotide sequence identity ranging from 97.1% to 99.1%. The remaining two American strains (OH851 and USA/IL20697/2014) belonged to G2-c (INDEL strain). These results show that the 32 strains differ genetically from the reference strains, as indicated by their distant location in the phylogenetic tree.

Fig. 1
figure 1

Phylogenetic analysis based on the S gene sequence of 32 PEDV strains obtained in the present study together with 31 PEDV reference strains from GenBank. The tree was generated by the neighbor-joining method using MEGA7 software. Bootstrap values are indicated for each node, based on 1000 replicates. The positions of the 32 PEDV strains from the present study are indicated by solid black circles. The recombinant strain SX/TY2/2017 is indicated in bold. GenBank accession numbers are shown in parentheses

On the basis of S sequences and comparison with reference strain CV777, S gene nucleotide sequence identity values ranged from 93.2% to 93.7% between the 32 PEDV strains and CV777. The 32 strains displayed two common aa deletions (163DI164) and, frequently, aa insertions (59QGVN62) in the S protein relative to the reference PEDV strains. It has been demonstrated that four major neutralizing epitopes (aa 499–638, 748–755, 764–771 and 1,368–1,374) are present on the surface of the PEDV S protein [23, 29]. In this study, amino acid sequence alignments of the S region with the CV777 strain indicated that these neutralizing epitope regions were highly variable, with major variable sites (Table 2). This suggests that these aa changes might have been involved in immune escape, reducing the effectiveness of the CV777-based vaccine and ultimately leading to outbreaks of severe diarrhea on Chinese swine farms. Elsewhere, it has been reported that a novel B-cell epitope, SE16, is the immunodominant region (722SSTFNSTREL731) of the PEDV S protein[20]. In our study, two aa mutations (F→L at 725; S→G at 727) were found in the strain CH/XX/2015 (KU237231), and experimental studies are needed to determine whether these mutations affect vaccines and diagnostic reagents. One aspartic acid (N) deletion at aa position 1,198 and five aa substitutions (D→H at aa 1,199, Y→F at aa 1,212, G→S at aa 1,220, E→D at aa 1,242, and S→P at aa 1,270) in the S protein were found in six strains (SX/TY1/2017, SX/TY2/2017, HN/KF/2017, HN/HX/2019, HN/XZ/2017 and HN/JY/2018) when compared with the other 26 strains in this study. The question of whether these mutations affect the antigenicity and pathogenicity of PEDV variants warrants further investigation.

Table 2 Amino acid mutations in neutralizing epitopes of the putative S protein of 32 PEDV strains from this study compared with CV777

Phylogenetic analysis of the ORF3 gene showed that the PEDVs could be divided into two main groups (Fig. 2). All of the PEDV strains in this study belonged to the G2-1 and G2-2, together with most Chinese strains and nine foreign strains, including five American strains (PC22A-P140.BI, OH851, USA/Colorado/2013, USA/OK10240-6/2017, and USA/IL20697/2014), three South Korean strains (DR13 Virulent, KNU-1305, and KNU-1807) and one Mexican strain (PEDV/MEX/QRO/02/2017). Group 1 included a Chinese strain (LZC), a South Korean strain (SM98), a German strain (Br1/87), and the CV777 vaccine strain. The phylogenetic tree showed that the 32 strains in this study had a close relationship to strains isolated in China, South Korea, and America but differed genetically from the other reference strains (SM98, Br1/87, LZC) and the CV777 vaccine strain.

Fig. 2
figure 2

Phylogenetic analysis based on the ORF3 gene sequence of the 32 PEDV strains obtained in the present study together with 31 PEDV reference strains from GenBank. The tree was generated by the neighbor-joining method using MEGA7 software. Bootstrap values are indicated for each node, based on 1000 replicates. The positions of the 32 PEDV strains from the present study are indicated by solid black circles. GenBank accession numbers are shown in parentheses

ORF3 is the only accessory gene of PEDV and is associated with virus production and virulence [31]. It can be used as a marker for identifying attenuated and virulent strains [26]. A previous study suggested that the ORF3 gene may be involved in cell tropism and the pathogenicity of the virus [12]. Park et al. [26] reported that a cell-culture-adapted PEDV with a smaller ORF3 gene exhibited reduced virulence compared with PEDV field strains. The PEDV strains in this study had intact ORF3 genes that encoded a 225-aa protein. A multiple sequence alignment illustrated that the ORF3 sequences of the 32 PEDV strains were 94.7% to 100% identical to each other and were 94.4% to 100% identical to those of the reference strains (except SM98). It has been reported that the ORF3 protein of wild-type PEDV CV777 contains helical transmembrane domains (TMDs) at residues 40–63 (TM1), 74–97 (TM2), 118–136 (TM3), and 151–172 (TM4), forming a tetrameric assembly that is important for potassium channel activity is associated with the virulence of PEDV [31]. In comparison to the CV777 reference strain, the following aa changes were observed: V→I at 54, V→I at 79, L→F at 92, and N→S at 166, all of which have been reported previously [35]. These mutations might have an impact on the virulence and adaptability of PEDV, and this needs to be explored further in animal experiments. Moreover, in 25 of the PEDV strains, an additional mutation was identified at position 80 (F→V), but this mutation was not found in seven strains (HN/JY/2018, HN/XC1/2018, HN/XC2/2018, SX/TY1/2017, SX/TY2/2017, HN/XX/2016 and SX/WS1/2018). Interestingly, one mutation at position 168 (D→N), which was found in HN/JY/2018, HN/XC1/2018, HN/XC2/2018, SX/TY1/2017, SX/TY2/2017, HN/XX/2016, and SX/WS1/2018, was also observed in Chinese PEDV strains (ZJCZ4, AJ1102, CH/GDGZ/2012 and PEDV JS-A). These strains were distributed within the G2-2 genotype cluster and had a close relationship to one another (Fig. 2). These results show that the ORF3 gene sequences in this study have a tendency toward variation.

Comparison of the S and ORF3 gene sequences of 32 selected PEDV field strains and the prototype CV777 strain showed that the field strains differed significantly from the early Chinese CV777 vaccine strain. The amount of genetic divergence within the 32 PEDV strains displayed a complex trend from 2015 to 2019, which might be related to the transfer of pigs and a lack of biosafety awareness on Chinese farms. Phylogenetic analysis indicated that the identified PEDV strains were non-S-INDEL strains and exhibited genetic diversity. The 32 PEDV strains in this study belonged to group 2 (pandemic variant strains) and evolved from the 2010 outbreak. Based on phylogenetic analysis of the S gene (classical grouping), the 32 PEDV strains could be classified into two subgroups (G2-a and G2-b), and 27 of 32 (84.38%) PEDV strains from Dengfeng (n = 3), Kaifeng (n = 2), Luyi (n = 1), Neihuang (n = 1), Pingdingshan (n = 1), Shangshui (n = 1), Xunxian (n = 4), Yuanyang (n = 2), Yuzhou (n = 2), Xinxiang (n = 1), Puyang (n = 1), Wenshui (n = 1), Xuchang (n = 3), Hebi (n = 2), Huixian (n = 1), and Taiyuan (n = 1) reported here were randomly distributed in the G2-a subgroup. In contrast, the other five PEDV strains (15.62%) were assigned to the G2-b subgroup. These observations suggest that G2-a subgroup strains were the dominant pandemic variant strains circulating in Henan and Shanxi provinces, which is consistent with recent reports [13]. These data will be helpful for research on PEDV vaccines and control.

Recombination analysis revealed one putative recombination event existed in strain SX/TY2/2017. The putative breakpoints were located at nt positions 1 and 1,116 of the S gene. The strains CH/GDGZ/2012 and CH/YZ1/2015 were identified as parental strains, which were confirmed using five programs (RDP, Av. P-Val = 4.894 × 10−5; MaxChi, Av. P-Val = 3.68 × 10−5; Chimara, Av. P-Val = 5.349 × 10−6; SiScan, Av. P-Val = 7.463 × 10−5; 3Seq, Av. P-Val = 4.21 × 10−4). Sequence comparisons of the recombinant and non-recombinant regions showed that the sequence of the potential major parental strain CH/GDGZ/2012 was 98.8% identical and that of the minor parental strain CH/YZ1/2015 was 98.9% identical (Fig. 3). In the phylogenetic tree based on the S gene, the recombinant strain (SX/TY2/2017) clustered closely with the major parental strain CH/GDGZ/2012, situated in group G2-b, whilst the putative minor parental strain CH/YZ1/2015 branched outside the group G2-b cluster. Therefore, the significantly discrepant topologies of the phylogenetic tree support the possibility of further recombination events in these strains. Indeed, previous studies have provided evidence that recombination among different strains is a major evolutionary factor that could allow the emergence of new strains with altered virulence and immunogenicity [4, 11]. Notably, the same subtype and recombination pattern in the S gene was found on different swine farms in different regions of Henan and Shanxi provinces, suggesting that the new subtypes of PEDV strains have spread and circulated in these provinces. Hence, the novel recombination event in the S gene of PEDV strains in Henan and Shanxi provinces further contributed to our understanding of the evolution and genetic diversity of PEDV. Recombination plays a pivotal role in the evolution of PEDV, and it allows the emergence of new variant strains with altered virulence, increasing genetic complexity, and requiring new strategies for vaccine development. Recently, pigs have been frequently transported and traded across China, especially in Henan and Shanxi provinces, which has resulted in the rapid transmission and accelerated mutation of PEDV in pigs in these provinces. Pig farmers should therefore strengthen their biosafety awareness and management, limit live pig transshipment, and attempt to reduce transmission of PEDV.

Fig. 3
figure 3

Identification of possible recombination events in SX/TY2/2017 using RDP4 (using six recombination detection programs: RDP, BOOTSCAN, MaxChi, Chimaera, SiScan, and 3Seq). One putative recombinant region (P < 0.01) was located at nt 1 to 1,116 in the S gene. The analysis was performed with an RDP distance model and a window size of 20

In summary, our results demonstrated that PEDV is one of the pathogens contributing to outbreaks of clinical diarrhea on swine farms, and the PEDV strains identified in our study are pandemic variants exhibiting genetic diversity. Recombination, insertions and deletions in the S and ORF3 genes might contribute to the genetic diversity of PEDV and might also be associated with vaccination failure and the emergence of new PEDV variants, making disease prevention more complicated. Collectively, our findings emphasize that the PEDV vaccines based on recent variants urgently need to be developed for the epidemic response to the current PED outbreak in China and that the pathogenesis and immunogenicity of PEDV field strains require further exploration in order to better prevent and control the epidemic of PEDV.