Emergence and evolution of highly pathogenic porcine epidemic diarrhea virus by natural recombination of a low pathogenic vaccine isolate and a highly pathogenic strain in the spike gene.

Outbreaks of a new variant of porcine epidemic diarrhea virus (PEDV) at the end of 2010 have raised interest in the mutation and recombination of PEDV. A PEDV strain (CN/Liaoning25/2018) isolated from a clinical outbreak of piglet diarrhea contained a 49-bp deletion in the ORF3 gene. This deletion is considered a genetic characteristic of low pathogenic attenuated vaccine strains. However, CN/Liaoning25/2018 was highly pathogenic. Complete genome sequencing, identity analysis, phylogenetic tree construction, and recombination analysis showed that this virus was a recombinant strain containing the Spike (S) gene from the highly pathogenic CN/GDZQ/2014 strain and the remaining genomic regions from the low pathogenic vaccine isolate SQ2014. Histopathology and immunohistochemistry results confirmed that this strain was highly pathogenic and indicated that intestinal epithelial cell vacuolation was positively correlated with the intensity and density of PEDV antigens. A new natural recombination model for PEDV was identified. Our results suggest that new highly pathogenic recombinant strains in the field may be generated by recombination between low pathogenic attenuated live PEDV vaccines and pathogenic circulating PEDV strains. Our findings also highlight that the 49-bp deletion of the ORF3 gene in low pathogenic attenuated vaccine strains will no longer be a reliable standard to differentiate the classical vaccine attenuated from the field strains.

1. Introduction

Porcine epidemic diarrhea virus (PEDV) is a single-stranded, positive-sense RNA virus belonging to the order Nidovirales, the family Coronaviridae, and the genus Alphacoronavirus. PEDV can cause severe diarrhea in infected piglets (Wang, Byrum, and Zhang 2014). Although PEDV was first identified over 30 years ago (Pensaert and de Bouck 1978), it was not until the end of 2010 that a new PEDV strain has appeared in China and the disease has begun to emerge on a large-scale worldwide (Lin et al. 2016). This variant PEDV strain causes 80–100 per cent morbidity and 50–90 per cent mortality in suckling piglets, leading to considerable economic loss to the global pig industry (Carvajal et al. 2015).

The genome of PEDV is approximately 28 kb and is composed of seven open reading frames (ORFs) in addition to a 5ʹ untranslated region (UTR) and a 3ʹ UTR (Kocherhans et al. 2001). Two-thirds of the 5′ region of the genome contain two overlapping ORFs, ORF1a and ORF1b, which encode nonstructural proteins that direct genome replication, transcription, and viral polyprotein processing. The remainder of the PEDV genome contains five ORFs encoding four structural proteins and an accessory protein in the following order: spike (S), ORF3, envelope (E), membrane (M), and nucleocapsid (N) proteins (Duarte et al. 1993). The PEDV S protein is a type 1 transmembrane envelope glycoprotein with S1 (amino acid 1–735) and S2 (736–the last amino acid) domains. S1 mediates the attachment of virus particles to the cell surface (via receptor interactions) and stimulates the induction of neutralizing antibodies in the natural host, and S2 mediates the fusion of viruses to host cells (Chang et al. 2002; Sun et al. 2008; Li 2015; Kim et al. 2018). Mutations or deletion in the S protein may change the pathogenicity and tissue tropism of PEDV (Marthaler et al. 2014; Zhang et al. 2015; Suzuki et al. 2018).

Based on a phylogenetic analysis, the complete PEDV genomes have been divided into two separate groups: genogroup I (GI; classical) and genogroup II (GII; variant or field epidemic). GI has evolved into two subgroups (GI-a and GI-b) and GII has evolved into three subgroups (GII-a, GII-b, and GII-c) (Guo et al. 2019). GII-c is also known as the S-INDEL subgroup.

Recombination is an important evolutionary factor in many RNA viruses (Simon-Loriere and Holmes 2011). PEDV recombination primarily occurs in seven main regions of the concatenated full genome, including nsp2, nsp3, nsp14–16, the S1 domain, and the nucleocapsid genes (Jarvis et al. 2016). There are two patterns of recombination within and between different virus subgroups/species. In one pattern, PEDV recombines with other coronaviruses, such as the chimeric virus SeCoV/GER/L00930/2012 the Italy/213306/2009 strain or PEDV/Belgorod/dom/2008 (Nefedeva, Titov, and Malogolovkin 2019), a recombinant form of PEDV and porcine transmissible gastroenteritis virus (TGEV). The other recombination pattern involves the recombination of PEDV with different PEDV strains, including recombination of different pathogenic strains from the same or different subgroups. Examples include the XM1-2 strain, derived from the recombination of pathogenic GYJ130330 and HNQX-3 strains of same subgroup GII-a (Chen et al. 2017), and the whole strains of subgroup GII-c derived from the recombination of the pathogenic strains of different subgroups (GI-a and GII-a) (Guo et al. 2019). Recombination has enriched the genetic diversity of PEDV and continues to pose challenges to the prevention and treatment of PEDV. Changes in the antigenic epitopes can lead to the failure of established vaccine responses (Servin-Blanco et al. 2016). Here we report the emergence of a new natural recombination event, in which the CH/Liaoning25/2018 strain was formed by the replacement of most of the S gene from the low pathogenic vaccine isolate, strain SQ2014, by that from the highly pathogenic CN/GDZQ/2014 strain through natural recombination. CH/Liaoning25/2018 was characterized by pathogenic features, antigen distribution, genetic characterization, phylogeny, and evidence of recombination.

2. Materials and methods

2.1. Rt-PCR diagnosis and virus isolation

In February 2019, outbreaks of fatal swine diarrhea appeared at a farm in the Liaoning Province of China, and spread to a nearby farm. PEDV, porcine deltacoronavirus (PdCV), TGEV, porcine respiratory coronavirus (PRCV), porcine rotavirus (PORV), African swine fever (ASF), classical swine fever (CSF), and swine acute diarrhea syndrome coronavirus (SADS-CoV) were detected by RT-PCR as described previously (Kim et al. 2000; Wang, Costantini, and Saif 2007; Cheng et al. 2008; Amimo, Vlasova, and Saif 2013; Oka et al. 2014; Zhai et al. 2016; Luo et al. 2017; Zhou et al. 2018).

Intestinal samples positive for PEDV were homogenized in DMEM (Gibco, Beijing, China), centrifuged at 4,000 rpm/min, filtered through a 0.22-μm filter (Millipore, MA) and used as inoculum. Vero cells (CCL-81, ATCC) were plated into 10 cm dishes one day before infection. Vero cells were inoculated with 1 ml of filtered inoculum for 1 h at 37 °C. After virus adsorption, cells were maintained in DMEM containing 8 μg/ml trypsin (Gibco, Beijing, China) and antibiotics at 37 °C with 5 per cent CO2. When obvious cytopathic effects were apparent, including the characteristic fused-syncytium cell formation, the viral supernatant was collected and used for RNA extraction for subsequent cloning and sequencing. The 49-bp deletion of the ORF3 gene was assessed for the differential diagnosis between vaccine and field strains for PEDV positive samples (Zhu et al. 2016; He et al. 2019). The primers used for the differential diagnosis and S gene cloning are shown in Supplementary Table S1.

2.2 Histopathology and immunohistochemistry (IHC)

Gastrointestinal tissue sections were formalin-fixed and paraffin-embedded by the standard methods, followed by hematoxylin and eosin (H & E) staining. Detection of the PEDV N antigen was performed using the Polymer Horseradish Peroxidase Immunohistochemical Staining Kit (Beijing Zhongshan Jinqiao Biotechnology Co., Ltd. Beijing, China) as described previously (Wang et al. 2013). Mouse anti-PEDV nucleocapsid (N) monoclonal antibody (2G3 mAb), a gift from Prof. Li Feng (Harbin Veterinary Research Institute, China), was used as the primary antibody (1:200 in phosphate-buffered saline—PBS). A duplicate of each section was incubated in PBS only, as a negative control. Pathological and IHC diagnoses of the tissue sections were performed and scored in a blinded manner by a veterinary pathologist. The intensity of PEDV positive cell staining was calculated using the Brightfield Immunohistochemistry Image Analysis Software (Leica, Wetzlar, Germany). Three to five microscopic fields were randomly selected from each section.

2.3. Genome sequences

Virion RNAs were extracted from PEDV-infected cell supernatants using the TRIzol LS reagent (Invitrogen, Shanghai, China). cDNAs were then synthesized using the Revert Aid H Minus First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, Shanghai, China). The extreme 5′ and 3′ termini of the CN/Liaoning25/2018 strain were cloned and sequenced using the SMARTer® RACE 5′/3′ Kit (Takara Bio, USA). Primers were designed based on the AF353511 and KM242131 accession numbers for the complete genome determination of PEDV (Supplementary Table S1). The RT-PCR products were purified in 1.5 per cent agarose gels and cloned into the pEASY®-Blunt vector or pEASYBlunt Simple vector (TransGen Biotech Co., LTD, China) with three to five clones chosen. Clones were sent to GENEWIZ, Inc (Jiangsu, China) or TIANYIHUIYUAN Biotech Co., Ltd (Beijing, China) for sequencing via the ABI 3730XL apparatus. All clone sequences were assembled with the Seqman package (DNAStar, Madison, WI). The vector sequences were removed and their 5′ and 3′ ends were trimmed based on a SeqMan’s trace quality evaluation. After all sequences were assembled, the conflicting bases were checked and confirmed according to chromatograms of no less than three cloned sequences. The complete genome sequence was obtained. Subsequent analysis showed that the strain was a recombinant strain, and the breakpoints of recombination coincided with the primer scheme. In order to avoid recombination caused by virus isolation after culture, or the amplification or assembly steps, three clinical intestinal tissues were selected, and the region of breakpoints (2,000–21,100 and 24,200–25,200, approximately 500 bp before and after breakpoints) was cloned into the pUC19 vector, sequenced and re-confirmed. The complete genome sequence (after validation) was submitted to the GenBank database.

2.4. Recombination and phylogenetic analyses

In order to more accurately determine the possible recombination events, all of the complete genome sequences of PEDV (n = 480) available up to 30 July 2019 were collected from the GenBank (Supplementary Table S2). Patents, artificial mutants, and repeated PEDV sequences were removed from the data set. To screen for recombination, a suite of recombination detection methods were used, including RDP, GENECONV, Chimaera, MaxChi, and 3Seq in RDP version 4.95 (Martin et al. 2015), followed by secondary scanning using BootScan and SiScan in RDP, with a cutoff of P = 0.01. SimPlot v.3.5.1 was used to confirm recombination events; nucleotide (nt) sequence similarity was also assessed (Lole et al. 1999).

To further investigate the recombination and the evolutionary relationship between CN/Liaoning25/2018 and the parental PEDV strains, phylogenetic analysis was performed. The complete genomes, nonrecombination regions, and recombination regions of the CN/Liaoning25/2018 strain were aligned with those of 65 representative PEDV genomes retrieved from the GenBank using MAFFT v. 7.037 (Supplementary Table S3). The bat coronavirus BtCoV/512/2005 sequence (Katoh and Standley, 2013) (GenBank accession number DQ648858) was used as an outgroup (Tang et al. 2006). A maximum-likelihood (ML) phylogenetic tree was constructed using IQ-TREE v.1.6.6 (Nguyen et al. 2015) following 5,000 bootstrap replicates. A separate Bayesian inference (BI) phylogenetic analysis was constructed using MrBayes v. 3.26 and the same best-fitting evolutionary model suggested by IQ-TREE (Katoh and Standley 2013). The robustness of the phylogenetic topology was tested using the Kishino–Hasegawa (KH) test and the Shimodaira–Hasegawa (SH) test from IQ-TREE (Nguyen et al. 2015).

2.5. Statistical analyses

The cell staining intensity data were analyzed using the Prism software, version 8.01 (GraphPad, La Jolla, CA). The Kruskal–Wallis test with Dunn’s post hoc analysis was applied, with a level of significance of P < 0.05 (Kruskal and Wallis 1952). The quantitative data are presented as mean ± SEM.

3. Results

3.1. Clinical symptoms and diagnosis

The main clinical symptoms included gray, yellow, or white watery diarrhea and vomiting (Fig. 1A). Piglets <3 days-of-age were more susceptible, with 90 per cent morbidity and 80–90 per cent mortality. The disease course was approximately one week. The small intestine, mainly the jejunum and ileum, was dilated. The wall of the small intestine was relaxed, lacked elasticity, and was thin and transparent. The intestinal contents were thin, yellow, and foamy. Mesenteric hyperemia as well as edema of the mesenteric lymph nodes was present (Fig. 1B and C).

Figure 1.

Clinical and necropsy findings in infected piglets. (A) Clinical presentation of infected piglets. (B) and (C) Typical pathological changes of the intestine.

Total RNA was extracted from intestine tissues of piglets with diarrhea and RT-PCR was used to determine the presence of various viral genes. Samples were only positive for PEDV. They were negative for PdCV, TGEV, PRCV, PORV, ASF, CSF, and SADS-CoV (data not shown). To further confirm the disease was caused by PEDV, sections of samples collected from the stomach, the small intestine, colon, cecum, mesenteric lymph nodes, tonsils, lungs, heart, liver, spleen, and kidneys were examined by H & E staining and immunochemistry. The histological examinations revealed lesions in the intestinal tissue, including atrophy of the intestinal villi in the duodenum, jejunum, ileum, and colon, which was accompanied by varying degrees of intestinal columnar epithelial cells’ vacuolization. The number of vacuoles was the highest in the jejunum and the lowest in the colon (Fig. 2A, D, G, and J, Table 1). No obvious pathological changes were observed in other tissues (data not shown).

Figure 2.

Histopathological features, PEDV antigen distribution, and PCR results. (A), (D), (G), and (J) Intestinal histopathological features of the affected piglets in H & E stained sections. Vacuolated epithelial cells are indicated by arrows. (B), (C), (E), (F), (H), (I), (K), and (L) show IHC staining to detect PEDV antigen in the different intestinal tracts of one infected piglet. The approximate locations of the first and third column tissues are indicated by the black squares in the second column tissue. (M) IHC analysis showing the staining density in different tissues, **P < 0.01, *P < 0.05. (N) PCR results: lane 1, marker; lane 2, field strain positive control; lane 3, vaccine strain positive control; lane 4, ORF3 gene detection of sample strain; lane 5, S gene detection.

Table 1.

Histopathological changes and IHC staining of intestinal columnar epithelial cells from infected piglets.

	Histopathology	IHC
	Vacuolationa	Staining intensityb	Staining densityc
Duodenum	−	+	+
Jejunum	+++	+++	+++
Ileum	++	++	++
Colon	+	+	++

−, no vacuolation of epithelial cells; +, 1%–25% of epithelial cells showed vacuolation; ++, 25%–50% of epithelial cells showed vacuolation; +++, 50%–75% of epithelial cells showed vacuolation; and ++++, >75% of epithelial cells showed vacuolation.

+, light and mild stain, ‘Minimal’; ++, light-medium stain, ‘mild’; +++, moderate stain, ‘Moderate’; and ++++, marked/strong stain, ‘Marked’.

+, 1%–25% of epithelial cells stained; ++, 25%–50% of epithelial cells stained; +++, 50%–75% of epithelial cells stained; and ++++, >75% of epithelial cells stained.

IHC for PEDV-N confirmed the presence of the PEDV antigen in the cytoplasm of villi enterocytes. No PEDV antigen was detected in other tissues. The staining intensity and density of positive cells from different tissues were manually scored and showed higher levels in the jejunum and lower levels in the duodenum (Fig. 2B, C, E, F, H, I, K, L, M, Table 1). To reduce subjective bias, an IHC analysis software was used to statistically analyze the staining density in different tissues. This analysis confirmed the previous manual scoring findings (Fig. 2M). Moreover, the cell density staining of the duodenum was significantly different from that of the jejunum (P < 0.01) and ileum (P < 0.05), and there was no significant difference in comparison with other tissues. More vacuolated epithelial cells and stronger staining intensity and density were observed in different tissues. These results showed that there was a higher intensity and density of stained cells when the tissues contained more vacuolated cells. On the contrary, there was less staining when the tissues contained fewer vacuolated cells. The number of vacuolated cells was positively correlated with the intensity and the density of stained cells, indicating that PEDV was the main factor causing the vacuolization of intestinal epithelial cells.

Since an attenuated live vaccine (CV777) was used for immunizations before the PEDV outbreak, we investigated whether PEDV was caused by a wild strain or a vaccine strain. The CV777 vaccine strain was characterized by a 49-bp deletion in the ORF3 gene, which is usually used in the differential diagnosis between wild and vaccine strains (Zhu et al. 2016). RT-PCR revealed that a 270-bp band (49-bp deletion) of the ORF3 gene was detected (different from the observed for field strains—319 bp) (Fig. 2N). The S protein is an important determinant of PEDV virulence (Zhang et al. 2015; Suzuki et al. 2018). Therefore, the S gene of the isolated strain was cloned, sequenced, and compared with representative strains of the five PEDV subgroups (GI-a: virulent CV777, AF353511; GI-b: attenuated CV777, KT323979; GII-a: AH2012, KC210145; GII-b: AJ1102, JX188454; GII-c: OH851, KJ399978) (Guo et al. 2019). The S gene of this isolated strain shared only 93.2 per cent and 93.4 per cent identity with the CV777 vaccine and field strains, respectively, but shared the highest identity (98.7%) with the AH2012 strain (Table 2). In order to fully understand the genetic characteristics of this strain, its complete genomic sequence was cloned and sequenced.

Table 2.

Nucleotide identity of CN/Liaoning25/2018 compared with five PEDV reference strains.

Subgroup	Representative Strains	GenBank	Identity
			Complete genome	ORF1ab gene	S gene	ORF3_E_M_N gene
GI-a	Virulent CV777	AF353511	97.2	98.1	93.4	96.5
GI-b	Attenuated CV777	KT323979	98.8 a	99.8	93.2	99.6
GII-a	AH2012	KC210145	97.7	97.8	98.7	95.6
GII-b	AJ1102	JX188454	97.6	97.8	97.6	96.1
GII-c	OH851	KJ399978	97.2	97.8	95.5	95.8

Segments of the reference strains sharing the highest similarity to CN/Liaoning25/2018 are shown in bold.

3.2. Unique genetic characteristics of the emergent PEDV strain

The complete viral genome sequence (27,961 bp) of the isolated strain that was designated CN/Liaoning25/2018 was obtained and deposited in the GenBank (accession no. MK796238). A total of 38 nt differences were identified between the CN/Liaoning25/2018 strain and the representative strains of the five gene subgroups (Supplementary Table S4), which resulted in 17 amino acid substitutions. The nonstructural protein 1a (ORF1a) and M protein genes harbored four and one amino acid substitutions, respectively. The remaining 12 amino acid differences were all found in the S protein, including 11 amino acid substitutions and 1 amino acid deletion. The V601 and V705 substitutions were located in the S protein neutralization epitope regions (COE and S1D, respectively) (Fig. 3A) (Chang et al. 2002; Sun et al. 2008). Thus, the substitutions of V601 and V705 could change the epitope determinants for neutralization. We also found a 12-nt and 49-nt deletion in the ORF1a and ORF3 genes of the CN/Liaoning25/2018 strain, respectively, which were similar to those observed in the attenuated CV777 strain (Fig. 3B). This explained why the ORF3 gene band of the isolated strain (observed using RT-PCR) suggested the presence of the attenuated vaccine strain.

Figure 3.

Molecular characterization of the emergent PEDV strain. (A) Locations of the unique amino acid (aa) changes identified. nt and aa differences between the CN/Liaoning25/2018 sequences and the consensus sequences of representative strain genomes of five gene subgroups (GI-a: field CV777, GI-b: attenuated CV777, GII-a: AH2012, GII-b: AJ1102, and GII-c: OH851) with their positions are depicted. Con: PEDV consensus sequences; CN: unique sequences in the CN/Liaoning25/2018 strain with a gray background; ‘−’: amino acid deletion; unique amino acids and nts shared by the CN/Liaoning25/2018 strain in the epitope region are indicated by a red square. (B) A large segment deletion of the CN/Liaoning25/2018 strain genomic sequence. The same large segment deletion in the CN/Liaoning25/2018 and attenuated CV777 strains.

When the complete genome sequence and different parts of CN/Liaoning25/2018 of genes were compared with those of the five PEDV reference strains, the CN/Liaoning25/2018 strain shared the highest nt identity with the attenuated CV777 with respect to the complete genome (98.8%), the ORF1a, and -b (ORF1ab) gene segment (99.86%), and ORF3_E_M_N genes (99.6%). However, the CN/Liaoning25/2018 strain shared the highest nt identity (98.7%) with the AH2012 strain regarding the S gene (Table 2). Overall these findings suggested that the CN/Liaoning25/2018 strain could be a recombinant strain.

3.3. Recombination and phylogenetic analysis

To determine whether the CN/Liaoning25/2018 strain was evolutionarily involved in recombination events, we used the recombination detection program version 4 (RDP4) software and compared the CN/Liaoning25/2018 strain with the 480 PEDV strains retrieved from the GenBank. All seven methods in RDP4 strongly supported that the CN/Liaoning25/2018 strain was a natural recombinant virus with an almost complete S gene from the CN/GDZQ/2014 strain, and the remaining genomic regions from the SQ2014 strain (Fig. 4A, Table 3). To further characterize recombination events and accurately determine parents, genome-scale similarity comparisons were made between the CN/Liaoning25/2018 strain and the reference strains with SimPlot v.3.5.1. This analysis further confirmed the chimeric nature of the CN/Liaoning25/2018 strain (Fig. 4B and C). Furthermore, a statistically significant signal for phylogenetic incongruence in the CN/Liaoning25/2018 strain detected by RDP defined two recombinant sources of the PEDV genome: 1, positions 1 to 2,597 and 24,758 bp to the end of the SQ2014 strain; 2, positions 20,597–24,758 bp from the CN/GDZQ/2014 strain—near the S gene of the CN/Liaoning25/2018 strain. Collectively, these results showed that the CN/Liaoning25/2018 strain evolved from a recombination event: the S gene was acquired from the CN/GDZQ/2014 strain and the remaining genomic regions were acquired from the SQ2014 strain (Fig. 4C).

Figure 4.

Recombination and similarity analysis of the CN/Liaoning25/2018 strain. (A) A rescaled structure of the CN/Liaoning25/2018 genome, inferred from the RDP analysis results. Light green: SQ2014 strain. Red: CN/GDZQ/2014 strain. The number indicates the approximate position of the melting point. The probable exchange region is shaded blue. (B) Bootscanning results for the putative recombinant CN/Liaoning25/2018 strain and its parental sequences. The parameters used were 200-bp window size, 20-bp step size, 1000 bootstrap replicates, Kimura two-parameter evolutionary distance model, and the neighbor-joining tree model. The arbitrary recombinant threshold was 70 per cent. (C) Similarity plot analysis with other PEDVs for reference. Phylogenetic analyses of the complete genome (D), nonrecombination regions (E) and recombination regions (F) of PEDV strains by BI. The same phylogenetic topologies were shown by both the ML and BI trees; only the BI tree is displayed. The nonrecombination region was at approximately 1–020,597 bp and 24,758 bp to the end, while the recombination region was at approximately 20,597–24,758 bp. Subgroups GI-a, GI-b, GII-a, GII-b, and GII-c are indicated using different colors. A black solid circle at the node indicates a posterior probability of 95% and above. The CN/Liaoning25/2018 strain and parental virus are indicated in red.

Table 3.

Summary of recombination events in the CN/Liaoning25/2018 strain identified using the RDP4 software.

Parental virus		Region	Breakpointsa		Score for the seven detection methods embedded in RDP4b
Major	Minor		Begin	End	RDP	GENECONV	BootScan	MaxChi	Chimaera	SiScan	3Seq
SQ2014	CN/GDZQ/2014	S gene	20,597	247,58	8.1 × 10−93	6.0 × 10−90	9.4 × 10−90	3.5 × 10−40	7.7 × 10−39	4.0 × 10−40	2.2 × 10−15

The breakpoints were based on the locations in the genome of CN/Liaoning25/2018.

The P-value cutoff was set at 0.01.

The phylogenetic tree showed that SQ2014 belonged to GI-b subgroup, which also includes the attenuated CV777 and DR13 strains; SQ2014 was closer to the attenuated DR13 strain, as per the evolutionary distance (Fig. 4D). Importantly, in parallel, SQ2014, and attenuated CV777 and DR13 strains all had the identical 49-bp deletion in ORF3, and no recombination signals were detected in the SQ2014 sequence. The data indicated that SQ2014 was derived from the low virulence vaccine isolates, consistently with a previous report (Zhang et al. 2017). CN/GDZQ/2014 belongs to the GII-a subgroup (Fig. 4D) and is a highly pathogenic strain (Song et al. 2014). Since the CN/Liaoning25/2018 strain was a highly pathogenic recombinant strain, the virulence was restored when the almost complete S gene from the low pathogenic vaccine SQ2014 strain was replaced by that of the highly pathogenic CN/GDZQ/2014 virus through recombination.

The nonrecombinant region, including positions 1 to 20,597 and 24,759 bp to the end of CN/Liaoning25/2018 fell within the GI-b subgroup and was close to SQ2014 (Fig. 4E). However, the recombinant region from 20,597 to 24,758 bp of the CN/Liaoning25/2018 strain fell into the GII-a group and was closer to CN/GDZQ/2014 (Fig. 4F).

4. Discussion

With the emergence of the variant PEDV in China in the latter part of 2010 (Li et al. 2012), the classical vaccine against PEDV could no longer provide satisfactory cross-protection against the new PEDV strains. This led to the rapid spread of the virus from China to the rest of Asia, the Americas, and even Europe, and has caused significant economic loss to the pig industry (Lee 2015). As with other coronaviruses, the diversity of PEDV is driven by genetic recombination (Guo et al. 2019). Several recombinant PEDV strains had been reported, but all recombination events have occurred between different highly pathogenic strains (Vlasova et al. 2014; Jarvis et al. 2016; Li et al. 2016; Chen et al. 2017; Nefedeva, Titov, and Malogolovkin 2019). Here, we report a new high pathogenic recombinant strain, CN/Liaoning25/2018, resulting from the recombination of the S gene from the high pathogenic CN/GDZQ/2014 strain with the remaining genes of the low pathogenic attenuated vaccine SQ2014 strain (Song et al. 2014; Zhang et al. 2017). CH/HNQX-3/14 strains reportedly arose because of the natural recombination between the attenuated vaccine strains (CV777, AF353511; DR13, JQ023161) and the circulating wild-type strain CH/ZMDZY11 (Li et al. 2016). However, using phylogenetic analysis we found that CV777 (AF353511) and DR13 (JQ023161) did not belong to vaccine strains, but instead were field strains (Fig. 4D), indicating that the CH/HNQX-3/14 strains also arose via recombination of two highly pathogenic strains. The emergence of the CN/Liaoning25/2018 strain led us to re-consider the clinical use of attenuated live vaccines. In clinical practice, attenuated live vaccines of PEDV are usually used to prevent PEDV infection. However, their often overlooked disadvantages include the possibility of the virulence recovery via mutation or recombination events. The vaccine isolate SQ2014, which is the parental virus of the CN/Liaoning25/2018 strain, is usually not pathogenic to piglets. In fact, the vaccine often protects piglets from PEDV infection. However, when the S protein of the SQ2014 strain was replaced by the S protein form the highly pathogenic CN/GDZQ/2014 strain through recombination, the recombinant CN/Liaoning25/2018 strain became highly pathogenic to piglets. The present findings demonstrate that recombination of PEDV live attenuated vaccine strains with a highly pathogenic PEDV field strain can occur in the natural environment, especially when the classical PEDV live vaccines do not provide complete protection (Crawford et al. 2016; Gerdts and Zakhartchouk 2017). The use of live attenuated vaccines for the prevention of PEDV should be considered more cautiously, because the live attenuated vaccines may play a role in the generation of new recombinant strains. This poses new challenges for the prevention and treatment of PEDV.

The current findings indicate that the S protein of PEDV is an important determinant of viral virulence (Zhang et al. 2015; Suzuki et al. 2018). In line with this, a study reported that partial or complete replacement of the S gene of attenuated CV777 by the corresponding S gene from the high-virulence OKN-After1/JPN/2013 strain led to increased virulence of the recombinant strain, similar to that of the high-virulent OKN-1/JPN/2013 strain (Suzuki et al. 2018). The natural recombinant CN/Liaoning25/2018 strain further confirmed this finding.

The 49-bp deletion in the ORF3 gene of classical vaccine attenuated strains (CV777 and DR13) compared with that of field strains is often utilized to distinguish classical vaccine from field strains. However, the present data suggest that if the S gene is changed by recombination (particularly with a pathogenic strain), the deletion of the ORF3 gene can no longer be used to distinguish between the classical attenuated and field strains.

We also wanted to estimate the time to the most recent common ancestor of the CN/Liaoning25/2018 strain. Unfortunately, this information could not be obtained because there were insufficient temporal signals in the data set when we used a regression of the root-to-tip genetic distances against the year of sampling in TempEst or a date-randomization test in the TipDatingBeast package (Rambaut et al. 2016; Rieux and Khatchikian 2017). In addition, the previously reported substitution rates could not be used in molecular clock calibration, as the temporal signal of their data was not tested. Perhaps, this is also the reason for the difference between their estimated substitution rates (Homwong et al. 2016; Jarvis et al. 2016; Jang et al. 2018). Furthermore, due to the inability to obtain the SQ2014 and CN/GDZQ/2014 strains, the recombinant strain in this study could not be reconstructed in the lab.

The emergence of the CN/Liaoning25/2018 strain suggests that the use of live attenuated vaccines against PEDV should be carefully considered; the same concern should be applied to other live attenuated vaccines, such as those to against TGEV, emerging SADS-CoV, PdCV, and even human coronavirus, such as middle-east respiratory syndrome-related coronavirus, severe acute respiratory syndrome coronavirus, and severe acute respiratory syndrome coronavirus 2. Additionally, our study also confirmed that the S protein is an important pathogenic protein of PEDV; that PEDV can cause intestinal epithelial cell vacuolation; and that the deletion of the ORF3 gene can no longer be used to distinguish PEDV vaccine from field strains.

Click here for additional data file.