Genomic analyses provide insights into breed-of-origin effects from purebreds on three-way crossbred pigs

Genome sequencing data of pigs

In this study, we first downloaded genome sequencing data of 30 pig individuals (including 11 Duroc, nine Landrace and 10 Yorkshire pigs) with mean genome coverage of ~15.47× for each individual (Table S1). We then sequenced a DLY pig family, including two grandparents (a male Landrace pig and a female Yorkshire pig), two parents (a male Duroc pig and a female LY crossbred pig) and six offspring (three male and three female DLY crossbred pigs) using Illumina HiSeq 4,000 platform (Table S2). In total, we generated ~962.46 Gb paired-end 100-bp (PE100) high quality sequencing data for 10 pig individuals, with mean genome coverage of ~38.51× for each individual (Table S2).

SNP calling

The genome sequencing data of pigs in this study were first mapped to reference pig genome (v.11.1, see “URLs”) using BWA (v.0.7.8) with default options (Li & Durbin, 2009). We used ‘MarkDuplicates’ module in package Picard (v.1.48, see “URLs“) to remove duplicated reads. The module ‘HaplotypeCaller’ in Genome Analysis Toolkit (GATK; v.3.7) (McKenna et al., 2010) was used to call SNPs, which were next filtered by following criteria: QUAL < 30.0, QD < 2.0, MQ < 40.0, FS > 60.0. Sex chromosomes (X and Y) were excluded for SNP calling. To further improve the accuracy of SNP calling, we empirically depleted ~2% SNPs located in left and right tails, based on the distribution of SNP depth (Fig. S2). Finally, we identified a total of ~12.56 million (M) SNPs of three pig breeds (Table S1) as well as ~13.29 M SNPs of the DLY pig family (Table S2). We used Illumina’s porcine 60K Genotyping Bead-Chip (v.2) to validate the accuracy of SNP calling for each individual of the sequenced DLY pig family (Table S3). PLINK software (Purcell et al., 2007) was used to calculate the mendelian error to confirm the accurate kinship of the DLY pig family.

Calculation of the probability of heterozygous SNP

To calculate the probability of heterozygous SNP (PHS) of simulated offspring generated by different crossbreeding systems, we first calculated the frequency of three genotypes (homozygous identical with reference, homozygous distinct from reference, heterozygous) for each SNP locus based on the ~12.56M population-scale SNPs for Duroc (11 individuals), Landrace (nine individuals) and Yorkshire (10 individuals) respectively. We then simulated crossbreeding between each two pig breeds and recalculated the frequency of the same three genotypes for simulated offspring. Based on this method, we could obtain PHS for each SNP locus of the simulated offspring generated by different crossbreeding systems.

Principal component analysis and function enrichment analyses

We performed principal component analysis with software GCTA (v.1.91.3 beta, see “URLs”) based on the ~12.56M population-scale SNPs. Function enrichment analyses were performed using the online toolkit ‘Metascape’ with parameters ‘Min Overlap = 3, P Value Cutoff = 0.01, Min Enrichment = 1.5’ (Zhou et al., 2019).

Identification of ASE genes

For each DLY individual, RNA was extracted from adipose, heart, liver and skeletal muscle. RNA-Seq libraries were prepared from total RNA using poly(A) enrichment of the mRNA to remove ribosomal RNA (rRNA). All RNA-Seq libraries were sequenced using an Illumina HiSeq4000 platform, generating a total of ~111.52 Gb transcriptome data (Fig. S3). After quality control, high-quality transcriptome data were mapped to reference pig genome (v11.1, see “URLs”) using STAR software (v.2.6.0) (Dobin et al., 2013) with parameters “—outSAMattributes NH HI NM MD—alignEndsType EndToEnd” (Fig. S4A). The heterozygous SNPs of each DLY individual were phased via trio genotypes and only the phased exonic SNPs were used for further analyses. Reads covering multiple SNPs were trimmed to leave one SNP locus so that a read only can be counted once, a necessary step for statistical analysis. SAMtools (Li et al., 2009) was then used to generate a pileup file, which were used to calculate the total number of reads aligning to each allele (Fig. S4A). A total number of reads mapped to each gene for different alleles were summed up and we required a minimal count of 10 for further analysis (Andergassen et al., 2015; Andergassen et al., 2017). A binomial test was then used to assess the signifcance of deviation of the observed allelic biases from the expected 1:1 distribution for ASE for each sample (Fig. S4A). P-values were adjusted by bonferrioni method (Fig. S5). Each gene was then classified to three categories, including ASE (adjust-P < 0.05, binomial test), NA (not significant for ASE), no heterozygous SNP (Fig. S4A). A high-confidence ASE gene should include two ASE categories across the six DLY individuals (Fig. S4B).

URLs

Reference pig genome, http://hgdownload.soe.ucsc.edu/goldenPath/susScr11/bigZips/susScr11.fa.gz; Picard, http://picard.sourceforge.net; GCTA, http://www.gcta-ga.org.

Distribution of heterozygous SNPs in three-way crossbred pigs

To estimate the distribution of heterozygous SNPs in three-way crossbred pigs, we calculated the PHS for simulated offspring generated by three kinds of three-way crossbreeding systems, namely, DLY, Landrace × (Yorkshire × Duroc) (LYD) and Yorkshire × (Duroc × Landrace) (YDL) (Figs. 1A–1F). This was based on the ~12.56M SNPs of the publicly available genome sequencing data of 30 pig individuals from three purebreds (including 11 Duroc, nine Landrace, and 10 Yorkshire pigs) with mean genome coverage of ~15.47× for each individual (Table S1). Compared with the other two kinds of three-way crossbreeding systems, DLY significantly exhibited more high-probability heterozygous SNPs (Fig. 1G, mean Pearson’s r = −0.99, P < 0.01). This was most likely contributed to the closer genetic relationship between L and Y (Fig. 1H). As a result, compared with the other two kinds of two-way crossbreeding systems (Fig. S1A), LY significantly exhibited fewer high-probability heterozygous SNPs (Fig. S1B, mean Pearson’s r = −0.93, P < 0.01). Taking together, we suggest that DLY is a more effective strategy for three-way crossbreeding system among Duroc, Landrace, and Yorkshire breeds as it could yield more stably inherited heterozygous SNPs.

Validation of the accuracy of PHS

To validate the accuracy of PHS, we sequenced a DLY pig family, including two grandparents (a male Landrace and a female Yorkshire), two parents (a male Duroc and a female LY crossbred individual) and six offspring (three male and three female DLY crossbred individuals) with mean genome coverage of 38.51× for each individual (Table S2). A total of ~13.29M SNPs were identified, which covered ~74.11% of the SNPs identified from the three pig purebreds as well as ~98.22% of the homozygous SNPs and ~98.36% of the heterozygous SNPs identified from Illumina’s Porcine 60K Genotyping Bead-Chip (v.2) for each individual (Tables S2 and S3). As expected, the six DLY pigs exhibited a larger number of heterozygous SNPs (ranged from 6.12M to 6.26M) than the other four pigs (ranged from 4.44M to 5.62M) (Table S2). The percentage of heterozygous SNPs of each pig individual was calculated for 20 equal intervals (based on PHS of simulated DLY offspring, from 0% to 100%, at intervals of 5%) (Fig. 1I). We found that the distribution of simulated DLY offspring was significantly similar with the six DLY pigs (mean Pearson’s r = 0.98, P < 0.01), but significantly different with the other four pigs (mean Pearson’s r = −0.39, P < 0.01) (Fig. 1I). This result strongly confirmed the accuracy of PHS and supported our findings above.

Function analyses of high-probability heterozygous SNPs in DLY pigs

Annotation of the 14,848 high-probability heterozygous SNPs in DLY pigs (PHS > 0.9) using ANNOVAR software (Wang, Li & Hakonarson, 2010) and S. scrofa ENSEMBL gene annotation (v.92) revealed that more than half (~59.68%) of the loci were enriched in intergenic regions (Fig. 2A), suggesting that a large fraction of the high-probability heterozygous SNPs were acted in regulatory regions of the genome (Carneiro et al., 2014). The majority (~93.72%) of the left loci were enriched in intronic regions, remaining ~6.28% located in upstream, downstream, UTR and exonic regions (Fig. 2A). We found none of the coding SNPs was a nonsense, suggesting that gene loss may not play a major role for breed-of-origin effect on DLY pigs. Only six genes (AHI1, AKAP9, C8H4orf54, ENDOU, USP20, ZNF507) were identified to harbor a nonsynonymous SNP locus (Fig. 2B), indicating that very few loci had the power to completely alter protein-coding genes.

We further tested whether the six genes may have a close association with the crossbred performance. AHI1 is required for both cerebellar and cortical development and has been previously shown to be associated with fat development and obesity via regulating insulin signaling (Niu et al., 2012). AKAP9 is a member of structurally diverse proteins that have the common function of binding to the regulatory subunit of protein kinase A (Jo et al., 2016; Venkatesh et al., 2016). Mice with AKAP9 knockout displayed decreased body fat and body weight, hematopoietic abnormalities, and an atypical plasma chemistry profile, suggesting that this gene may be important for fat and body weight development (Gardin & White, 2011; Van Der Weyden et al., 2011).

Identification of ASE genes

To investigate the breed-of-origin effects from purebreds on gene expression pattern during the development of DLY pigs, we sequenced four representative tissues, including adipose, heart, liver and skeletal muscle (Fig. S3) and sought to identify ASE genes for each tissue. Firstly, a large fraction (ranged from ~88.59% to ~89.59%) of heterozygous SNPs for each DLY pig individual were phased via trio genotypes, that is, when at least one member of the trio was homozygous for the reference or alternate allele (Iliadis, Anastassiou & Wang, 2012; Lajugie et al., 2013) (Figs. 3A–3G). ASE genes were then identified based on these phased SNPs across the six DLY pigs for each tissue (Figs. S4A and S4B). In total, we identified 95, 79, 132 and 42 ASE genes in adipose, heart, liver and skeletal muscle, respectively, with a median value of allelic ratio ranged from 0.68 to 0.82 (Fig. 4A; File S2). Notably, we observed that ~78.09% (196/251) of ASE genes were only detected in one tissue (Fig. 4B), indicating that ASE genes were mainly tissue-specific, consistent with the results of principal component analysis (Fig. 4C) and unrestricted cluster analysis (Fig. 4D) based on gene expression of the 251 nonredundant ASE genes.