Detection of genomic regions associated malformations in newborn piglets: a machine-learning approach

Introduction

The mortality of piglets represents a severe loss for the pig farming industry (Pig Progress, 2014; Walters, 2016). The losses predominantly occur during the last days of pregnancy, delivery, and the first three days of life (Woollen, 1993). Many factors contribute to the loss of piglets: the age of sows and their maternal qualities, the birth duration, the litter size, the parity number, farm management, and other factors (Olsson, Botermans & Englund, 2018). Pigs have more frequent congenital malformations than any other domestic animal (Woollen, 1993), and congenital malformations account for a significant proportion of perinatal losses.

The profitable production of pigs requires that a farm consistently produces healthy piglets (Staarvik et al., 2019). However, the piglet mortality rate is positively correlated with the number of piglets born (Olsson, Botermans & Englund, 2018), making it a limiting factor for farm growth. Therefore, identifying and addressing the causes of congenital malformations in pigs is essential for pig farming.

Congenital diseases often pose a diagnostic problem, having unknown etiology, pathology, and prognosis. This is also true beyond animal husbandry: up to 79% of human congenital diseases have unknown etiology (Papatsiros & Others, 2012).

The frequency of occurrence of specific malformations varies from farm to farm. Breeds with small populations have smaller genetic diversity and are prone to more defects. The most common defect, affecting between 2.3% and 5.0% of all newborn piglets, is the piglet splay leg (PSL) (Papatsiros & Others, 2012). PSL is a condition when a newborn piglet cannot hold its hindlimbs together (in severe cases, the forelimbs are also affected) and cannot stand or walk after birth (Partlow et al., 1993) and may, therefore, starve. There are also cases of inguinal, scrotal, and umbilical hernia, cryptorchidism, and atresia. The prevalence of this condition is between 2.3% and 5.0%.

Genetic factors are believed to play an essential role in the development of congenital malformations in piglets. Studies found genetic drivers of PSL (Hao et al., 2017) and cleft palate (which is a lethal anomaly in pigs) describing chromosomal imbalances in affected offspring (Grahofer et al., 2019). All affected offspring were carriers of partial trisomy of chromosome 14, including the FGFR2 gene associated with various dominant hereditary syndromes of craniofacial dysostosis in humans. Hao et al. (2017) presented seven significant SNPs and four genes related to muscle development, glycogen metabolism, and mitochondrial dynamics, identified as potential candidate genes for PSL. However, only limited research aimed at finding SNPs and candidate genes associated with congenital malformations in pigs has been conducted (Bermingham et al., 2014; Ji et al., 2016; Bakoev et al., 2020).

When analyzing extensive genomic data with traditional parametric models, the problem of a small number of observations and a large number of predictive variables (typically, tens or even hundreds of thousands of SNPs are interrogated in thousands or, sometimes, hundreds of samples) is challenging to solve due to the high dimension or highly correlated data structure (Li et al., 2018). Traditional GWAS analysis focuses on a univariate hypothesis and assumes the presence of independent explanatory variables. These methods suffer significantly from a lack of power and accuracy when dealing with the complexity of multiple interactions or correlations between predictors (e.g., SNP-SNP and SNP- covariate interactions) (Lettre, Lange & Hirschhorn, 2007; Zheng et al., 2007; So & Sham, 2011; Adams, Bello & Dumancas, 2015). Machine learning approaches, and Gradient Boosting Machines (GBM), in particular, are suitable for such complex problems (Li et al., 2018). Various machine learning applications were developed to build a nonparametric regression or classification model from data. One strategy is to build a model from theoretical considerations and then define its parameters based on the observed data. The main idea behind the GBM algorithm is to create new baseline learners that are maximally correlated with the negative gradient of the loss function associated with the entire ensemble. There is a wide choice of existing loss functions that can be adjusted depending on the specific task at hand. This flexibility makes GBMs highly customizable for any specific data-driven task. It introduces much freedom in model design, making it a matter of patience since it takes trial and error to select the most appropriate loss function.

We have conducted a pilot GWAS to identify potential genomic regions and genes associated with congenital diseases and anomalies in newborn piglets. We used two common pig breeds for our study: Landrace and Large White. The Landrace pigs originated from the free-breeding, non-pedigreed stock of swine, the regional landrace native to Denmark. The Large White pig was produced in the 18th century by crossing the large indigenous white pig of North England with the smaller, fatter, white Chinese pig. These breeds have different genetics but share the problems with congenital malformations. We use the gradient boosting machine (GBM) approach since it is flexible and relatively easy to implement. Moreover, GBM has demonstrated considerable success in practical applications and various machine learning and data extraction tasks (Bissacco, Yang & Soatto, 2007; Hutchinson, Liu & Dietterich, 2011; Pittman & Brown, 2011; Johnson & Zhang, 2014). We have evaluated the applicability of GSM for searching for genomic areas and genes associated with congenital anomalies and compared its performance to standard approaches.

Materials & Methods

Results

We have conducted admixture analysis of both breeds for K = 2...11. Both breeds appear to be genetically diverse. Both Admixture and PCA (Fig. S1) plots showed that LA and LW pigs belong to two different sample groups. However, for all K (Fig. 1), the admixture profiles of Large White pigs are a mosaic of almost all components, while a smaller number of admixture components represents the Landraces. Therefore, the sampled Landrace pigs were genetically more homogeneous than the Large White pigs. To determine the optimal number of components, we have analyzed cross-validation error and accuracy of clustering.

As a function of the number of admixture components, the cross-validation error drops to a plateau value at K = 4. Next, we used NbClust (Charrad M., Ghazzali N., Boiteau V., Niknafs A., 2014) kmeans procedure in R to partition the samples using 26 different scores to define the optimal number of clusters. The number of clusters for K>3 was between 6 and 8.

GPS (Elhaik et al., 2014) leave-one-out analysis validated the clustering assignment correctness. For K = 2 and 3, nearly 10% of the samples were misclassified (the rest was placed to correct clusters). For the higher values of K, GPS was placing 100% of samples into correct clusters. The fraction of affected animals per cluster ranged from 11% to 42% for splay leg phenotype and 13% to 83% for other congenital anomalies. For example, for K = 11, in cluster LA_2, all six piglets were affected; one had a splay leg and five other congenital anomalies. In cluster LW_3, having 23 members, three piglets were born with splay legs, and no other anomalies were detected.

We have calculated Pearson’s correlation coefficient between piglets’ status and values of admixture components. For example, for K = 18, in both breeds, the presence of congenital anomalies is negatively correlated with the admixture component K8 (LW: r =−0.26, t-test p-value =0.0002). This effect is especially pronounced in Large White piglets: 3% of piglets with K8>0.3 (1 out of 38) have congenital defects, as compared to 31% of piglets with K8<0.3 (29 out of 94). Cluster LW_7, where K8 is the main ancestry, is the healthiest. Congenital anomalies in Large White piglets are positively correlated with the admixture component K5 (r = 0.21, t-test p-value = 0.007). 41% of piglets with K5>0.3 (15 out of 37) have congenital defects, as compared to 16% of piglets with K5<0.3 (15 out of 95). Therefore, it appears that certain pig lines are more susceptible to the development of congenital anomalies.

F_st between Landrace and Large White pigs was estimated by the method of Weir & Cockerham (1984) to be 0.165. Hierarchical clustering (Fig. 2A), PCA analysis (Fig. 2B) showed the separation of Landrace and Large White pigs into two distinct groups. Note that the PCA plot had the arch and horseshoe effects for Large White pigs. These effects are common in ecological and population studies when the sampling units cover a long ecological gradient. Those samples at each end of the gradient have few species in common. These artifacts produce visually beautiful but false patterns that are not present in the original data. Using the SNPRelate package in R (Zheng et al., 2012), we have calculated relatedness (Fig. 2C) and IBS (Fig. 2D). Within- and between-breed relatedness have been inferred by the IBD using the Plink method of moments (Purcell et al., 2007). We have calculated the Cotterman coefficients of relatedness K0 and K1 (Cotterman, 1940). These plots show that Large White samples mainly were from unrelated individuals, while the Landrace samples range from unrelated to first-degree relatives. Clustering (A) and IBS (D) plots also showed population stratification within the Large White pigs: the samples could be separated into two main clusters.

Next, we used the rrBLUP (Endelman, 2011) R package and Plink (Chang et al., 2015) to find significant SNPs associated with the splay leg and congenital malformations phenotypes separately for the Large White and the Landrace breeds. Plink performs five tests: basic Allelic test, Cochran-Armitage trend test, Genotypic test, Dominant gene action test, and Recessive gene action test.

Using the p-value cut-off of 10⁻⁵, there were no significant SNPs in Landrace pigs identified by Plink (p-value <10⁻⁵ in at least one of the five Plink tests) and also by rrBLUP (Fig. 3). In Large White pigs, one SNP (MARC0062453 a.k.a. rs80912121) was found significant by two Plink tests for piglet splay leg (its rrBLUP p-value = 1.334 ×10⁻⁵, hence, it nearly misses the cut-off). It is an intron variant for CRISP3 (cysteine-rich secretory protein 3) on chromosome 7. For other congenital anomalies of Large White pigs, rrBLUP has identified two SNPs, intergenic variant ALGA0086405 and downstream gene variant WU_10.2_13_216912741, the first of each was significant in four out of five tests, according to Plink. SNP ALGA0086405 (rs81454026) is located on Chromosome 15 position 92,733,414. There are three genes near this SNP: GULP, TFPI, and CALCRL. Gene GULP PTB (domain-containing engulfment adaptor 1) is located on Chromosome 15: 92,838,715–93,251,019. Gene TFPI (tissue factor pathway inhibitor) is located on Chromosome 15: 92,345,075–92,409,905. Gene CALCRL (calcitonin receptor-like receptor) is located on Chromosome 15: 92,233,841–92,332,363. These genes may be associated with the disease since, due to artificial selection and small effective population size in pigs, the LD in pigs is much larger than in humans (Badke et al., 2012).

SNP WU_10.2_13_216912741 (located downstream of a lncRNA located on Chromosome 13: 207,036,621-207,042,808 reverse strand) was not significant according to Plink. Plink has identified four SNPs (ALGA0097461, ASGA0009945, WU_10.2_2_80100147, and WU_10.2_6_139448200), that were not significant according to rrBLUP. ALGA0097461 is in an intron of a lncRNA on Chromosome 18: 23,230,190–23,234,168, forward strand. ASGA0009945 is an intron variant of gene growth arrest-specific 2 gene GAS2 on chromosome 2: 36,713,106–36,840,407 reverse strand. WU_10.2_2_80100147 is an intron variant of TBC1 domain family member 9B, located on Chromosome 2:78548370 (forward strand). WU_10.2_6_139448200 is an intergenic variant on Chromosome 6:151776796 (forward strand).

Next, we used the GBM approach (see Table 2, Table 3, Table S5, and Table S6. A complete list of GBM SNPs is in Table S7). GBM has identified 205 variants combined from separate analyses of two breeds and two conditions. The Ensembl Variant Effect Prediction tool classified the SNPs into the following categories: intron variant: 71%; non coding transcript variant: 9%; intergenic variant: 8%; downstream gene variant: 5%; upstream gene variant: 4%; 3-prime UTR variant: 2%; non coding transcript exon variant: 1%; synonymous variant: 1%. To compare the GBM and Plink performance, we have calculated the confusion matrices (Table S8). GBM was more accurate than Plink in labeling the affected piglets as affected, but this came at the cost of a higher false-positive rate. After training and testing the models, the primary performance statistics were calculated (Table 2).

GBM has found 49 SNPs (23 SNPs in LA pigs and 26 SNPs in LW) linked with congenital anomalies. Twenty-two of these SNPs were located in genes (Table 3). In total, 156 significant SNPs were identified for the splay leg condition. In Landrace pigs, twenty-eight SNPs were identified; fifteen were located in genes (Table S5). In Large White pigs, 128 SNPs were identified (63 of the SNPs were in fifteen genes, Table S6).

The overlap between GBM, Plink (selecting SNPs significant in at least one of the five Plink tests), and rrBLUP identified SNPs is shown in Figs. S2 and S3. Only five SNPs were detected by all three methods: ALGA0075476, MARC0030657, WU_10.2_6_138715609, CASI0008977, and ALGA0021723 (Fig. 3).

ALGA0075476, significant in Landraces for splay leg, is on the 14, position 12790416 in the intron of homeobox-containing gene 1, HMBOX1. According to the EBI GWAS catalog, SNPs in this gene’s human ortholog are associated with obesity-related traits (https://www.ebi.ac.uk/gwas/variants/snp:rs2221894).

MARC0030657 (also a significant marker in Landraces for the splay leg) is on chromosome 1, position 112210717. This SNP is located in the Intron 5–6 of lncRNA ENSSSCG00000048340. In general, lncRNAs serve as molecular signals to regulate transcription.

WU_10.2_6_138715609 (significant in Large White pigs for the splay leg) is an intergenic variant on chromosome 6, position 150780123. Two more intergenic variants (CASI0008977 at chr2:107190210 and ALGA0021723 chr3:129009157) were identified as significant in Landrace pigs for the splay leg.

SNP WU_10.2_18_48204887 was identified as significant by GBM for the splay leg in Large White pigs, for other congenital anomalies in Landrace pigs by Plink and rrBLUP. WU_10.2_18_48204887 (on chromosome 18, position 43798925) is in the intron of gene TRIL (TLA4 interactor with leucine-rich repeats). In humans, SNPs in this region are associated with anomalies in the urogenital system.

The reason for the small number of SNPs in the overlap between the three methods is in the approaches used by the programs: the first two programs use parametric methods (negatively affected by small sample size and a large number of SNPs), while in the GBM algorithm, a nonparametric approach is taken.

In Landraces, one of the identified genes (ACSL3, Acyl-CoA synthetase long-chain family member 3) encoded a metabolite interconversion enzyme from the long-chain fatty acid coenzymes A-ligases. All members of this family convert free long-chain acids into esters but differ in tissue-specific expression. Complex fatty acyl-COA esters are substrates for numerous fatty acid metabolic pathways, including mitochondrial β-oxidation and phospholipids, and triacylglycerol synthesis. Changes in ACSL enzymes can contribute to various physiological processes, including apoptosis, gene transcription, signaling, and vulnerability to oxidative stress (Johnson et al., 2012).

SNP rs344559304 in BTB16 gene (zinc finger and BTB domain-containing 16 as PLZF and ZFP145) received the top score in PSL in Large While pigs. The protein it encodes belongs to the Krüppel-like family of zinc finger proteins. It has been previously shown that ZBTB16 acts as an essential transcriptional repressor or activator (Kolesnichenko & Vogt, 2011) and participates in diverse biological processes (Dick & Doulatov, 2009), spermatogenesis (Buaas et al., 2004; Costoya et al., 2004), the formation of hind limbs (Barna et al., 2000), and the apoptosis of cells (Cheung et al., 2010). Impressive results were obtained in studies of the effect of inactivated ZBTB16 on the formation of the skeletal structures of the leg (Wei et al., 2018): ZFP145 -/- mice showed morphological defects affecting the hindlimb’s skeletal structures, whereas the forelimb’s skeletal defects occurred with lower frequency. Also, ZBTB16 was shown to regulate the expression of the HOX and BMP genes (Wei et al., 2018), influencing leg defect formation.

For PSL in LA pigs, the top SNP was rs81225364, inside an extended non-coding RNA region. LncRNAs participate in many biological processes, such as the regulation of epigenetic modifications (Roberts, Morris & Weinberg, 2014), embryo development (Zhang et al., 2014), and the development of skeletal muscle (Shi et al., 2020).

We have identified SNP rs327546322, located in the UBAC2 gene (ubiquitin-associated domain-containing protein 2), to be of the highest significance for the congenital anomalies in Landrace pigs. Polymorphisms in the UBAC2 gene are associated with a genetic predisposition to Behcet’s disease in humans (Sawalha et al., 2011; Yamazoe et al., 2017), a rare disorder causing blood vessel inflammation throughout the body: in joints, vascular system, lungs, gastrointestinal tract, central nervous system, and epididymis (Kaklamani, Vaiopoulos & Kaklamanis, 1998). The symptoms include recurrent eye inflammation, oral and genital sores, and skin lesions. No information on polymorphism of the UBAC2 gene in pigs is available in the literature; however, based on the association of the gene with Behcet’s disease and its manifestation in humans, the UBAC2 gene is a promising candidate for further investigation of congenital anomalies in pigs.

For congenital anomalies in LW pigs, the top-scoring SNP is rs81374120 in the AUTS2 gene. AUTS2 is one of 10 genes with the highest intron RNA score in the fetal brain (Sultana et al., 2002; Bedogni et al., 2010; Oksenberg & Ahituv, 2013). (Ameur et al., 2011) showed that alternative splicing programs control the expression of intronic RNAs in the fetal brain. In recent years, the role of the AUTS2 gene in developing autism spectrum disorders has been extensively studied in humans. Also, variants of this gene are associated with brain malformations and congenital disabilities (Beunders et al., 2013). Therefore, it can be hypothesized that the identified SNP rs81374120 in the AUTS2 gene in pigs is associated with nervous system development disorders and leads to congenital anomalies in piglets.

Discussion

Currently, assessing genetic load is one of the primary problems in pig farming. Various anomalies affect pig breeding’s economic efficiency, and the search for genetic factors associated with their manifestation is an urgent task for researchers. It is vital to simultaneously reduce the number of carriers of congenital defects in the Sus scrofa population and improve productivity indicators. It is not possible to solve these problems using classical breeding without using modern genomic selection approaches.

The traditional strategy of culling piglets with phenotypic defects is effective in only a small part of hereditary pathologies caused by dominant or sex-linked genes. However, most Sus scrofa genetic diseases have a recessive or polygenic type of inheritance; therefore, these culling methods are insufficient.

Therefore, the task of finding new methods to identify animals that carry hereditary pathologies is so essential. Nowadays, congenital anomalies cause significant economic losses in the pork industry. The development of genomic technology has given new meaning to the search for solutions.

When analyzing large genomic datasets using conventional parametric models, small observations and large numbers of predictive variables present an obstacle. This problem is challenging to the high dimensionality or strongly correlated data structure (Chen & Ishwaran, 2012). GWAS analysis, based on traditional statistical methods, focuses on a one-dimensional hypothesis and assumes independent explanatory variables; these methods suffer significantly from a lack of power and accuracy when dealing with the complexity of multiple interactions or correlations between predictors (e.g., SNP–SNP and SNP–covariate interactions) (Lettre, Lange & Hirschhorn, 2007; Zheng et al., 2007; So & Sham, 2011; Chen & Ishwaran, 2012; Adams, Bello & Dumancas, 2015).

The etiology and pathogenesis of the splay leg are complex and still poorly understood. Histomorphological studies have described PSL as myofibrillar hypoplasia, but this condition has also been detected in clinically normal piglets (Ducatelle et al., 1986). Ultrastructural analysis has shown that piglets with PSL have higher muscle glycogen content than normal piglets (Antalíková, Horák & Matolín, 1996). Maak et al. (2009) compared gene expression in posterior leg muscles (M. adductores, M. gracilis, M. sartorius) from diseased and healthy piglets. As a result, they proposed four genes (SQSTM1, SSRP1, DDIT4, and MAF) related to splay legs in piglets. Hao et al. (2017), in a study on 185 animals from four populations (Yorkshire, Duroc, Landrace, and crossbred Yorkshire and Landrace), identified four genes (HOMER1 and JMY on SSC2, ITGA1 on SSC16, and RAB32 on SSC1) associated with muscle development, metabolism, and mitochondrial dynamics that were proposed as candidate genes for PSL.

Our work has identified 79 genes that could be considered candidate genes for the splay leg. Even though the genes identified by Hao et al. (2017) were not identified in our work (neither by plink or rrBLUP analysis), it may be noted that in LW pigs, seven of the identified genes (EXT2, NDUFA10, MTMR12, HS6ST3, DBT, ACOXl, CHST11) belong to the protein class “metabolite interconversion enzyme” (PC00262). Metabolic changes have been documented in multifactorial diseases in humans (DeBerardinis & Thompson, 2012). Metabolism is carried out through enzymatically catalyzed biochemical reactions and involves the mutual conversion of small molecules (metabolites) that play a crucial role in various cellular functions, from energy production to complex biosynthesis macromolecules (Pey et al., 2013). The EXT2 gene (Exostosin Glycosyltransferase 2) encodes one of the glycosyltransferases involved in the chain extension stage of heparan sulfate biosynthesis. EXT2 gene polymorphism is associated with humans with exostosis (or osteoma, benign growth of new bone on top of the existing bone) (Wuyts et al., 1998).

Conclusions

Our pilot study showed that machine learning approaches could identify genomic loci associated with congenital malformations in piglets. We have identified lists of SNPs and candidate genes associated with congenital anomalies and piglet splay leg for Large White and Landrace pigs. We have also highlighted the need to study farm animals to determine population characteristics and identify genotypes associated with significant selection indicators of malformations and productive qualities.

The small sample size used here is a limitation of our pilot study. We will extend this approach to larger datasets since it is vital for ensuring efficient and ethical handling of farm animals.

Supplemental Information