Evolutionary and genotypic analyses of global porcine epidemic diarrhea virus strains.

Porcine epidemic diarrhea virus (PEDV), which re-emerged in China in October 2010, has spread rapidly worldwide. Detailed analyses of the complete genomes of different PEDV strains are essential to understand the relationships among re-emerging and historic strains worldwide. Here, we analysed the complete genomes of 409 strains from different countries, which were classified into five subgroup strains (i.e., GI-a, GI-b, GII-a, GII-b, and GII-c). Phylogenetic study of different genes in the PEDV strains revealed that the newly discovered subgroup GII-c exhibited inconsistent topologies between the spike gene and other genes. Furthermore, recombination analysis indicated that GII-c viruses evolved from a recombinant virus that acquired the 5' part of the spike gene from the GI-a subgroup and the remaining genomic regions from the GII-a subgroup. Molecular clock analysis showed that divergence of the GII-c subgroup spike gene occurred in April 2010, suggesting that the subgroup originated from recombination events before the PEDV re-emergence outbreaks. Interestingly, Ascaris suum, a large roundworm occurring in pigs, was found to be an unusual PEDV host, providing potential support for cross-host transmission. This study has significant implications for understanding ongoing global PEDV outbreaks and will guide future efforts to develop effective preventative measures against PEDV.

INTRODUCTION

Porcine epidemic diarrhoea (PED) is a devastating enteric disease in pigs that results in severe diarrhoea, vomiting, and dehydration, with very high mortality observed in suckling pigs (Pensaert & de Bouck, 1978). Porcine epidemic diarrhea virus (PEDV), which is the causative agent of PED, belongs to the genus Alphacoronavirus within the family Coronaviridae and is an enveloped single‐stranded positive‐sense RNA virus (Woo et al., 2012). The PEDV genome consists of seven open reading frames (ORFs) organized in the order ORF1a, ORF1b, spike (S) glycoprotein gene, ORF3 hypothetical protein gene, envelope (E) gene, membrane (M) gene, and nucleocapsid (N) gene (Li et al., 2016). Among the proteins encoded by the ORFs, the S glycoprotein is located on the envelope of the virus in the large surface projections of the virion and plays an important role in the attachment of viral particles to host cell receptors (Lee, Park, Kim, & Lee, 2010). Thus, the S gene is considered important for understanding the genetic relatedness and epidemiological status of PEDV field isolates, as well as for advancing vaccine development (Chen, Liu, Lang, et al., 2013).

The PED disease was first discovered in pig farms in Belgium and the United Kingdom in 1976 (Pensaert & de Bouck, 1978), with reports of its occurrence in China as early as the 1980s (Li et al., 2012). With the emergence of new PEDV strains, however, serious disease epidemics have been observed in China since October 2010 (Sun et al., 2012). Beyond China, the disease has rapidly spread to more than 38 states in the USA following its first outbreak in May 2013, affecting more than 4,000 farms accounting for more than 7 million piglets (Cima, 2014). Japan, Canada, Mexico, and Colombia have also experienced successive outbreaks, with considerable economic losses to the global pig industry (Lara‐Romero et al., 2018; Ojkic et al., 2015; Takahashi, Okada, & Ohshima, 1983; Valko et al., 2017).

Major global outbreaks since 2012 have renewed concerns about the potential changes in the mode of PEDV transfer (Chen, Liu, Lang, et al., 2013; Li et al., 2012; Sun et al., 2012). Although increasing evidence suggests that PEDV routinely undergoes significant changes, especially in spike proteins (Lara‐Romero et al., 2018; Stott et al., 2017), the prevalence and evolution of PEDV strains is not well‐defined and limited knowledge is known regarding the ways in which PEDV subgroups circulate among themselves and how they might influence the evolution of PEDV. To better understand the molecular epidemiology and genetic diversity of PEDV field isolates, we investigated the genetic characterization, origin, and evolution of emergent PEDV strains worldwide, which will provide much needed information for the effective prevention and control of this disease.

MATERIALS AND METHODS

Sample collection

The complete PEDV genome was selected for genetic analysis. To clarify the evolution of PEDVs, we obtained four complete genome sequences of PEDV from our own lab (i.e., ZL29, AJ1102, HuB1‐2017, and HuB7‐2017) (Bi, Zeng, Xiao, Chen, & Fang, 2012) and 405 sequences (obtained 1 January 2018) from the National Center for Biotechnology Information (NCBI, https://www.ncbi.nlm.nih.gov) during the different time periods. Repeated sequences and culture‐attenuated PEDV sequence were removed from the data set. The final data set contained 409 sequences originating from Asia, Europe, and America. Details on the data set are summarized in Supporting information Table S1.

Phylogenetic analyses

Firstly, we performed multiple sequence alignment of the 409 complete PEDV genomes, as well as the ORF1ab, S, ORF3‐E‐M‐N genes, and applied the bat coronavirus BtCoV/512/2005 (GenBank accession no. DQ648858) sequence as an outgroup (Tang et al., 2006). A maximum‐likelihood (ML) phylogenetic tree was constructed using IQ‐TREE v.1.6.5 (Nguyen, Schmidt, von Haeseler, & Minh, 2015), with the best fitting evolutionary model suggested by the program following 1,000 bootstrap replicates. The phylogenetic tree was rooted against the PEDV‐related bat coronavirus, with removal of the long‐branch leading to greater resolution of the viruses of interest. Nucleotide and deduced amino acid sequences were aligned using MAFFT v.7.402 (Katoh & Standley, 2013). The resulting tree was visualized using iTOL v.4 (Interactive Tree of Life, http://itol.embl.de/).

Recombinant analyses

We preliminarily screened the PEDV sequence data set for recombination using RDP, GENECONV, Chimaera, MaxChi, and 3Seq, followed by secondary scanning and recombination using BootScan and SiScan in Recombination Detection Program version.4.95 (RDP v.4.95) (Martin, Murrell, Golden, Khoosal, & Muhire, 2015). Sequences with significant signals for recombination determined by more than two methods were analysed in greater detail. Nucleotide sequence similarity 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. All data were analysed using GraphPad Prism software (v.5.03, San Diego, CA, USA).

Bayesian phylogeographic analyses

A Bayesian framework was applied to reconstruct the spatiotemporal diffusion history of the PEDVs. In brief, the spatial diffusion of the time‐scaled genealogy was modelled as a standard continuous‐time Markov chain (CTMC) process over discrete sampling locations. Analysis was performed using BEAST v.1.8.2 (Baele, Lemey, Rambaut, & Suchard, 2017) under the assumption of a molecular clock model with uncorrelated lognormal distribution. The analysis was run for 200 million chains, with sampling every 20,000 generations.

Comparison of S protein sequences from different subgroups

All S protein sequences from the PEDV sample strains were analysed using the meta data‐driven comparative analysis tool (meta‐CATS) (Pickett et al., 2013), with a p‐value threshold of 0.05 (this threshold was the maximum probability level for the likelihood that the position differed among groups simply by chance), to identify significantly different sites between the five subgroups.

RESULTS AND DISCUSSION

Phylogenetic analyses of complete genomes of global PEDV strains

To characterize the genetic diversity of PEDVs circulating globally, we constructed a phylogenetic tree using IQ‐TREE based on the 409 complete PEDV genomes (see Materials and Methods 2.2). Consistent with our previous research (Wang, Fang, & Xiao, 2016a), the phylogenetic tree indicated that the complete PEDV genomes evolved into two separate genogroups, GI (classical) and GII (variant), as presented in Figure 1a. Furthermore, genogroup GI evolved into two subgroups (GI‐a and GI‐b) and genogroup GII evolved into three subgroups (GII‐a, GII‐b, and GII‐c). The GI‐a subgroup mainly included the earlier PEDV strains found in Europe and Belgium (virulent CV777) as well as some classical strains (LZC and SM98). The GI‐b subgroup included most cell culture‐adapted vaccine strains (attenuated CV777 and DR13) and other pandemic classical strains (AH‐M, SD‐M, and SQ2014) from China and South Korea. Since 2010, identification and sequencing analyses have shown that PEDV strain variants are highly prevalent worldwide. These prevalent variants, which have led to huge economic losses (Wang et al., 2016a), were mostly located within the GII genogroup. The GII‐a subgroup included strains from the USA as well as from countries that reported US‐like PEDVs (e.g., AH2012, HuB1‐2017, and HuB7‐2017) (Figure 1a). In contrast, the GII‐b subgroup predominantly consisted of strains from China, South Korea, and Thailand (AJ1102 and CH/SD2014) (Supporting information Figure S1a). The phylogenetic tree further indicated a new GII subgroup (GII‐c), which consisted of S‐INDEL strains from America (e.g., OH851) and Europe (e.g., GER/L00862/2014) with high degrees of sequence similarity to the ZL29 strain from China (this study, Supporting information Figure S1b), suggesting that S‐INDEL strains from Europe may have originated from a common ancestor with strain ZL29.

Figure 1

Genotyping and origin of the 409 PEDV strains based on full‐length genomic sequence analyses. (a) Phylogram was tested by 1,000 bootstrap replicates, branch lengths were measured by the number of substitutions per site (see scale bars). Names of strains, years, places of isolation, GenBank accession numbers, genogroups, and subgroups are shown. (b) Line chart shows the number of PEDV sequences obtained by gene subgroup and year of sampling. Yearly percentages of samples positive for PEDV are indicated by different coloured lines respectively. Data are indicated below sampling years. (c) Subgroup distribution of all available complete or partial PEDV genome sequences from countries reporting PEDV infections (n.a., sequence not available). In the bar charts, counts are shown by country or region. Data are indicated below bar charts [Colour figure can be viewed at http://wileyonlinelibrary.com]

Geographical and temporal analysis of PEDVs

We identified the geographical and temporal distributions of the PEDV strains to clarify the evolution of the virus. As shown in Figure 1b, only sporadic outbreaks of PEDV were reported before 2010, with the pathogens involved in these outbreaks found within the GI genogroup. The PEDVs were primarily located in the GI‐a (virulent CV777) and GI‐b (attenuated DR13) subgroups and were predominantly from Asia (Figure 1b). However, considerable PEDV outbreaks were reported in Asia and the United States after 2010, even for vaccinated piglets (Lin, Saif, Marthaler, & Wang, 2016). Based on our examination and assembly of public data, we identified that the PEDV strains were primarily from the GII genogroup.

Interestingly, GII‐b subgroup strains were reported predominantly in 2011, whereas GII‐a strains were reported more prevalently after 2011 and occupied a larger proportion of strains. Moreover, reports of the newly discovered GII‐c strains showed a significant increase after 2012 due to further sequencing from Europe. The above results indicate that the epidemic strains from different periods were from the five different subgroups. Based on the geographical distribution of PEDVs (Figure 1c), the GI genogroup (classical and cell culture‐adapted vaccine strains) largely originated from the earlier PEDV‐threatened areas, such as China, South Korea, and Europe. The different subgroups from the GII genogroup also showed characteristic geographical distribution. While most GII‐a subgroup strains were from the Americas, a small number were from China and Japan or from sporadic outbreaks in a few other isolated areas. The GII‐b subgroup strains were mostly endemic to Asia, especially China, South Korea, and Japan. The GII‐c subgroup strains were primarily from Europe, with some from the USA and China. Our study showed that PEDV strains from different subgroups were prevalent within the same areas, implying that the coincident “hot spots” in PEDV‐endemic areas (e.g., China and South Korea, Figure 1c) are critical for determining the sources of some PEDV variations. These “hot spot” areas have the potential to be important reservoirs for the genetic variation of PEDVs, resulting in recombination between different PEDV subgroups.

We also examined the potential hosts of the 409 PEDV strains. Results showed an unusual PEDV strain (GenBank accession no. KX883635) hosted by Ascaris suum (Supporting information Table S1) (Shi et al., 2016), a large roundworm in pigs, thus providing novel insight into the possible epidemiology of PEDV infection. Indeed, parasites have long been regarded as a harmful factor to the pig industry as sources for a variety of infectious agents (Jesudoss Chelladurai et al., 2017). For example, Metastrongylus larvae are considered a reservoir for various porcine viral pathogens, primarily swine fever virus and swine flu virus (Sen, Kelley, Underdahl, & Young, 1961). However, whether Ascaris suum plays a critical role as a PEDV reservoir requires further investigation.

Global emergence of recombinant PEDV

Previous studies have reported on the occurrence of recombination in different PEDV strains (Chen et al., 2017, Li et al., 2018; Boniotti et al., 2016). Here, the aligned complete genomes were all scanned for recombination using seven different algorithms implemented in RDP v.4.95. The RDP4 results were used to determine if the recombination parents of GII‐c were from the GI‐a or GII‐a subgroups. To further characterize recombination events and accurately determine parents, we performed genome‐scale similarity comparisons between the GII‐c and other subgroups with SimPlot v.3.5.1, as demonstrated in Figure 2a. This analysis confirmed the chimeric nature of the GII‐c subgroup. Furthermore, a statistically significant signal for phylogenetic incongruence in the GII‐c subgroup defined two recombinant sources of the PEDV genome: that is, (i) positions 1 to 20,551 and 21,551 to the end from the GII‐a subgroup, and (ii) positions 20,551 to 21,551 from the GI‐a subgroup. Collectively, these phylogenies showed that the GII‐c subgroup evolved from a recombination event, with the 5′ part of the S gene acquired from GI‐a and the remaining genomic regions acquired from GII‐a (Figure 2a). Molecular clock data analysis indicated that the GII‐c subgroup S gene likely originated in April 2010 (Figure 2b). As the re‐emergence of PEDV in China was reported in October 2010, the recombination event of the ancestral virus strain in the GII‐c subgroup appears to have transpired before the PEDV re‐emerging outbreaks.

Figure 2

Recombination and origin of the PEDV GII‐c subgroup. (a) SimPlot analysis for possible recombination events of the GII‐c subgroup genome compared to the other four subgroups. (b) Time‐scaled phylogeny of the PEDV S gene from the GII‐c subgroup. Time‐scaled Maximum Clade Credibility (MCC) trees were estimated from the complete S protein from the GII‐c subgroup, as indicated in Figure 2c, with tip times reflecting time of sampling (x‐axis). Node age estimates are shown, and the time line is indicated under the tree [Colour figure can be viewed at http://wileyonlinelibrary.com]

Phylogenetic analyses of ORF1ab, S, and ORF3‐E‐M‐N genes

To further explore the evolution of PEDVs, we constructed three phylogenetic trees based on the ORF1ab, S, and ORF3‐E‐M‐N gene sequences of the 409 PEDV strains. The ORF1ab and ORF3‐E‐M‐N gene alignments confirmed that the GII‐c subgroup was deeply nested within the GII‐a subgroup (Figure 3a and b). Strikingly, the phylogeny of the S gene suggested an entirely different evolutionary history for the GII‐c subgroup compared to the other subgroups (Figure 3c). All GII‐c subgroup strains showed inconsistent topology in the S gene phylogenetic tree, differing from the ORF1ab gene phylogenetic tree. This inconsistent topology, in which outlier sequences were found between two well‐defined subgroups in a phylogenetic tree, was attributed to genomic recombination. These results further demonstrated that the GII‐c subgroup evolved from the recombination of different subgroups.

Figure 3

Genotyping and origin of the 409 PEDV strains based on different genes. (a–c) Phylogenetic trees were constructed by the maximum‐likelihood (ML) method based on (a) ORF1ab, (b) ORF3‐E‐M‐N, and (c) S gene sequences, with 1,000 bootstrap replicates. Names of strains, years, places of isolation, GenBank accession numbers, genogroups, and subgroups are shown. (d) Model for the recombination and origin of the S gene in the five different PEDV subgroups [Colour figure can be viewed at http://wileyonlinelibrary.com]

Several viruses, such as virulent DR13 (isolated in 2009), Italy/7239/2009 (isolated in 2009), and KUPE21/2001 (isolated in 2001), which showed inconsistent topologies in the subgenomic trees, suggesting a more varied recombination history, were not classified within the five main subgroups. After multiple sequence alignment of the complete PEDV genome, virulent DR13 showed strong similarity with KUPE21 (95%–99% genetic identity to the ORF1ab and ORF3‐E‐M‐N genes). In contrast, the genetic identity between the S gene sequences of these viruses was only 85% and showed strong similarity with virulent CV777 (95% gene identity). In the S gene phylogeny, virulent DR13 did not cluster with KUPE21 but instead with virulent CV777, suggesting that the PEDV S gene was subjected to relatively frequent recombination, even between divergent subgroups (Figure 3c and Supporting information Figure S2a). Moreover, the Italy/7239/2009 PEDV strain (Chen, Liu, Shi, et al., 2013) may have originated from similar recombination events as those between KUPE21 and virulent CV777 (Supporting information Figure S2b).

Based on the phylogenetic and recombination analyses, the evolution of the PEDV S gene through time was reconstructed (Figure 3d), showing the influence of recombination. This model highlighted the genetic diversity expansion of the S gene, initially through antigenic drift and more recently due to widespread intergenotypic recombination. Since the emergence of Italy/7239/2009 in 2009, recombination has influenced all subsequently identified strains in the GII‐c subgroup. Intragenogroup recombination provides a mechanism for amalgamation among these distinct subgroups and increases the genetic repertoire of co‐circulating PEDV strains.

Analysis of amino acid sequences of complete S protein

The S protein attaches to the cellular receptors of a host, resulting in virus entry by membrane fusion, and contains the domain that stimulates the production of neutralizing antibodies. Variations in the S protein are important for understanding the genetic relatedness of PEDV field strains (Chang et al., 2002; Lara‐Romero et al., 2018; Wang et al., 2016b). To determine the significant S protein sequences among the five different PEDV subgroups, meta‐CATS analysis was performed for spike protein sequences of all 409 PEDV strains. We identified 257 amino acid positions with significant variation among the isolates from the five subgroups (Supporting information Table S2). As shown in Figure 4, the GII genogroup contained 11 distinct patterns of aa mutations (I116T, I356T, E365Q, T549S, G594S, N724S, A959V, S1044A, G1173D, S1232R, and R1298Q), distinguishing it from isolates in the GI genogroup. Although the GII‐c subgroup S genes shared some aa substitutions with the GI‐a and GII‐a subgroups (two parents of recombination), they also exhibited three unique patterns (L76, A/S92, and H/T113), which clearly distinguished these isolates from those in the other subgroups, suggesting that the GII‐c subgroup S gene evolved gradually through antigenic drift.

Figure 4

Amino acid variants of S proteins in different PEDV subgroups. (Upper) Amino acid changes in S proteins of the two gene groups, with 11 amino acid changes found to distinguish the two gene groups. (Lower) Amino acid substitutions in the S protein neutralization epitope regions of the five subgroups [Colour figure can be viewed at http://wileyonlinelibrary.com]

It is well established that the S protein of coronaviruses induces high levels of neutralizing antibodies. Four neutralizing epitopes (499–638, 748–755, 764–771, and 1,368–1,374 amino acids) have also been identified in the PEDV S protein (Chang et al., 2002). To explore whether the GII genogroup acquired substitutions in neutralizing epitopes characteristic of the GI genogroup, we mapped the significantly different positions to the equivalent positions in representative sequences (Figure 4). Results identified seven substitutions in the neutralizing epitopes of the S protein in the GII genogroup (L521H, S523G, V527I, T549S, G594S, A605E, and L612F), which may explain why traditional inactivated vaccines and attenuated vaccines against the GI genogroup cannot effectively protect piglets threatened by a pandemic strain from the GII genogroup.

In summary, this study revealed the genetic diversity and evolutionary dynamics of PEDV strains. Our genetic analyses showed that the PEDV strains could be categorized into two groups, namely, GI (classical) and GII (variant). We also discovered a new subgroup (GII‐c) with novel genetic, molecular, and phylogenic characteristics. The GII‐c subgroup evolved from a recombination event between the GI‐a and GII‐a subgroups, and we further found recombination in two relatively early strains: virulent DR13 and Italy/7239/2009. These recombination events occurred prior to the re‐emergence of PEDV in 2010. Additionally, to explore the potential link between S protein amino acid sequence variations and recombination, we performed a series of comparative analyses of the PEDV S protein sequences. We found 10 positions that were localized in a well‐known neutralizing epitope and revealed several unique amino acids that could easily distinguish the different subgroups. This study provides critical information to help trace the sources of PEDV variants and identify the evolutionary mechanisms involved. Furthermore, this research will hopefully facilitate the development of diagnostic kits, vaccines, and new therapeutic strategies, which are expected to turn the tide in the prevention of pandemic outbreaks of PEDV.

CONFLICT OF INTEREST

None.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.