Molecular characterization of the spike gene of the porcine epidemic diarrhea virus in Mexico, 2013-2016.

In Mexico, the first outbreaks suggestive of the circulation of the porcine epidemic diarrhea virus (PEDV) were identified at the beginning of July 2013. To identify the molecular characteristics of the PEDV Spike (S) gene in Mexico, 116 samples of the intestine and diarrhea of piglets with clinical signs of porcine epidemic diarrhea (PED) were obtained. Samples were collected from 14 farms located in six states of Mexico (Jalisco, Puebla, Sonora, Veracruz, Guanajuato, and Michoacán) from 2013 to 2016. To identify PEDV, we used real-time RT-PCR to discriminate between non-INDEL and INDEL strains. We chose samples according to state and year to characterize the S gene. After amplification of the S gene, the obtained products were sequenced and assembled. The complete amino acid sequences of the spike protein were used to perform an epitope analysis, which was used to determine null mutations in regions SS2, SS6, and 2C10 compared to the sequences of G2. A phylogenetic analysis determined the circulation of G2b and INDEL strains in Mexico. However, several mutations were recorded in the collagenase equivalent (COE) region that were related to the change in polarity and charge of the amino acid residues. The PEDV strain circulating in Jalisco in 2016 has an insertion of three amino acids (232LGL234) and one change in the antigenic site of the COE region, and strains from the years 2015 and 2016 changed the index of the surface probability, which could be related to the re-emergence of disease outbreaks.

Introduction

Porcine epidemic diarrhea (PED) is characterized by severe watery diarrhea, vomiting, and dehydration in all swine productive stages. In fattening pigs, an increase in the mortality rate of 2.5% has been registered [1]; however, mortality can reach 100% in piglets [2]. PED was first described on a pig fattening farm in England in 1971 [3]. From there, it disseminated to several countries in Europe with sporadic outbreaks, such as Italy in 2006 [4]. However, from 2014 to 2016, several emerging cases have been reported in other countries, including Germany, Belgium, Ukraine, France, Italy, and Austria [5-10]. In Asia, the first outbreak was reported in Japan in 1982; since then, PED has been considered a re-emerging disease [11]. The first outbreak in the Americas appeared in April 2013 in the USA, which disseminated rapidly throughout the country and caused severe economic losses [2]. In Canada, the first outbreak of PED was identified in January 2014, however in 2017, re-emerging outbreaks occurred in the province of Manitoba [12, 13]. The General Direction of Animal Health, SENASICA, in May 2014, recognized the presence of PED in Mexico, confirming its isolation and molecular detection in pigs with suggestive signs of the disease in 2016 [14]. The etiological agent of PED corresponds to the porcine epidemic diarrhea virus (PEDV), which was characterized for the first time in Belgium, in 1978 [15]. The virus is an enveloped, positive sense RNA virus that belongs to the Nidovirales order, Coronaviridae family, Coronavirinae sub-family, and Alphacoronavirus genus [16]. The PEDV genome has an approximate size of 28 kb and contains seven open reading frames (ORFs). ORF1a and ORF1b are located in region 5′; they occupy two-thirds of the genome and encode two polyproteins (pp1a and pp1b). The last third of the genome encodes four structural proteins: the spike protein (S), the envelope protein (E), the membrane protein (M), the nucleocapsid protein (N), and a non-structural protein encoded by ORF3 [17]. Similar to other coronaviruses, the S protein of PEDV plays a fundamental role in the interaction with the cellular receptor. Additionally, the S protein induces neutralizing antibodies in the natural host. This protein contains two domains, S1, which mediates the binding of the virus to the cellular receptor, and S2, which is involved in the fusion process of the virus to the cell. Within these two domains, there are four neutralizing epitopes: three of them are located in domain S1 (COE residues 499–638, SS2 residues 748–755, and SS6 residues 764–771), and the fourth is in domain S2 (2C10 residues 1368–1374) [18-20]. Furthermore, the S1 region contains the N-terminal domain (NTD) and the C-terminal domain (CTD) [21]. Phylogenetic analyses have grouped PEDV into two different genogroups (G1 and G2). Molecular analyses have identified the circulation of new viral strains; some of them are insertions or deletions of the S gene, called INDEL strains, and they are related to a diminished virulence of PEDV [22]. The objective of the present study was to perform the molecular characterization of the S gene of PEDV to identify the circulating strains from 2013 to 2016 in Mexico.

Materials and methods

Sample collection in PEDV outbreaks

We obtained 116 intestinal and diarrheal samples from piglets with clinical signs suggestive of PEDV infection (vomiting, watery diarrhea, dehydration, and anorexia) from 14 farms of multi-site intensive production systems located in six states of Mexico (Jalisco, Puebla, Sonora, Veracruz, Guanajuato, and Michoacán) from 2013 to 2016. The samples were collected in liquid nitrogen and then stored at −70 °C. RNA extraction of the samples was performed with commercial RNeasy Mini Kit (Qiagen) kits following the manufacturer’s protocol, and the real-time RT-PCR test described by Wang et al. [23] was performed to amplify a fragment of domain S1 of the S gene. This process, aside from determining the presence of genetic material of PEDV, also indicates whether the positive samples correspond to INDEL or non-INDEL strains. The samples corresponding to the years of 2013 and 2014 were analyzed with a differential test for transmissible gastroenteritis virus, porcine rotavirus A, and porcine deltacoronavirus (PDCoV), detecting positivity only for PDCoV in 35% of the samples analyzed (data not shown).

cDNA synthesis and PCR

Of the total samples positive in the real-time RT-PCR, we selected those presenting the highest amounts of genomic equivalents according to state and year of sampling. The cDNA was synthesized from extracted RNA using RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific) following the manufacturer’s protocol using random hexamer primers. Amplification of the S gene was performed through seven PCR reactions using the primers depicted in Table 1. For the amplification, we used GoTaq Green Master Mix (Promega); the amplification conditions were as follows: 95 °C for 10 min, 35 cycles of 95 °C for 30 s, annealing temperature dependent on primers, 72 °C for 60 s and a final extension of 72 °C for 10 min.

Table 1

Primers used for the amplification of the S gene of PEDV by RT-PCR

	Primer name	Sequence 5′–3′	Tm ( °C)	Position (whole genome)	Size (bp)
1	Fw-PEDV 2Fa	TAAGTTGCTAGTGCGTAATAATGA	56	20,568–21,326	758
	Rv-9Rb	GCAATTAGCTGTACAGGGTTCA
2	Fw-492 Fwc	GTTGTCGGCATAACATGGGA	56	21,055–21,865	810
	Rv-1302 Rvc	GTCAGTGCCATGACCAGTGAA
3	Fw-987 Fwc	GTCATTCTTGCTGAAGGCTC	56	21,550–22,305	755
	Rv1742 Rvc	ACAGATTGCAAGGTGAAAGG
4	Fw-1376 Fwc	CCATTCAGCGTATTCTTTATTG	56	21,934–22,778	839
	Rv-2215 Rvc	GGTAGAAGAAACCAGGCAACTC
5	Fw-S2-Ud	CTGATTCTGGACAGTTGTTA	60	22,631–23,457	826
	Rv-10Rb	GTACATGTGAAGCTTCTCAGCGT
6	Fw-2599 Fwc	AGTTCGTTTAATGGTGATGG	60	23,162–24,156	994
	Rv-3593 Rvc	GCGTAAACAAGACTAAGCCAGG
7	Fw-pspike 8a	AGCCATTTCTAGTTCTATTG	58	23,816–24,871	1055
	Rv-S2-ILd	TTGGACAGCATCCAAAGACA

a[34]

b[35]

cPresent study

d[36]

DNA sequencing

The amplified products were purified by means of commercial columns, QIAquick Gel extraction kit (Qiagen), following the manufacturer’s protocol. The purified products were sequenced in both directions using an ABI Prism 3100 xL Genetic Analyzer with the BigDye® Terminator v3.1 Kit (Applied Biosystems, Foster City, CA) in the Instituto de Biología of the Universidad Nacional Autónoma de México. The obtained sequences were compared to the GenBank database using the NCBI’s BLAST tool (https://blast.ncbi.nlm.nih.gov/Blast.cgi).

Phylogenetic analysis

The amino acid sequences were aligned and analyzed with ClustalW using MEGA 7.0.21 software. The antigenic index was calculated using the Jameson–Wolf algorithm [24], the surface probability profiles were obtained by the Emini method [25], and the hydrophilicity analysis was performed using the Kyte–Doolittle method [26]. With the amino acid sequences obtained from the present study and sequences of PEDV reported in the Americas, Asia, and Europe, a phylogenetic tree was constructed using the maximum likelihood based on the JTT matrix-based model and Gamma distribution with 1000 bootstrap replicates using the MEGA 7.0.21 software.

Database accession numbers

To accomplish the phylogenetic analysis, we used the S gene sequences from the different PEDV genogroups with the following access numbers to GenBank for G2: USA/Colorado/2013 (KF272920), MEX/104/2013 (KJ645708), MEX/Mexico329/2014 (KR265766), USA/Indiana/17846/2013 (KF452323), USA/Indiana34/2013 (KJ645641), PEDV/USA/Minnesota124/2015 (KY982981), COL/Cundinamarca/2014 (KU569509), IA1 (KF468753), AH2012 (KC210145), AJ1102 (JX188454), LC (JX489155), ZJCZ4 (JX524137), and CH/FJZZ-9/2012 (KC140102). For G1, we included sequences: CH/S (JN547228), SD-M (JX560761), attenuated DR13 (JQ023162), SM98 (GU937797), LZC (EF185992), and CV777 (AF353511). The INDEL strains were: L00721/GER/2014 (LM645057), 914 2015 AUT (KT895908), OH851 (KJ399978), 1842/2016 ITA (KY111278), and CH/HBQX/10 (JX501318). All the sequences obtained in the present study were also deposited in GenBank (Table 2).

Table 2

Characteristics of the PEDV strains analyzed in this study

Strain	Strain type	Region	Year of collection	Access number
PEDV/MEX/MICH/01/2013	No-INDEL	Michoacán	2013	KY828999
PEDV/MEX/SON/01/2014	No-INDEL	Sonora	2014	KY828992
PEDV/MEX/VER/01/2014	INDEL	Veracruz	2014	KY828991
PEDV/MEX/JAL/01/2014	No-INDEL	Jalisco	2014	KY828993
PEDV/MEX/MICH/01/2015	No-INDEL	Michoacán	2015	KY828995
PEDV/MEX/MICH/02/2015	INDEL	Michoacán	2015	KY828996
PEDV/MEX/PUE/01/2016	No-INDEL	Puebla	2016	KY829000
PEDV/MEX/GTO/02/2016	No-INDEL	Guanajuato	2016	KY828994
PEDV/MEX/JAL/01/2016	No-INDELa	Jalisco	2016	KY828997
PEDV/MEX/JAL/03/2016	No-INDELa	Jalisco	2016	KY828998

aStrains with an insertion of three amino acids (232LGL234)

Results

PEDV sequencing and analysis

Of the 116 obtained samples, 103 (88.8%) were positive for PEDV based on the real-time RT-PCR test. Of the positive samples, 87.4% corresponded to samples of the high virulence PEDV (G2 non-INDEL) and 12.6% to INDEL PEDV strains. We selected 10 samples for amplification and sequencing, which represents 74% of the farms sampled (10/14). The obtained sequences correspond to the six analyzed states of Mexico (Table 2). Alignment was performed using the previously reported sequences pertaining to G1a, G1b, G2a, G2b, and INDEL to compare the sequences obtained in this study. A neutralizing epitope analysis was performed through the alignment of amino acids. In the COE region (Fig. 1), seven of the ten analyzed USA/Colorado/2013 strains possessed mutations (Table 3); however, epitopes SS2, SS6, and 2C10 (Fig. 2) remain conserved. The insertion of three amino acids (232LGL234) was detected in strains PEDV/MEX/JAL/01/2016 and PEDV/MEX/JAL/03/2016 (Fig. 3). The alignment identified changes in specific residues according to the genogroups.

Fig. 1

Amino acid alignment of the PEDV sequences of the COE region. The sequences are grouped based on the results of the phylogenetic analysis. Circles indicate the sequences of the present study. The variable regions of the alignment are shown in yellow (Color figure online)

Table 3

Identification of amino acid changes in the COE region with respect to strain USA/Colorado/2013

Strain	Changes in region COE	Strain	Changes in region COE
PEDV/MEX/MICH/01/2013	L504P L618V	PEDV/MEX/MICH/02/2015	No changes
PEDV/MEX/SON/01/2014	H524R K566T	PEDV/MEX/PUE/01/2016	No changes
PEDV/MEX/VER/01/2014	No changes	PEDV/MEX/GTO/02/2016	T503I
PEDV/MEX/JAL/01/2014	S611R	PEDV/MEX/JAL/01/2016	G526D Y564F K569N G615V
PEDV/MEX/MICH/01/2015	K566T	PEDV/MEX/JAL/03/2016	G526D Y564F K569N G615V

Fig. 2

Amino acid alignment of the PEDV sequences. The neutralizing epitopes SS2, SS6, and 2C10 are show in red. Circles indicate the sequences of the present study. The variable regions of the alignment are shown in yellow (Color figure online)

Fig. 3

Insertion region of strains PEDV/MEX/JAL/01/2016 and PEDV/MEX/JAL/03/2016. Circles indicate the sequences of the present study. The insertion of three residues (232LGL234) is shown in red. The variable regions of the alignment are shown in yellow (Color figure online)

Analyses of potential antigenicity sites, hydrophilicity, and surface probability

The Spike protein was analyzed by means of the Jameson–Wolf method that predicts potential antigenicity sites; these sites are the result of changes in polarity and charge in the amino acid sequence of the different assessed strains. The prototype strains we used were the parental strain CV777 (G1), USA/Colorado/2013 as the high virulence strain (G2), and OH851 as the INDEL strain. The INDEL strain showed differences in antigenic sites compared to the high virulence strains in the approximate regions of residues 55, 180, and 225; however, strain PEDV/MEX/PUE/01/2016 showed the same difference in residues 180 and 225, despite being classified as a virulent strain. A diminished antigenicity was observed at the insertion site 232LGL234 of strain PEDV/MEX/JAL/01/2016. Conversely, in the COE region, an additional antigenic peak exists in the PEDV/MEX/SON/01/2014 and PEDV/MEX/JAL/01/2016 strains compared to those observed in the other assessed strains. The site found shows a variation in all analyzed strains at approximately residue 400 (Fig. 4). The surface probability analysis identified changes in the NTD region of the S1 domain in the Mexican strains from 2015 and 2016 compared to the strains of previous years that were analyzed in the present study. The hydrophilicity analysis reported a hydrophobic region we identified in the insertion region of the PEDV/MEX/JAL/01/2016 strain (Fig. 5).

Fig. 4

Antigenic analysis of the amino acid sequences of the Spike protein. The sites of the neutralizing epitopes COE, SS2, SS6, and 2C10 are shown. Red arrows indicate noticeable differences in the potential antigenicity sites. Circles indicate the sequences of the protein of the present study (Color figure online)

Fig. 5

Surface probability a and hydrophilicity b analysis of the amino acid sequences of the NTD S1 region of the Spike protein. Red arrows indicate differences between the sequences. Circles indicate the sequences of the protein of the present study (Color figure online)

Phylogenetic analysis of the S protein of PEDV

A phylogenetic analysis of the complete amino acid sequence of the S protein was performed to obtain information on the evolution of PEDV. The standard sequences identified the clades of the different genogroups. G1 comprises the strain CV777 and the Asiatic strains reported before 2010. Strains PEDV/MEX/VER/01/2014 and PEDV/MEX/MICH/02/2015 were classified within clade INDEL, which is integrated by strains reported in Europe and the American strain OH851. Strain PEDV/MEX/PUE/01/2016 is grouped with the strain USA/Indiana34/2013, whose complete genome is grouped in clade 1 of the USA clades. G2a is integrated by the Chinese strains reported in re-emerging PEDV cases in 2012. The Mexican strain previously reported and eight strains of the present study were grouped in clade G2b, which is composed of the ancestral strain AH2012 and the strains of the Americas (Fig. 6).

Fig. 6

Phylogenetic analysis of the amino acid sequences of the Spike protein of PEDV. The evolutionary history was inferred using the Maximum Likelihood method based on the JTT matrix-based model. A discrete Gamma distribution was used to model evolutionary rate differences among sites (5 categories). The tree is drawn to scale, with branch lengths corresponding to the number of substitutions per site. The analysis involved 34 amino acid sequences. All positions containing gaps and missing data were eliminated. There were a total of 1354 positions in the final dataset. Circles indicate the sequences of the present study (Color figure online)

Discussion

The first outbreaks of PED in Mexico occurred in July 2013. Since then, the disease has caused severe economic losses to the Mexican swine industry because of its high mortality, which reached 100% in piglets in the first outbreaks; in the re-emergence of outbreaks, the reported mortality has been 75% in newborn piglets. Until 2017, in Mexico, there was no vaccine available, only the feedback immunization strategy was implemented. Currently in Mexico, is authorized the killed virus vaccine based with the USA/Colorado/2013 (KF272920) strain. The PEDV/MEX/MICH/01/2013 strain originated in the state of Michoacán, where the first cases of the disease occurred; the S gene of that strain differs by only four nucleotides compared to the USA/Colorado/2013 strain, which was the first complete sequence of PEDV reported in the USA [27], indicating that the first PED outbreaks in México were caused by a related-highly virulent strain. Since 2013, the circulation of the first INDEL strain (OH851) was identified, that was closely related to the strain CH/HBQX/10 [22, 28]. The first Mexican INDEL strain was identified in Veracruz, Mexico and was named PEDV/MEX/VER/01/2014; it had a 99.8% identity with strain OH851, which is related to the reduced severity of the clinical disease in the piglets [22].

In PEDV, the receptor-binding domain (RBD) is located in the C-terminal end of region S1 of the Spike protein [29]; this region, named COE, is located in residues 499–638 of the strains of gene group 1 and INDEL and in G2 strains of high virulence, and in no-INDEL in residues 502–641. The COE region has been identified as a neutralizing epitope and as the region that recognizes the porcine amino peptidase N (pAPN), allowing entrance of the virus to the susceptible cell [18]. The analysis of this region of the protein identified diverse mutations in the obtained Mexican strains. These mutations result in changes of charge and polarity, which, in the case of the strains of Sonora and Jalisco (PEDV/MEX/SON/01/2014, PEDV/MEX/JAL/01/2016 and PEDV/MEX/JAL/03/2016) assessed by means of the analysis of the potential antigenicity sites, demonstrated a very evident change that corresponds to an additional antigenic peak compared to diverse PEDV strains. The 2016 strains from Jalisco were obtained from a farm where a PEDV re-emerging outbreak was reported with severe clinical signs of diarrhea in growing and fattening pigs. Additionally, in these strains, a three amino acid insertion (232LGL234) was found, which has not been reported previously at a worldwide level. This insertion contains two leucines; therefore, this region of the protein increases its hydrophobicity, as observed in the hydrophilicity analysis. The hydrophobicity, together with the changes in surface index in the strains of years 2015–2016 and mutations in the COE region, could be related to re-emerging outbreaks, tropism, and the presentation of the clinical PED signs [30]. Changes in the S gene have been related to recurring cases, as in Japan, where a co-infection of PEDV strains with the complete S gene and strains with deletions of amino acids 194–216 in region S1 were identified [31]. However, the neutralizing regions SS2 and 2C10 remain highly conserved, as reported previously in studies performed with strains from Japan in 2014 and China in 2015 [11, 32]. Region SS6 does not show mutations in relation to strains of G2b; however, there are changes compared to G1. These analyses of neutralizing epitopes together with the analysis of antigenicity sites provide basic information on the antigenic structure of the Spike protein. The phylogenetic analysis showed that the North American strains are located in clade G2b, with a common ancestor, strain AH2012, as reported previously, whereas the study of the complete genome of the North American strains reported until 2014 places them in three clades (North American clade I and II and USA INDEL), with two strains from Mexico in clades I and II [28]. Regarding the circulating strains of G2b in Mexico, the most distant with respect to time and substitutions of amino acids are PEDV/MEX/MICH/2013 and PEDV/MEX/JAL/01/2016, with a 99% identity. However, in 2015 in China, the S gene of diverse strains was assessed, yielding a homology between 96.6 and 100% [32]. This result indicates that the strains from this study have not undergone a high rate of mutations, although the antigenic sites have changed. Trujillo-Ortega et al. reported an analysis of a fragment of the S gene of four strains from a PED outbreak in Mexico in 2014 [14]. These strains seem to have a common ancestor, strain MEX/104/2013 (KJ645708), reported by Vlasova [28], and a 99% similarity with strains of the USA from 2013 and 2014, and even with a strain from South Korea from 2013. The PEDV/MEX/PUE/01/2016 strain was grouped outside the clade of G2b, due to probable recombination events, as described by Song et al. with strains CH/S and DR13 [33], so this recombination will be confirmed in the future analysis of complete genome sequences. However, in the present study, the phylogenetic analysis of the S protein demonstrated that the strains circulating in Mexico have different ancestors.

In conclusion, in this work, ten PEDV strains were identified from outbreaks in six regions of Mexico from 2013 to 2016. Of the four analyzed neutralizing epitopes, the Mexican strains presented mutations only in the COE region. In the 2016 strains of Jalisco, we identified a new insertion of three amino acids in PEDV, which has not been reported before at a worldwide level. In addition, the present mutations result in changes in the surface index and the antigenic sites of the COE region; therefore, the clinical presentation of the disease and the virulence can change with respect to the other circulating strains. Lastly, the phylogenetic analysis determined the circulation of PEDV INDEL strains of high virulence in the porcine population of Mexico. The results of this study help identify the molecular characteristics of PEDV, which will lead to improving disease control measures.