Genomic divergence between populations

The results of the MDS and the estimation of individual ancestries with maximum likelihood are presented in Figures 1 and 2, respectively. Both approaches indicate that the genomic data clearly identified the genetic structure of the populations. The representation of the MDS with the first two dimensions distributed the individuals in six groups that correspond to the pure (EE, RR and TT) and crossbred (ER, ET and RT) populations. Graphically, the crossbred were located between the purebred populations that generated the cross. In addition, the estimation of individual ancestries after the assumption of three populations distributed the purebred individuals into three different ancestries and assigned crossbred individuals with approximately half ancestry from each parental population.

Average ± standard deviation of the F_ST results between the three varieties of Iberian pig was 0.069 ± 0.060, which is lower than the estimate reported by Fabuel et al. (2004), who used a set of 36 microsatellites (F_ST=0.129). However, they are not directly comparable since microsatellites have a higher mutation rate and, therefore, F_ST estimates are not in the same scale. In addition, Faubel et al. (2004) included five varieties (Entrepelado, Retinto, Lampiño, Torbiscal, and Guadyerbas), and the Guadyerbas population had the greatest genetic distance from the other populations. Therefore, a greater estimate of F_ST was expected in their analysis. We did not include Guadyerbas because it is almost entirely restricted to conservation programs. The mean ± standard deviation of F_ST statistics between pairs of populations was 0.045 ± 0.055 between EE and RR, 0.049 ± 0.059 between EE and TT, and 0.057 ± 0.072 between RR and TT. The results of paired t-tests were highly significant (p<1e-8) for all comparisons. The results confirmed the closeness between Entrepelado and Retinto, as it was previously reported by Fabuel et al. (2004).

The results of the genomic scans through the porcine autosomal genome of the single SNP F_ST statistic for all populations and after averaging them in sliding windows of 1, 2, and 3 Mb centered at each SNP are presented in Figure 3. The distribution along the genome of F_ST statistics calculated with a single SNP was extremely noisy. Therefore, it was not possible to extract a clear pattern of the genomic differentiation between populations and confirmed the need to average estimates of the F_ST in wider genomic regions. The number of markers included in each window was 18.37 ± 6.63, 33.77 ± 11.36, and 48.97 ± 15.68 SNP for sliding windows of 1, 2, and 3 Mb, respectively. The results from the genomic scan with sliding windows of 1 Mb were somewhat noisy, and the results based on sliding windows of 2 Mb and 3 Mb were very similar. Therefore, to achieve a compromise between noise reduction and the genomic size of the windows, we decided to explore the results based on 2 Mb sliding windows.

The genomic scans for each pair of populations (EE and RR, EE and TT and RR and TT) with sliding windows of 2 Mb are presented in Figure 4. The distributions of the F_ST estimates along the autosomal chromosomes were clearly different among the three genomic scans. However, we were able to detect some degree of similarity since the correlations between the F_ST estimates obtained from EE and RR and EE and TT, EE and RR and RR and TT and EE and TT and RR and TT were 0.218 (p<0.001), 0.214 (p<0.001), 0.317 (p<0.001), respectively.

Biological processes and candidate genes between Entrepelado and Retinto

The genomic regions that had an average F_ST(ER) greater than the 95^th percentile (0.082) contained 651 genes, which were used in an ORA for the GO for biological process. Among them, 569 and 157 genes were annotated to functional categories in the Homo sapiens and Sus scrofa databases, respectively. The enriched GO terms that had a p-value < 0.0001 are presented in Table 1. We identified up to 15 and 3 terms within the human and the porcine databases, respectively. All of the GO terms identified with the Homo sapiens reference (anterior/posterior pattern specification, skeletal system development, limb morphogenesis, appendage morphogenesis, pattern specification process, regionalization, embryonic skeletal system development, appendage development, limb development, skeletal system morphogenesis, embryo development, embryonic limb morphogenesis and embryonic appendage morphogenesis) and two identified with the Sus scrofa (skeletal system development and anatomical structure morphogenesis) were associated with embryogenesis and morphogenesis. In addition, a biological process associated with regulation of gene expression was significant with the Sus scrofa database.

Table 1. Enriched gene ontology (GO) terms for biological processes (p < 0.0001) with genes located within the genomic regions that have an average F_ST over the 95^th percentile between Entrepelado and Retinto Iberian pig populations based on the Homo sapiens and Sus scrofa databases.

Those results are reinforced by the genes within the genomic regions greater than the 99.9^th percentile (0.157) which were located at SSC8 (56983392-60232132 bp), SSC15 (80515676-86990138 bp) and SSC17 (42130987-42480736 bp) as presented in the Table S1 [suppl.]. Among them is a family of HOXD (Homeobox protein) genes that encode a family of transcription factors that play a crucial role in morphogenesis (Myers, 2008), jointly, with the tightly linked EVX2 (Even-Skipped Homeobox 2) (Hérault et al., 1997), and the SP9 (Sp9 Transcription Factor) (Kawakami et al., 2004). In addition, the NEUROD1 (Neurogenic differentiation 1) is a transcription factor involved in regulatory networks in embryonic stem cells (Marchand et al., 2009), the OLA1 (Obg Like ATPase 1) gene plays a role on the attachment of cells to the extracellular matrix (Jeyabal et al., 2014), and the ATF2 (Activating transcription factor-2) which has been found to affect skeletal growth (Vale-Cruz et al., 2008). Some other interesting genes located within them that are related with morphogenesis are the CHN1 (Chimerin 1) and the PRKRA (protein kinase, interferon inducible double stranded RNA dependent activator). The CHN1 is mostly expressed in the brain and it is associated with the early development of the nervous system (Lim et al., 1992) and PRKRA has been related with the development of the cerebellum (Yong et al., 2015). An additional evidence of the relationship of those genomic regions with the embryological development is that, in pigs, they have been associated with QTL for teat numbers in SSC8 (Velardo et al., 2016), number of mummies (Onteru et al., 2012) and stillbirths (Schneider 5, 2015) in SSC15.

Moreover, some of those genes such as the HOXD family are also associated with the regulation of gene expression, but it is noteworthy that the ACRT5 (ARP5 Actin Related Protein 5 Homolog) is involved in the INO80 complex that contributes to transcription, DNA repair, and DNA replication (Conaway & Conaway, 2009). Therefore, they can be related with the biological process last identified (regulation of gene expression) in the Sus scrofa database.

Biological processes and candidate genes between Entrepelado and Retinto

The genomic regions that were within the 95^th percentile (0.092) of the average F_ST(ET) in sliding windows of 2 Mb contained 886 genes, of which 784 and 228 were annotated in the Homo Sapiens and Sus scrofa databases, respectively. The results of the ORA with these genes are presented in Table 2. One of the enriched GO terms for biological processes was also associated with morphogenesis (animal organ morphogenesis) and the other terms were associated with the global metabolism of the individuals through metabolic or catabolic processes (collagen catabolic process, multicellular organismal catabolic process, and multicellular organism metabolic process), regulation of hormone secretion (negative regulation of hormone secretion), or the ubiquitination of proteins (positive regulation of ubiquitin-protein transferase activity). Furthermore, six genomic regions were detected that had an average F_ST(ET) greater than the 99.9^th percentile (0.178) and were at SSC2 (60177754-61303696, 71588075-71588075 and 92958395-93175114 bp), SSC6 (104376948-105459825 bp), SSC10 (36887963-37186755 bp), and SSC14 (45509383-45509383 bp). A search of the porcine QTL database (https://www.animalgenome.org/QTLdb/pig/) indicated that these regions have been associated with copying behavior (Ponsuksili et al., 2015) and hemoglobin content (Zhang et al., 2014) in SSC2, intramuscular fat content (Cepica et al., 2012) and vertebra number (Rohrer et al., 2015) in SSC6 and the ratio to non-productive days in SSC14 (Onteru et al., 2011). In addition, some of the genes within or near those genomic regions (see Table S2 [suppl.]) are associated with the same biological processes highlighted above, and are good candidates for had being affected by selection or adaptation; e.g., INSR (insulin receptor) or ADCYAP1 (Adenylate Cyclase Activating Polypeptide 1). The INSR gene has been proposed as a candidate gene for intramuscular fat content by Cepica et al. (2012) while ADCYAP1 is a member of the glucagon superfamily of hormones that are involved in in growth, metabolism, and immune response (Sherwood et al., 2000). In addition, NWD1 (NACHT and WD Repeat Domain Containing 1) modulates androgen receptor signaling (Correa et al., 2014), and MN1 (Transcriptional activator MN1) is involved in the development of craniofacial traits (Pallares et al., 2015). Further, another interesting gene is the CD209 (dendritic cell-specific intercellular adhesion molecule-3-grabbing non-integrin, DC-SIGN), that functions as an important pattern recognition receptor (PRR) in immune defense and plays a role in the immune modulation during pathogen infection (Soilleux et al., 2002).

Table 2. Enriched gene ontology (GO) terms for biological processes (p < 0.0001) with genes located within the genomic regions that have an average F_ST over the 95^th percentile between Entrepelado and Torbiscal Iberian pig populations based on the Homo sapiens and Sus scrofa databases.

Biological processes and candidate genes between Retinto and Torbiscal

The genomic regions that had an average F_ST(RT) within the 95^th percentile (0.101) in sliding windows of 2 Mb contained 596 genes, of which 530 and 146 were annotated in Homo sapiens and Sus scrofa databases, respectively. The results of the ORA analyses identified statistically significant (p<0.0001) GO terms when crossed with the Homo sapiens database, only (Table 3). There were up to four GO terms associated with development of the neuronal system (generation of neurons, neurogenesis, neuron differentiation and neuron projection development), and another term was associated with the general development of the organism (regulation of multicellular organismal development). Furthermore, two GO terms were associated with the responses of the organism to stress (regulation of response to stress) and lipopolysaccharides (cellular response to lipopolysaccharide), and another was associated with microtubule activity (regulation of microtubule motor activity). Some of the genes within the genomic regions that had an average F_ST that was in the 99.9^th percentile (0.189) confirmed these results (see Table S3 [suppl.]). They were located at SSC1 (77056462-77472299, 90748703-90896710, and 141989498-142042139 bp), SSC6 (104376948-105459825 bp), SSC8 (89447995-89679243 bp), and SSC12 (28731262-28831639 bp). The genomic regions of SSC6 was also identified within the most divergent genomic regions between EE and TT, which strengthens the argument that, among other, INSR and AD-CYAP1 genes may be good candidates to have been affected by selection or adaptation processes that have influenced the genetic configuration of the TT population. Moreover, other genes of note include FYN (FYN Proto-Oncogene, Src Family Tyrosine Kinase), which is involved in the early stages of neurogenesis (Yagi et al., 1994), CDK19 (Cyclin Dependent Kinase 19), a regulator of the p53 network for cellular response to stress (Audetat et al., 2017), and TRAF3IP2 (TRAF3 interacting protein), which plays a central role in innate immunity in response to pathogens (Wu et al., 2013). In addition, COL12A1 (Collagen, type XII, alpha-1) is involved in bone formation (Izu et al., 2011) and CA10 (Carbonic anhydrase-related protein CA10) has been identified as an important neurexin ligand (Sterky et al., 2017). Finally, it must be remarked that the genomic regions of SSC1 (77056462-77472299) were associated by Fontanesi et al. (2017) with a large QTL for intramuscular fat and that the genomic regions of SSC6 and SSC12 has been related with the genetic variation in vertebra numbers (Rohrer et al., 2015).

Table 3. Enriched gene ontology (GO) terms for biological processes (p < 0.0001) with genes located within the genomic regions that have an average F_ST over the 95^th percentile between Retinto and Torbiscal Iberian pig populations based on the Homo sapiens and Sus scrofa databases.

Final remarks

The main conclusion of this study is that the processes of differentiation among Iberian pig varieties have heterogeneously affected the autosomal genome. A first approximation has identified potential candidate genes, most of which are associated with morphogenesis, neuronal development, regulation of metabolism, or cellular response to stressors. The presence of those candidate genes is coherent with the recent evolution of the Iberian pig populations, they have evolved through adaptation to harsh environmental conditions. Furthermore, producers have subjected the Iberian pig to “empirical” selection in which adipogenic capacity and morphological traits have played an important role. Nevertheless, further research must be done to confirm these results.