Porcine epidemic diarrhea virus (PEDV) was first discovered in Europe and the first strain, CV777, was isolated in 1976 from pigs during an outbreak in Belgium [1]. PEDV continued to exist and spread in European countries in the 1970s and 1980s and afterwards. The number of recorded outbreaks has gradually decreased in this area [2]. In Asia, at this period, PEDV also appeared in some countries, such as Japan [3], China [4], and Thailand [5], causing severe economic loss. In the spring of 2013, the PEDV began to cause disease in the US, leading to the death of more than 8 million newborn piglets in just one year of the outbreak [6], which was named the US-like strain and soon spread to Canada, Mexico, and Columbia in the Americas [7]. Subsequently, US-like PEDV strains were reported in South Korea in late 2013 [8], in Germany in May 2014 [9], in Taiwan in late 2013 [10], in Japan in October 2013 and between 2013 and 2014 [11], and have soon become widespread, causing an epidemic wave in many Asian countries, such as Japan, China, Vietnam, and Thailand [7, 11–13]. These emerging PEDVs with their major-genetically modified novel descendants have initiated the second PED epidemiologic flow worldwide [7, 14]. Thus, the PEDVs were divided into two historic lineages: the “classical” strains since the 1970s and the “virulent” strains that emerged after the 2010s. The strains for these lineages have substantiated changes in the S region and in the genome [7].
The PEDV is an enveloped, positive-sense, single-stranded RNA virus that belongs to the genus Alphacoronavirus of the family Coronaviridae [1]. The genome of PEDV is about 28 kb in length, and contains seven open reading frames (ORFs) encoding for four major structural proteins (spike, S; envelope, E; membrane, M; and nucleocapsid, N), two nonstructural proteins, and one accessory protein, ORF3 [15]. Similar to other coronaviruses, the S protein of PEDV is the main envelope glycoprotein and a surface antigen that plays an important role in interacting with the host cell glycoprotein receptor during infection. Therefore, it contains the antigenic determinant regions that are targets of neutralizing antibodies in natural hosts [16]. The S protein has 1,383–1,386 amino acids (aa) with two subunits based on its homology with S proteins of other coronaviruses: S1 from position aa 1–789 and S2 from position aa 790–1383; and S1 can be further subdivided into five, e.g. S10 (aa 1–219), S1A (aa 435–485), S1B (aa 510–640), S1C, and S1D (aa 638–789) structural domains [17–19]. The S1 region is responsible for viral-host recognition and receptor binding, whereas the S2 region is for viral fusion and internalization. Mutations occurring in the S gene may result in amino acid changes leading to alterations in pathogenicity, transmission, antigenic properties, and virus-neutralization reactivity [20]. The S1 region functions as a receptor-binding domain (RBD) and has the ability to recognize a variety of host receptors, including proteins and sugars, that could trigger humoral and cellular responses. Meanwhile, the S2 subunit mediates the membrane fusion process during the release of viral RNA from infected cells [11, 20]. The S1 contains two main domains: the S1-N-terminal domain (S1-NTD) and the S1-C-terminal domain (S1-CTD) [17]. The S1 protein contains the highly conservative neutralizing epitopes CO26K equivalent (COE) (aa 499–638) [21] within the S1-CTD, which could produce PEDV neutralizing antibodies and serve as epitope candidates for vaccines against PEDV [22]. Therefore, the S1 subunit is very important to the host defense against PEDV.
The S gene of PEDV is subjected to frequent mutations under the pressure of vaccination and herd immunity, with an estimated substitution rate of 1.683 × 10−4 substitutions/site/year and 2.239 × 10−3 substitutions/site/year at nucleotide and aa levels, respectively [23]. As the S protein contains various immunodominant epitopes, this high level of mutation may affect the vaccine efficacy. Vaccines derived from classical strains in the G1 group are highly effective against “classical” PEDV strains but do not provide adequate protection against highly “virulent” strains in the G2 group [20]. In particular, the S1 subunit even has a higher mutation rate than the overall rate of the full S gene, at 1.5 × 10−3 substitutions/site/year [23]. Molecular studies of the S1 region are beneficial for the understanding of immune responses to genetic variation.
The PEDV was detected in Vietnam in 2008 and rapidly spread throughout the country, causing severe economic loss [24], and the PEDVs are now persistent and circulating in the swine populations [25, 26]. Studies on the S gene sequences of PEDV strains isolated during 2012–2016 showed that these strains had undergone remarkable changes [11, 27, 28], which might explain the decrease in vaccine efficiency. Moreover, PEDV continues to predominate in the pig populations, even in the vaccinated herds, and has become a persistent endemic nationwide [25, 26]. In Vietnam, imported attenuated PEDV and locally produced inactivated vaccines have been used to immunize pigs and sows. These include the imported “Avac PED live” vaccine (AVAC company, Vietnam) and attenuated “PED Pig VAC” vaccine imported from Daesung Microbiological Labs (Korea) that are all produced from the Korean G1a-SM98 strain (GU937797/KOR/SM98/2010) [29]. Besides, the inactivated vaccines were also used, including the combined “Provac TP” produced from the SM98P and 175L strains, respectively, against porcine epidemic diarrhea and transmissible gastroenteritis imported from the KOMIPHARM International Co., Ltd (Korea), and the local “PED vaccine” manufactured by the HANVET company in Vietnam (http://hanvet.com.vn/vn/Scripts/default.asp). Although vaccinated, severe outbreaks of PED have been continually occurring in the pig herds throughout the whole country. The reasons for vaccine failure may be sought through insights into the genetic variations within the genomes of the currently circulating PEDV strains that continue to generate new descendants with incompatible antigenicity [27, 28, 30, 31]. Molecular analysis of the clinical PEDV strains recently or currently circulating in Vietnam is therefore needed. In this study, we report the molecular characteristics of the S1-CTD domain with the variations in the core neutralizing epitope (COE) of the S1 protein from 26 PEDV strains isolated in seven provinces of Northern Vietnam in 2018 and 2019. The molecular data obtained can be used to select appropriate vaccine strains for controlling the PEDVs that are circulating in Vietnam.
Fifty samples, including small intestine tissue and feces, were collected in 7 provinces from June 2018 to January 2019. Three clinical samples were collected from Ha Noi (21°01′42.5″N, 105°51′15.0″E), two from Hai Duong (20°56'27.56"N, 106°19'58.87"E), eight from Hung Yen (20°38′46.93″N, 106°03′4.03″E), three from Nam Dinh (20°26'19.7628''N, 106°9'43.578''E), four from Thanh Hoa (19°48'16.19"N, 105°46'20.99"E), three from Tuyen Quang (22°10' 21.6156''N, 105°18'47.2248''E), and three from Vinh Phuc (21°19'59.99"N, 105°34'0.01"E) (Table 1). The collected samples were transported to the laboratory and stored at –20°C until used.
Fecal or interstinal samples were homogenized in sterile water and centrifuged at 13,000 rpm for 15 min. Total RNA was extracted using Trizol reagent (Merck, Darmstadt, Germany) according to the manufacturer’s instructions. The RNA was suspended in DEPC-treated water and stored at –80 °C. The first-strand cDNA synthesis was performed using the Maxima Reverse Transcriptase Kit (Thermo Fisher Scientific Inc, Waltham, MA, USA) following the manufacturer’s protocols. A pair of oligonucleotide primers were used to screen for PEDV positive samples, including PEDV_F: 5’-TTCTGAGTCACGAACAGCCA-3’ (S gene, positions 1475–1494) and PEDV_R: 5’-CATATGCAGCCTGCTCTGAA-3’ (S gene, positions 2106–2125), amplifying a fragment of the S1-CTD region (651 bp in length). PCR was performed on a C1000 Thermal Cycler (Bio-Rad, California, USA) with one cycle of initial denaturation at 95 °C for 10 min, followed by 35 cycles of denaturation at 95 °C for 30 s, annealing at 55 °C for 30 s, extension at 72 °C for 1 min, and final extension at 72 °C for 10 min.
The PCR products (651 bp) were purified with the GeneJET PCR Purification Kit (Thermo Fisher Scientific Inc.) and sequenced directly or after cloning into the pCR2.1-TOPO TA-cloning vector (Invitrogen, USA) by a service company. The resulting sequences were assembled using BioEdit software (http://www.mbio.ncsu.edu/BioEdit/bioedit.html) and used for searching using BlastN and BlastX tools in the NCBI database. The partial S1 sequence of 630 nucleotides obtained from each sample in this study and reference sequences from previous studies and from the GenBank database were used for molecular analysis.
An alignment of 61 partial S1 nucleotide sequences (about 630 bp each of the S1-CTD region), including 26 partial S1 sequences of PEDV isolated in North Vietnam (Table 1) and 35 reference strains of G1a (four strains), G1b (three strains), recombinant/G1c (12 strains), G2a (six strains), and G2b (eight strains) (Table S1), was performed using the GENEDOC 2.7 program (http://iubio.bio.indiana.edu/soft/molbio/ibmpc/genedoc-readme.html) for sequence comparisons and phylogenetic analysis. For phylogenetic reconstruction, MEGA X (www.megasoftware.net) was used with the maximum likelihood method and 1000 bootstrapping replications [32]. Amino acid sequences of 26 Vietnamese and one G1a (CV777, 1a-AF353511/BE/CV777/2001) reference strains were also aligned for genotypic and antigenic comparisons.
Phylogenetic analysis was performed based on the alignment of 61 PEDV S1 partial sequences (about 630 bp of the S1-CTD region, position: 1480–2110), including 35 representative strains of genotype or genogroup 1 (G1) and genotype or genogroup 2 (G2), and some prototype vaccine strains and 26 sequences from this study (listed in Tables 1 and S1). The (sub) genotype-representative strains included in the alignment were the prototype strain Belgium/CV777/2001 (AF353511) and the vaccine strain KR/SM98-5P/1998 (KJ857455) for subgroup G1a, the attenuated CV777 vaccine strain (JN599150) for subgroup G1b, the virulent CH/HNAY/2015 strain (KR809885) for subgroup G2a, and the vaccine strain AJ1102 (JX188454) for subgroup G2b. In addition, some recombinant strains reported recently in the United States (for example, strain USA/Iowa/2013), in China (such as ZL29, CH/SCZY44/2017) and Taiwan (TW/Yunlin550/2018, suggested as G1c) as US-like and/or recombinant strains were also included for reference. The ML-phylogenetic tree shown in Fig. 1 indicated that all of the strains were divided into two major groups: G1 (classical strains) and G2 (virulent pandemic strains) [7]. A cluster, which was formed from seven strains of the Vietnamese PEDVs in this study, was placed in between the G1 and G2 genogroups, with distinct amino acid changes.
The phylogenetic tree classified the 26 Vietnamese PEDV strains into G2b (19/26 = 73.1%) and a distinct subgroup (7/26 = 26.9%). Nineteen strains collected in 2018 belonged to G2b. They included the strains isolated in Thanh Hoa province (n = 4; MW345261/VN/TH/2018; OM681601/VN/TH-1/2018; OM681602/VN/TH-2/2018; and OM681603/VN/TH-3/2018), in Vinh Phuc (n = 3; MW345259/VN/VP/2018; OM681606/VN/VP-1/2018; and OM681607/VN/VP-2/2018), in Nam Dinh (n = 3; MW345262/VN/ND/2018; OM681599/VN/ND-1/2018; and OM681600/VN/ND-2/2018), in Hung Yen (n = 8; MT198679/VN/IBT/2018; OM681592/VN/HY-1/2018; OM681593/VN/HY-2/2018; OM681594/VN/HY-3/2018; OM681595/VN/HY-4/2018; OM681596/VN/HY-5/2018; OM681597/VN/HY-6/2018; and OM681598/VN/HY-7/2018), and in Hai Duong (strain OM681587/VN/HD-1/2018). They were placed in a group with the reference strain MK584552/CN/AJ1102(F12)/2011 and six previously reported Vietnamese and three Chinese G2b strains (Table S1; Fig. 1). Nine of these strains belonged to the same branch G2b as the AJ1102 vaccine strain, while ten belonged to a different cluster that included previously identified Vietnamese G2b strains isolated between 2014 and 2016 [13, 28]. None of the currently studied strains were placed in the G2a and G1c (or recombinant) clusters [33]. Notably, seven strains fell into the basal position, which is bracketed between G1b and the major group of the genotype G2a, G2b, and G1c (recombinant), separate from both G1a and G1b with a high bootstrap value (98%). This branch included three strains isolated in Ha Noi (n = 3, OM681589/VN/HN-1/2018; OM681590/VN/HN-2/2018; and OM681591/VN/HN-3/2018); three strains in Tuyen Quang (n = 3, MW345260/VN/TQ/2019; OM681604/VN/TQ-1/2018; and OM681605/VN/TQ-2/2018) and one strain in Hai Duong (OM681588/VN/HD-2/2018). We named this cluster a distinct or novel sub-genogroup and designated it as the “sG” (referred to as a distinct sub-genotype) (Table 1; Fig. 1).
To evaluate the antigenic similarity with the vaccine strain currently used in Vietnam, the S1-CTD amino acid sequences spanning from amino acid position 495–700 of the 26 studied strains were compared with the classical strain CV777 (AF353511) (Fig. 2). This region includes the COE domain (sites 499–638) [21], and a part of the S1D (sites 638–733) [34]. Table 2 showed a total of 19 amino acid substitutions for all studied Vietnamese strains in comparison with the CV777 strain. Thirteen positions of variation, including 513T>S, 520G>D, 521L>(Y/H), 523S>G, 525N>D, 527V>(I/L/M), 549T>S, 567S>G, 591L>F, 594G>S, 605A>E, 612L>F and 635I>V, fell into the COE domain, and other six (665S>A, 665I>F, 668L>F, 669A>(S/P), 672Y>H, and 665V>I, were found in the S1D region. Among these, five of them were potential, namely 521L, 523S, 527V, 605A, and 635I as previously identified neutralizing epitopes [35].
Additionally, to examine the conserved, variable, and neutralizing amino acids in the S1-CTD region, an alignment of sequences from 61 strains, including 26 from this study and 35 representative genotype-reference strains, e.g., three G1a, four G1b, seven G2a, eight G2b, and 12 recombinant (G1c) strains (for sequences, see Table 1 and S1), was also performed (Fig. S1). As result, all sequences from G1b, 2a, 2b, 2c, and the distinct sub-group (our studied strains) had the amino acid change in three positions in the COE (521L>(Y/H), 523S>G, 605A>E, and 635I>V), and one in the S1D region (667I>F) from G1a, whereas 527V>I, 549T>S, 594G>S, and 612L>F were found between the majority of the G2 (2a, 2b) and G1 (1a, 1b) strains. Notably, some studied Vietnamese G2b strains had distinguishing aa changes at 521L>I, 525N>D, 567S>G, 668L>F, and 672Y>H, and seven strains of the distinct subgroup had unique aa changes at 513T>S, 520G>D, 527V>(L/M), 591L>F, 669A>(S/P), and 691V>I. All sites of amino acid substitutions that occurred in the studied strains were noted and summarized (Figs. 2 and S1; and Table 2). These results suggest that the Vietnamese PEDV strains in this study might have changed remarkably in comparison with the vaccine strains being used on pig farms.
Among the immunodominant regions of the S protein, the neutralizing epitope COE (aa 499–638) is the most conserved domain and is capable of inducing PEDV-neutralizing antibodies and suitable for genogroup-phylogenetic analysis [21, 33, 34]. Therefore, in the present study, we selected the sequence in the S1-CTD domain that contains the COE to find genetic variation and predict a new PEDV variant that may emerge for strains circulating in Northern Vietnam (2018–2019). Phylogenotypic analysis has been conducted with 26 sequences from clinical samples and 35 reference strains representing all genogroups available to date (Fig. 1). With strong nodal support (98%), the topology showed the majority (73.1%) of strains (19/26) belonged to the subgroup G2b and seven strains clustered into an independent, “distinct subgenogroup”, separate from both G1a and G1b and other G2a, 2b, and 1c (recombinant strains).
Seven of our PEDV viruses were classified into a distinct subgroup of which the amino acids have been found changed at six antigenic sites in the S1-CTD protein sequences investigated. This new antigenic variant subgroup revealed the continuous evolution in the S1 region of the currently circulating Vietnamese PEDV strains [13, 28], which emphasizes the need for frequent updates of the vaccines for effective protection. Vaccination is an effective method to reduce the loss of morbidity and mortality in piglets infected with PEDV, and in particular, vaccination for pigs during pregnancy will increase the ability to protect newborn pigs [36]. Subunit vaccines using whole or partial spike proteins, particularly the COE domain, are safe and effective in triggering immunization [36]. The occurrence of various mutations in the COE domain due to antigenic drift as well as INDEL or recombination affects the antigenicity of the strain, and may be selected as potential markers for identification, molecular epidemiology, pathogenicity, phylogenetic analysis, molecular genotyping, antigenic studies, and evaluation of vaccine efficacy [20, 22, 29, 36, 37].
In comparison with the prototype strain CV777 (AF353511), we found a total of 19 aa substitutions in this domain for the 26 isolated strains, among which aa changes at six positions were observed in all strains (521L, 523S, 527V, 605A, 635I, and 667I). These substitutions were also summarized in another study on global PEDV strains [34]. Particularly in the COE domain, three substitutions were observed in all Vietnamese G2b strains (527V>I, 549T>S, and 594G>S). Similar changes at positions 549 and 594 were found in the PEDV strains isolated after 2011 in China, Korea, and the US [37], but most strains isolated in Korea from 2002–2009 had an R residue at position 549 [38]. The substitution sites in four potential epitopes observed in the studied strains (521L, 523S, 527V, 605A, 635I) were also reported in the study of [35]. These variations might not generate novel epitopes or alter the efficacy of existing epitopes. These variations belonged to the coiled region or at the end of β-sheets [35], not exposed to the surface, hence not affecting the neutralizing activity. It is noteworthy, however, that in comparison with the attenuated CV777 (JN599150) and AJ1102 (MK584552) vaccine strains, apart from the aa changes that were found in previous studies [13, 28], Vietnamese PEDV strains from different geographical regions in this study also possessed some distinct mutations, such as 513T>S, 520G>D, 521L>Y 525N>D, 527V>(L/M), 567S>G, 591L>F in the COE, and 669A>(S/P), and 672I>H in the S1D region (Table 2). Changes in a range of epitopic residues of the immunodominant regions may contribute to redistribution of the capsid surface with new conformational protruding epitopes and to alter the conformational specificity for immunity and receptor binding, leading to the development of novel antigenic variants that help the virus escape the host immune system. The molecular data indicated that the Vietnamese PEDV strains might have experienced dramatic antigenic drift, which ultimately led to the reduced protection by classical strain-based vaccines and resulted in a relatively high infection rate among pig farms.
According to previous evolutionary analyses based on the nucleotide sequences, PEDV is divided into two groups, the classical G1 group and the pandemic G2 group. Each group can further be divided into at least two subgroups, including G1a, G1b, G2a, and G2b. The G1a subgroup includes the CV777 (JN599150/CN/attenuated-CV777/1994) and DR13 (DQ862099/KR/DR13/2005) cultured vaccine strains. The G1b subgroup includes new variant strains discovered first in China, then in the United States, Korea, and more recently in European countries [39, 40]. The G2 strains appeared due to point mutations from G1a strains. While most of the G2a subtype originated mainly in the United States and a small part of China, Japan, or individual regions with small outbreaks, the G2b subgenogroup is the major cause of the pandemic outbreaks in Asian countries, especially China, Japan, and Korea [14, 34, 39]. Recently, many recombinant PEDV strains have been reported globally. The G1c subgroup contains the recombinant strains between the G1a and G2a subgroups, which include the strains identified in the US (such as USA/Iowa106/2013) and in Europe, with high sequence identity to the Chinese ZL29/2015 strain [34]. Meanwhile, the strains in the S-INDEL G1c subgroup resulted from the recombination events between a G2b strain and a wild-type G1 strain (such as strain TW/Yunlin/2018) [33] or a variant G1 strain (such as strain CH/SCZY44/2017) [35]. These recombinant strains were also included in our study as references. However, as the C-terminal region of S1 PEDV G1c strains was highly similar to that of G2b strains [33], the phylogenetic tree based on the S1-CTD region used in this study placed these recombinant strains close to other G2b strains. Nonetheless, all the recombinant variants were closely grouped together.
Regarding our 26 PEDV strains, 19 (73.1%) of them landed in the G2b subgroup, while the remaining seven PEDV strains formed an independent branch in the G1 group or between G1 and G2, clearly separate from the G1a and G1b reference strains. Before 2014, in Vietnam, the strains of G1b were predominant, as reported in previous studies [13, 28, 29]. In 2015–2016, the G2 group totally dominated worldwide with the presence of both subgroups G2a and G2b [28]. As depicted on the phylogenetic tree, the Vietnamese G2b PEDV strains isolated from 2014 onwards appeared to have undergone a certain degree of evolution and were located in separate clusters from the original AJ1102 G2b strains isolated in 2011–2012. Therefore, the G2 strains now appear to be the major threat of PEDV infection in Vietnam, as the classical vaccine strains provide limited protection against the non-S INDEL strains [41, 42]. Furthermore, ten G2b strains in this study, together with some other Vietnamese PEDV strains from previous reports, evolved in a different cluster far from the AJ1102-based vaccine strains (Fig. 1). This highlights the need for frequent vaccine updates to ensure protection efficacy despite the continuous evolution of PEDVs.
In this study, a remarkable finding was that seven strains (MW345260, OM681604, OM681605, OM681588, OM681589, OM681590, and OM681591) formed a separate branch in the G1 group or were basally placed between the G1 and G2 subgroups (Fig. 1). Besides, alignment of the S1-CTD sequences with the representative G1, G2 and recombinant strains revealed some unique patterns of variation in these seven strains. Previous studies on the PEDV S gene mostly focused on the four eminent neutralizing epitopes (COE aa 499–638, SS2 748–755, SS6 764–771 and 2C10 1368–1374). The analyses in our study included the region spanning from position 499–700. We found substitutions at two positions (669A>(S/P and 691V>I) that were absent in the reference strains. These two positions are located within the broad epitope region of an identified neutralizing monoclonal antibody 8A3A10 [18]. Besides, seven strains in our study had other aa variations at four antigenic sites (e.g., 513T>S, 520G>D, 527V>(L/M), and 591L>F) in the COE domain (Figs. 1 and S1). These amino acid changes at the neutralizing epitopes in the COE domain have not been reported in any strains studied so far, and thus they are restricted to the emerging Vietnamese PEDVs, representing a new variant antigenic sites in the S1 spike protein. In brief, the genetic polymorphisms observed in our PEDV strains all occurred at the important sites of the S protein. Based on the combined analyses of phylogeny, nucleotide and amino acid sequences, we would like to propose these seven strains as the members of a new PEDV subgroup (a distinct subgenogroup). Further studies should be conducted to gain a better understanding of the evolution of the Vietnamese PEDV strains and investigate the genetic variations for the selection of highly protective vaccine designs.
In conclusion, phylogenetic analysis and comparisons of 26 Vietnamese PEDV-S1-CTD sequences from samples collected from seven northern provinces from June 2018 to January 2019 revealed the presence of two genogroups of virulent strains in Vietnam: diverse G2b and an independent, “distinct subgenogroup”. They are separate from G1a and G1b and other G2a, 2b, and 1c (recombinant strains) and possess distinct changes of 513T>S, 520G>D, 527V>(I/L/M), 591L>F, 669A>(S/P), and 691V>I in the COE and S1D regions that were only restricted to these Vietnamese strains. PEDVs in Vietnam continue to evolve, which may cause vaccination failure in the pig herds immunized with the classical vaccines. This emphasizes the need for frequent updates of the vaccines for effective protection in Vietnam and other countries.