Genomic investigation of piglet resilience following porcine epidemic diarrhea outbreaks.

Porcine epidemic diarrhea virus (PEDV) belongs to the Coronaviridae family and causes malabsorptive watery diarrhea, vomiting, dehydration and imbalanced blood electrolytes in pigs. Since the 1970s, PED outbreaks have become a source of problems in pig producing countries all over the world, causing large economic losses for pig producers. Although the infection in adults is not fatal, in naïve suckling piglets mortality is close to 100%. In this study, we investigated genome-wide differences between dead and recovered suckling piglets from commercial farms after PED outbreaks. Samples from 262 animals (156 dead and 106 recovered) belonging to several commercial lines were collected from five different farms in three different countries (USA, Canada and Germany) and genotyped with the porcine 80K SNP chip. Mean Fst value was calculated in 1-Mb non-overlapping windows between dead and recovered individuals, and the results were normalized to find differences within the comparison. Seven windows with high divergence between dead and recovered were detected-five on chromosome 2, one on chromosome 4 and one on chromosome 15-in total encompassing 152 genes. Several of these genes are either under- or overexpressed in many virus infections, including Coronaviridae (such as SARS-CoV). A total of 32 genes are included in one or more Gene Ontology terms that can be related to PED development, such as Golgi apparatus, as well as mechanisms generally linked to resilience or diarrhea development (cell proliferation, ion transport, ATPase activity). Taken together this information provides a first genomic picture of PEDV resilience in suckling piglets.

Porcine epidemic diarrhea virus (PEDV) belongs to the Coronaviridae family; this family includes a wide variety of viruses that affect humans and other animals, causing respiratory and gastroenteric diseases (Weiss & Navas‐Martin 2005). As with many coronaviruses, PEDV has a limited host range and tropism, infecting only the small intestine of pigs (Song & Park 2012; Jung et al. 2014). This infection causes acute intestinal disease during which the infected enterocytes rapidly develop necrosis, leading to a villous atrophy (Debouck et al. 1981; Ducatelle et al. 1982; Jung et al. 2014) and development of malabsorptive watery diarrhea, vomiting, dehydration and imbalanced blood electrolytes (Jung & Saif 2015). The major route of transmission for PEDV infection is fecal–oral (Turgeon et al. 1980; Utiger et al. 1995; Riley 2007) even if aerosolized PEDV remains infectious (Alonso et al. 2014). It is known that PEDV, as well as other swine coronaviruses, such as porcine TGEV (transmissible gastroenteritis virus), gains entry in the host cells through interaction with alanyl (membrane) aminopeptidase, encoded by the ANPEP gene (Oh et al. 2003; Li et al. 2007). This receptor is largely expressed by enterocyte cells, and its density has been correlated to the replication rate of the virus in vitro (Nam & Lee 2010). Since discovery of the virus during the 1970s in Europe (Pensaert et al. 1978), PED outbreaks have become an increasing source of problems in many swine breeding countries all over the world, causing severe problems as well as large economic losses for pork producers (Song & Park 2012; Snelson 2014). It is widely understood that, although the infection in adults may not be fatal, PEDV infection in unprotected piglets under 3 weeks of age results in severe disease with mortality reaching 100% on many farms (Pensaert et al. 1978; Shibata et al. 2000). The degree to which the small percentage of naïve suckling piglets recover during PED outbreaks in the wider industry is unknown, as is the biological mechanisms involved, but two main hypothesis can be suggested: (i) survival can be related to variation in the intestinal receptor used by PEDV to gain entry to intestinal epithelial cells (the ANPEP gene) or (ii) as summarized by Schneider & Ayres (2008) and Ayres & Schneider (2012), survival and recovery can be related to particular host immune responses that enhance viral clearance, increase cell epithelial regeneration rate or reduce viral replication rate.

In this study, our aim was to investigate genome‐wide differences between dead and recovered suckling piglets from commercial farms during PED outbreaks. The study was conducted using piglets younger than 3 weeks of age from farms during the acute phase of a PED outbreak (as soon after herd diagnosis as physically possible). As depopulation is one of the main approaches used to deal with outbreaks, it is difficult to obtain good information on survival rates, or indeed samples, to test these hypotheses. To address this, we collected samples from 262 animals belonging to several commercial terminal crossbred lines from five different farms in three different countries (USA, Canada and Germany). Among these animals, 156 animals died as a result of severe diarrhea and dehydration and 106 animals successfully recovered (Table S1). Specific piglets were collected to ensure that balanced numbers of dead and recovered piglets were collected from each farm, with the exception of one farm in the USA from which the numbers were not balanced. Moreover, in three of the five farms (one farm located in USA and the farms located in Germany) dead and recovering piglets belonged to the same litters (1–2 dead and 1–2 recovered per litter). Each animal was genotyped with the GeneSeek Genomic Profiler 80K SNP chip (Geneseek–Neogen), which contains 68 517 SNPs. For the analysis, only SNPs mapped to the autosomes and with a call rate >90 were considered. In total, 57 422 SNPs passed the filtering steps, and missing SNPs were imputed using beagle 3.3.3 (Browning & Browning 2007). admixture software (Alexander et al. 2009) was used to investigate population stratification. Following the software guidelines, the cross‐validation (CV) procedure was applied to choose the best K for the model, where K is the number of (sub)populations that was assumed for the analysis. Within a range of 10 K values, 1–10, K = 8 showed the lowest CV error and therefore was considered a sensible model. The admixture plot considering K = 8 is reported in Fig. S1. As expected, the samples appeared to be stratified, with more similarity between the two US farms. Due to the expected low number of recovering piglets from each farm, detailed analysis by farm or within breed or population strata could not be considered. Therefore, we focused our attention on a comparison of all dead with all recovered, attempting to detect common windows of divergence. A multidimensional scaling plot was obtained using plink 1.7 (Purcell et al. 2007), comparing dead and recovered piglets (Fig. S2). The plot clearly shows a lack of clustering between the two groups. To find windows that were more divergent within the two groups, F st was calculated for each SNP, adapting the formula reported by Karlsson et al. (2007): where k is the SNP marker k, with frequency , ;

Then, a mean F st value (F stm) was calculated in 1‐Mb non‐overlapping windows, and the result was Z transformed using a derived formula provided by Rubin et al. (2010): where i is the 1‐Mb window, μF stm is the mean and σF stm is the standard deviation.

The distribution of the frequency for the normalized values is reported in Fig. S3. Values lower than 6 represent 0.2% of the empirical distribution of all the normalized values and were the most divergent between the two groups and, therefore, were considered as significant. The normalized plot is reported in Fig. 1. A total of seven windows had values that were higher than the threshold: five located on chromosome 2 (positions 66–67, 67–68, 75–76, 77–78 and 79–80 Mb), one window located on chromosome 4 (position 53–54 Mb), and one on chromosome 15 (position 105–106 Mb). The number of SNPs included in the selected windows ranged from eight to 14 (Table S2). The ANPEP gene, known to be the receptor of PEDV and a possible candidate for PEDV resistance (Nam & Lee 2010), is located on chromosome 7 at position 60 240 145–60 262 914. The window that contains this gene had a ZF stm score of 0.2, which is far below the threshold chosen for significance. Therefore this gene, which has previously been suggested to be linked to the replication rate of the virus (Nam & Lee 2010), did not appear to be related to PEDV recovery in this study. More specific analyses need to be conducted to establish if this gene does play a role in the recovery of suckling piglets after acute infection, because the number of markers and the linkage disequilibrium in the considered windows can affect the significance above all when only one gene in the window is involved.

Figure 1

Normalized F st plot. ZF st values are reported on the Y axis; on the X axis, chromosomes are labelled in different colors (as indicated by the legend). The black line across the plot indicates the fixed threshold.

After the completion of the genomic analysis, immunoglobulin G (IgG) in the serum of all animals was measured using either an indirect ELISA (Swinecheck PED; Biovet Inc.) or an indirect immunofluorescence assay (Iowa Veterinary Diagnostic Laboratory), depending on the farm (data not shown). Half (52/100) of the German pigs, including both dead (n = 11) and recovered (n = 41), tested ELISA positive. Because all dead pigs were 8 days old or younger at sampling, positive ELISA tests were indicative of maternally derived antibodies. In spite of the fact that the recovered pigs were exposed to PEDV during their first week of life, serum PEDV IgG was not measured in German pigs until 14–21 days of age. At that age, the majority of the piglets likely had naturally seroconverted (Opriessnig et al. 2014); hence, it was not possible to distinguish whether a positive response was the result of maternal antibodies or an adaptive antibody‐mediated immune response. Furthermore, it is known that the Swinecheck PEDV ELISA cross‐reacts with porcine deltacoronavirus (personal communication, A. Ambagala, June 24, 2016). Re‐analysis of the genomic data with and without the ELISA positive pigs still identified the most interesting regions on the same chromosomes (chromosomes 2: 55–56 Mb and 86–90 Mb, 4: 53–54 Mb and 15: 105–106 Mb; for more details see Fig. S4), and hence, the results obtained appear reliable. Moreover, for this specific case, the imbalance between dead and survived and across farms (Germany contributed with 39 dead piglets and nine recovered) could influence the results. As both dead and recovered pigs in these cohorts were ELISA positive, it still suggests they responded differentially to the infection during the acute phase of the disease and are useful in this genomic analysis.

The seven selected windows were then screened to find annotated genes located 200 kb before and after the windows, according to Sscrofa 10.2 and the annotation provided by Ensembl (http://www.ensembl.org/). To avoid ambiguous windows, we considered only protein‐coding genes (not pseudogenes or miRNAs) that were officially named, uniquely orthologous in mammalian species or named in the Uniprot database. A total of 152 genes fit these requirements, with only one gene located in the window on chromosome 4 and the majority (40 genes) located on chromosome 2: 74.8–76.2 Mb. The list of the annotated genes for each window is reported in Table S2. The detected genes were analyzed to investigate their role during virus infection and proliferation.

All genes were screened to find evidence of gene perturbation and differential gene expression after virus–host protein interaction using the online tool enrichr (Chen et al. 2013). This analysis was performed because it is known that host targets of viral proteins reside in networks in proximity to products of disease susceptibility genes (Goh et al. 2007; Zhong et al. 2009; Gulbahce et al. 2012). Then, enrichr and the gorilla tool (Eden et al. 2007, 2009), with Homo sapiens as a reference, were used to find Gene Ontology (GO) terms of functions, processes or components for the considered genes. Only GO terms with P ≤ 0.001 were considered, which is the minimum threshold allowed by gorilla. Several genes that have been found to be underexpressed or overexpressed in many viral infections, including coronaviruses (such as SARS‐CoV), are shown in Table S3. A total of 32 genes were included in one or more GO terms (Fig. 2). The GO term outputs are related mainly to regulation of cell proliferation, ion transport, Golgi apparatus, ATPase activity, hydrolase activities and pyrophosphatase activity. Interestingly, several of these GO terms could be directly linked to PEDV proliferation and diarrhea development, as follows: (i) Cell proliferation has been linked to the ability of the organism to overcome virus or bacterial tissue destruction (Schneider & Ayres 2008); (ii) ion channels, particularly K+ and Ca++ channels, are necessary for intestinal homeostasis, and their functions are altered during diarrhea in several species including humans (Field 2003) and pigs (Moeser & Blikslager 2007); (iii) ATPase activity is required for the function of several ion channels (reviewed by Gouaux & Mackinnon 2005); and (iv) the assembly of the PED virus occurs by budding through intra‐cytoplasmic membranes, such as the endoplasmic reticulum and Golgi apparatus (Ducatelle et al. 1982).

Figure 2

A Gene Ontology (GO) term table with the genes included in each GO term detected with gorilla (a) and enrichr (b); P‐values ≤0.001. In the plots, output for biological processes (yellow), functions (green) and components (blue) are reported.

In conclusion, we did not find evidence of a direct association between the ANPEP gene and PEDV resilience in naïve piglets. Instead, the picture provided by the F st analyses seems to support the associations between recovery and host responses that could influence cell epithelial regeneration rate, virus replication rate and the general consequences of virus infection.

Conflicts of interest

The authors declare they do not have any conflicts of interest. Data reported in this work will be available one year after the publication of the manuscript at the following link: http://www.animalgenome.org/repository/pub/ISU2016.0617/.

Figure S1 Admixture plot considering K = 8 for the analyzed data.

Click here for additional data file.

Figure S2 Multidimensional scaling plot of all dead and recovered piglets considering the first two dimensions.

Click here for additional data file.

Figure S3 Distribution of the frequency of the ZF stm values in 1‐Mb windows calculated comparing the dead with the recovered suckling piglets.

Click here for additional data file.

Figure S4 Normalized F st plot without the samples that were positive based on the ELISA test.

Click here for additional data file.

Table S1 Number of dead and recovered piglets from the PEDV outbreak in each considered farm.

Click here for additional data file.

Table S2 Information about the windows with ZF stm > 6.

Click here for additional data file.

Table S3 Virus changes, upregulated and downregulated from the Gene Expression Omnibus (GEO).

Click here for additional data file.