The datasets analyzed during the current study are available from the NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra/) under project PRJEB1683 (Groenen et al., 2012), PRJEB29465 (Grahofer et al., 2019), PRJEB9922 (Frantz et al., 2015), PRJNA186497 (Li et al., 2013), PRJNA213179 (Ai et al., 2015), PRJNA231897, PRJNA238851 (Wang et al., 2015), PRJNA254936, PRJNA255085 (Ramirez et al., 2015), PRJNA260763 (Choi et al., 2015), PRJNA273907, PRJNA305081, PRJNA305975, PRJNA309108 (Li et al., 2017), PRJNA314580, PRJNA320525 (Bianco et al., 2015), PRJNA320526, PRJNA320527, PRJNA322309, PRJNA369600, PRJNA378496 (Zhao et al., 2018), PRJNA398176 (Zhu et al., 2017), PRJNA438040, PRJNA488327 (Yan et al., 2018), PRJNA488960 (Zhang Y. et al., 2020), PRJNA524263 (Zhang W. et al., 2020), and PRJNA550237 (Chen et al., 2020).

A total of 730 samples were included with Asian, Western, domestic and wild backgrounds (See Table S1). Raw reads were aligned to the Sscrofa11.1 reference genome (Warr et al., 2020) using the bwa-mem algorithm (Li and Durbin, 2009). Samtools-v1.8 (Li, 2011) was used for sorting, merging, and marking potential PCR duplications. Finally, haplotype-based variant detection was conducted with freeBayes-v1.1 (--min-base-quality 10 --min-mapping-quality 20 --min-alternate-fraction 0.2 --haplotype-length 0 --pooled-continuous--ploidy 2 --min-alternate-count 2) (Garrison and Marth, 2012). After SNP calling, SNP loci were screened and retaining with a quality value greater than 20 (vcffilter -f “QUAL> 20”). Further quality control was conducted with the following criteria: minor allele frequency (MAF) > 0.01, missing rate <0.01, call rate >90%, sequencing depth of sample >4. Individuals and loci satisfying the above criteria were retained for futher analyses, and assigned to their specific background (193 Chinese indigenous pigs, 30 Asian wild boars, 13 Yucatan mini-pigs, 40 Western wild boars, 298 Western commercial pigs and 18 wild suidae; Supplementary Table S1). Finally, phasing and imputation were performed based on this data set with Beagle 5.1 (Browning et al., 2018; Browning et al., 2021) (window = 20 overlap = 4 gp = true ap = true).

T-SNE dimensionality reduction was first conducted by sklearn. manifold.TSNE (n_components = 2, perplexity = 24) in scikit-learn-0.23.1 python package on high-quality Sus. Scrofa samples. To construct the Neighbor-joining tree (NJ-tree), we calculated the IBS distance matrix by plink-1.9 on phased SNP data with default parameters. Then the NJ-tree was constructed by fastME-v2.15 (-D 1 -m N -b 10000 -T 10 -s -I) (Lefort et al., 2015) with Sus cebifrons as the outgroup. The tree was plotted by the iTOl-v5 online tool (Letunic and Bork, 2021). Model-based global ancestry estimation was conducted with Admixture-1.3 (-B10 -c10) (Alexander et al., 2009) with cross-validation to assess the best fitting K-value.

Western wild boars and Yucatan minipigs were combined as the Western haplotypes background for the introgression study. Chinese groups were set as the donor population for every commercial line. Then putative introgression segments were detected for every donor-recipient combination with relative Identity-by-descent (rIBD) method using whole-genome sequencing data (Bosse et al., 2014). Identity-by-descent (IBD) detection was performed with the refinedIBD algorithm (length = 0.1 trim = 0.01 lod = 1) (Browning and Browning, 2013). These parameters were adjusted to detect not only segments that are identical, but segments with similar origins (i.e., Western or Asian) that show higher similarity than expected between Chinese and European ancestries.

The rIBD values were calculated on non-overlapping bins of 10 kb along the genome. For every bin, we calculated rIBD values with the following formula: r I B D = n I B D R , D − n I B D R , B . n I B D R , D denotes the normalized IBD (nIBD) value of the recipient-donor pair, n I B D R , B denotes the nIBD value of the recipient-background pair. n I B D = C o u n t _ I B D T o t a l _ I B D , Count_IBD = shared IBD counts between group1 and group2, T o t a l _ I B D = N 1 * N 2 . N₁ and N₂ are the sample size of group1 and group2, respectively. That way, rIBD >0 indicates that commercial breeds (recipient) shares more IBD traces with Chinese indigenous (donor) than Western background, and thus denotes introgression from the donor into the recipient population. In contrast, a negative rIBD value indicates that the number of haplotypes shared by the recipient and the background population is greater than that shared with the donor at that locus. We then performed a Z-transformation of the rIBD values with the mean and standard deviation values of overall IBD from all donor-recipient pairs. We independently selected the presumed introgression bins with a Z-rIBD threshold of μ+2σ for every pair, where μ and σ are the mean and standard deviation of Z-rIBD values. Positive significant Z-rIBD values are thus indicative of the presumed introgression from Chinese breeds into the Western commercial pig.

To measure the coincidence of significant positive/negative Z-rIBD fragments between different donor recipients, we calculated the overlapping ratio by the following formula: O v e r l a p p i n g   r a t i o C , A : B = 1 − N o n e o v e r l a p p e d   C o u n t s C , A + N o n e o v e r l a p p e d   C o u n t s C , B T o t a l   C o u n t s C , A + T o t a l   C o u n t s C , B Where A, B, C denote the three populations. When we compare the overlapping level between “C - > A” and “C - > B”, population C denotes one donor population while A and B denote two different recipients. To compare the overlapping level between “A- > C” and “B- > C”, population C denote on recipient while A and B denote two different donors. N o n e o v e r l a p p e d   C o u n t s C , A : the number of positive/negative fragments/bins shared between C and A but not shared with B. T o t a l   C o u n t s C , A : the total number of positive/negative fragments shared by C and A.

We performed a genome scan to detect recent adaptive introgression events using polymorphism data from the recipient populations only, using the VolcanoFinder-v.1.0 tool (Setter et al., 2020). The ancestral genome was constructed with 16-way Enredo-Pecan-Ortheu multiple alignments files, downloaded from the Ensembl v.103 databases (https://www.ensembl.org/). We obtained the allele frequency and the unnormalized site frequency spectrum files required for Volcanofinder, and performed the analysis according to the standard workflow from VolcanoFinder (https://doi.org/10.5061/dryad.7h44j0zr7) (Setter et al., 2020). Finally, variants were polarized by the ancestor alleles status and used as the input for VolcanoFinder-v1.0 (-big 30000, -1 1 1) (Setter et al., 2020).

To locate important introgressed segments, we merged the consecutive significant ZrIBD bins into one introgression segment. We ranked introgression segments by segment length as the first criterion and average rIBD value as the second. Then, we selected the segments with a length larger than 11 Kb and Log (10) likelihood ratio of selective sweep footprint >11 (the 0.95 quantile). After that, we computed the length of the introgressed segment, selective sweep footprints, average ZrIBD value and average minor allele frequency for every introgression segment and selected the segments that matched all criteria as top candidates for further analysis.

To trace the origin of haplotypes, alleles were first joined into a “FASTA” format sequence from the phased VCF file by an in-house python script. For a genomic region of interest, variants belonging to the same haplotype were joined to a sequence. These haplotypes thus consist of a string of variants derived from the phased VCF. Subsequently, the SNP distance matrix between haplotypes was calculated with SNP-dists v0.7.0 (https://github.com/tseemann/snp-dists). Finally, hierarchical clustering was conducted in R using the gplots package. Paterson’s D-statistics (Patterson et al., 2012) were computed by Dtrios (-j100) of the Dsuite v0.4 (Malinsky et al., 2021) tool package.

We refer to an allele as a “Chinese-derived allele” when it occurs in Duroc and Chinese pigs, but is nearly absent in European wild boars. So the “Chinese-derived allele” should match the following criteria: 1) allele frequency in Duroc pigs≥0.1.2). Allele frequency in any of the Chinese local pig groups≥0.1.3) allele frequency in European wild boars≤0.0125 (i.e., only one European wild boar among 40 boars has that allele and is heterozygous).

Chinese-derived variants were annotated by snpEff-v5.0 (Cingolani et al., 2012). To pinpoint potential causal variants with a high effect on the phenotype, the variants were then ranked using pCADD’s PHRED score. Briefly, the pCADD is the “pig combined annotation dependent depletion”, a model to score single nucleotide variants in pig genomes in terms of their putative deleteriousness, or effect on phenotypes, based on a combination of annotations, see (Gross et al., 2020). The pCADD model is a pig-specific variant of the original CADD model that was developed for human aimd aims to discriminate neutral variants from variants with high impact. Then, “Candidate variants” were selected by the following principles: 1) PHRED score >4.3 (the whole genome mean value). 2) missense variant, 3′UTR variant, or 5′UTR variant.

The LD level of the proxy SNPs with other variants from the sequence data was calculated by plink v1.90b6 (--ld-snp new14_97387849 --ld-window 3000 --ld-window-kb 3000 --r2 --ld-window-r2 0) in the Duroc population. The mean r2 values between proxy SNP and other variants in every block were used as the LD level of the proxy SNP and that block.

To be able to test phenotypic effects of the Chinese-derived introgressed haplotypes on chromosome 14, we wanted to expand our sample size by incorporating genotype data from commercial Duroc pigs. The genotype data was obtained from routinely screened pigs from Topigs Norsvin pigs that were genotyped by the (Illumina) Geneseek custom 50 K SNP chip with 50,689 SNPs (50 K) (Lincoln, NE, USA). The chromosomal positions are based on the Sscrofa11.1 reference assembly. In our set of re-sequenced Duroc pigs, SNPs were filtered using the following requirements: Each marker had a MAF greater than 0.01, a call rate greater than 0.85, and an animal call rate >0.7. SNPs with a p-value below 1 × 10⁻⁵ for the Hardy-Weinberg equilibrium exact test were also discarded. All pre-processing steps were performed using plink v1.90b3.30 (Chang et al., 2015). SNPs on chromosome 14 were retained for further LD analyses to identify the SNP in highest LD with the candidate variants on the introgressed haplotypes. We tested LD between the candidate SNPs from the sequence in the Asian derived haplotype and SNPs on the 50 K chip by usingPlink-1.9 (--ld-snp new14_97387849 --ld-window 3000 --ld-window-kb 3000 --r2 --ld-window-r2 0).

To estimate the impact of genotypes on production traits, we used the genotype data for the candidate SNP from 11,255 Duroc animals (not all animals have all phenotypes) to test the association of our introgressed allele with the following traits: daily gain from birth to Tstart (25 kg) for 9,921 animals, daily gain from Tstart to the Tend (25–120 kg) in 10,986 animals, backfat at 120 kg (Tend) in 7,192 animals, lean meat percentage and loin depth at the end (120 kg) in 7,688 and 7,192 animals respectively. The corrected phenotypes for all traits of each animal were obtained from the routine genetic evaluation by Topigs Norsvin. Then, for each trait, we conducted a Welch’s t-test (significance threshold p < 0.05) to test for differences in phenotypes of the different genotypes at our candidate SNP, that were assigned to either European or Asian background.

We collected 730 samples (Supplementary Table S1) from NCBI (https://www.ncbi.nlm.nih.gov/) and conducted SNP Calling with freeBayes-v1.1 (Garrison and Marth, 2012). After strict quality control, 19, 656, 271 SNPs and 592 samples were retained for further analyses, including 193 Chinese indigenous pigs, 30 Asian wild boars, 13 Western local pigs (Yucatan mini-pigs), 40 Western wild boars, 298 Western commercial pigs, and 18 samples from suidae in Southeast Asian islands (Supplementary Table S1).

We predefined the groups of Chinese pigs according to our previous analysis (Peng et al., 2022) and their geographical origins (Wang et al., 2011) (Supplementary Figure S1). The NJ-tree, t-SNE dimensionality and global ancestry analysis were used to dissect the genetic structure of our samples. The NJ-tree and ancestry inference separate Western and Chinese-derived pigs (Figures 1A–C). Chinese animals clustered into a monophyletic clade, and different Chinese origins clustered into sub-clades except for Chinese Northern pigs (Figure 1A). For Western pigs, every commercial line clustered into a monophyletic clade (Figure 1A) and breeds were clearly distinguished in the t-SNE plot (Figure 1B). Admixture analysis was consistent with this pattern, with increasing values of K above six indicating local ancestry for European pigs, and Asian substructure was best captured with K = 14, when the cross-validation reached a plateau (Figure 1C and Supplementary Figure S2). After we removed a few of the eastern Chinese samples that showed hybrid ancestral components in the admixture result, we assigned the China and Western pig breeds to specific clusters according to geographical sources and genetic relationships.

Finally, Duroc (DUC), Dutch Large White (LWHNL), French Large White (LWHFR), Dutch Landrace (LDRNL), and American Landrace (LDRUS) were recognized as five distinct Western commercial lines (WS). European wild boars (EUW) plus local European pigs (EUD) were defined as the Western background population (WB). Moreover, Southern (SCN), Eastern (ECN), Northern (NCN), and Southwestern (SWCN) Chinese pigs were defined as the four Chinese local groups (AB). We combined Chinese wild boars (CNW), Korean wild boars (KRW), and Thai wild boars (THW) into the Chinese background population (AB) (Supplementary Table S1).

We assessed local signatures of introgression in the Western pig genomes using an IBD haplotype sharing method. There are large introgressed fragments and introgression clusters from China to Western pigs (Figures 2A–D). On a genome-wide scale, the proportion and local regions of putative introgression are highly diverse between different donor-recipient pairs. The highest proportion of introgression into Western commercial breeds is NCN, followed by SCN (Figures 2A–D). The amount of introgression varied between the European breeds, with most putative introgression segments from Chinese pigs found in Large White breeds and the French Large White line in particular (Total length of Chinese-derived segments is 33.82 Mb for DUC, 25.57 Mb for LDRNL, 12.68 Mb for LDRUS, 47.29 Mb for LWHFR, 33.49 Mb for LWHNL. Figures 2A–D and Figure 3F). The putative introgression fragments also varied in length and number (Figures 2E, 3E). The longest introgression fragments reach ∼1.2 Mb between LWHFR and NCN (Total length: ∼38 Mb), but only ∼0.34 Mb between LDRNL and ECN (Total length: ∼1 Mb). For any Western commercial line, the average introgressed segment length from NCN is longer than from other Chinese populations (Figure 2E), suggesting a relatively recent genetic exchange between NCN and Western pigs.

We studied the number of genes affected by the introgression fragments for different donor-recipient pairs. Results are consistent with the total introgression length (Figure 3F). Most of the genes affected by introgression from local Chinese pigs into Western commercial pigs are specific for every donor-recipient pair (Figures 3A–E). Furthermore, most introgressed genes are from NCN to Large White (LWH), especially LWHFR (428 genes, Figure 3C).

We further compared the degree of overlap of positive/negative Z-rIBD segments among donor-recipient combinations to study the global introgression differences. We found that the overall positive Z-rIBD (introgression footprint, see method) patterns are less similar than negative Z-rIBD patterns in different donor-recipient pairs (Figure 4, And Supplementary Figure S3). For donor groups from different sources in China, the degree of overlap significant positive Z-rIBD segments between donor-recipient pairs range from 3.3% to 20.56% (Supplementary Table S2), but for negative Z-rIBD segments, it is 69.73%–93.51% (Supplementary Table S2). The degree of overlap of positive Z-rIBD segments is much lower than the negative Z-rIBD segments (Supplementary Table S2), suggesting specific introgression.

Here we take Duroc as an example. The difference in the positive peaks of Z-rIBD is noticeable (Figure 4). There are peaks at different locations or heights on chromosome nine for ECN, NCN, and SCN (Figures 4A–C), but there is no significant Z-rIBD peak for the SWCN-DUC pair (Figure 4D). On chromosome 11, there are peaks located at 34–39 Mb with different heights or widths (Figures 4A–D). Likewise, on chromosome 15, the highest peak is located at different positions for the Chinese groups (Figures 4A–D) except for NCN and SWCN (Figures 4B–D). Suggesting that pigs from different regions in China contributed differently to Western commercial pig breeds.

There is a large difference in the introgression patterns between specific Chinese groups and Western commercial lines (Figure 5, And Supplementary Figure S4). The degree of overlap of positive Z-rIBD segments for the specific Chinese local pigs to different European commercial pigs ranges from 0% to 34.83% (Supplementary Table S3). In contrast, for the negative Z-rIBD segments, it ranges from 27.48% to 49.81% (Supplementary Table S3). However, the differences in the introgression patterns within related breeds from a specific Chinese group are smaller. Dutch and French Large White breeds show a similar introgression pattern from North Chinese pigs compared with other commercial pigs (Figures 5D,E, and Supplementary Table S3). The degree of overlap of positive Z-rIBD for these breeds is as high as 34.83% (Supplementary Table S3). In contrast, this is only around 10% compared to the other Western commercial lines (Supplementary Table S3). A broad peak on chromosome 3 (chr3:48–52 Mb) was found in Dutch Large White and French Large White (Figures 5D,E), but not in the other breeds (Figures 5A–C).

In the Landrace breed, the degree of overlap of positive Z-rIBD is as high as 26.55% (Figures 5B,C, and Supplementary Table S3). A significant introgression signal on chromosome 17 (CHR17:17–18 Mb) is observed in both Dutch and American Landrace (Figures 5B,C) but not seen in the other breeds (Figures 5A,D,E). Besides, the Z-rIBD peaks mentioned above are different in the different pig lines. The observed difference in introgression signal from specific Chinese groups to related Western commercial lines reflects a difference in the extent of introgression. This suggests that gene flow occurred mainly in the early stages of commercial pig breeding rather than after the differentiation of the lines. However, the tendency of artificial selection caused changes in signal strength.

We observed a cluster of Duroc-specific introgression signatures spanning ∼2.65 Mb on chromosome 14 (chr14: 95.68–98.33 Mb) (Figure 6). Such a strong introgression and selection signal is not seen for the other commercial pigs at that region (Supplementary Figures S5, S6). This introgressed region appears to be a set of segments derived from Chinese pigs in the Duroc population. The Z-rIBD value for SCN-DUC is up to 5.65 for segment 3 (the mean Z-rIBD value is 2.48 for segment 1, 3.90 for segment 2, 2.19 for segment 3 and 2.77 for segment 4, Figure 6A). Except for SWCN-DUC in segment 1 (mean Z-rIBD = 5.89, Figure 6D) and segment 3 (mean Z-rIBD = 4.74, Figure 6D), the mean Z-rIBD values is highest in the SCN-DUC pair (Figures 6A–D). We also observed a lower minor allele frequency than expected by chance (0.04 for this region but 0.12 for whole-genome) in Duroc (Figure 6E). These signatures are located within a strong adaptive selection region (chr14: 92–101 Mb) on the Duroc genome (Figure 6F). Moreover, there are five candidate genes in this region: PCDH15, MBL2, DKK1, PRKG1, and CSTF2T (Figure 6G).

Additionally, PCA plots of Duroc and Chinese pigs from SNPs across the full genome and local SNPs in this region display a strong discondancy (Supplementary Figure S7). The clustering of the Duroc and Chinese pigs in this region hint at introgression, evident from the big difference between the global and local PCA analyses. Combining the above results, we suspect that this haplotype in the Duroc genome was inherited from SCN or SWCN pigs.

To trace the sources of the haplotype, we then calculated a distance matrix between individuals for every segment by SNP-dists v0.7.0 followed by hierarchical clustering in R-4.0 using the gplots package (Warnes, 2020). The results (Figure 7) show that most Duroc pigs clustered together with Chinese pigs (especially with ECN, SCN, and SWCN), in sharp contrast to LWH and Landrace (LDR). The LWH and LDR clustered with Western background populations on segment 1 and segment 4 (Figures 7A,D). In the clustering of segments 2 and 4, more SCN pigs are located witnin Duroc clades (Figures 7B,D), suggesting that fragment 2 is more likely derived from SCN pigs. The results of the ABBA-BABA test (D-statistics) highlights that Duroc shares more derived alleles with SCN than other Chinese pigs for segment 2, segment 3, and segment 4 (Table 1). Moreover, a high degree of linkage disequilibrium (LD) in this region (r2 = 0.56, Supplementary Figure S8) was found. According to the dist trees, D-statistics and the degree of LD, we believe that the haplotype (chr14: 95.68–98.33 Mb) in the Duroc genome is derived from SCN pigs.

We further investigated SNP allele frequencies in the four segments in different pig populations (Table 2, And Supplementary Table S4). There are many alleles with low (≤0.0125) frequencies in Western wild boars but high frequencies in Duroc pigs for each of the segments (123 variants in segment 1,118 variants in segment 2,436 variants in segment 3,383 variants in segment 4, Supplementary Table S4). Furthermore, the derived alleles in Duroc pigs at these loci seem to have undergone strong selection (Figure 6F). We believe these are candidate alleles derived from Chinese pigs due to their moderate allele frequencies in Chinese pigs (Table 2, And Supplementary Table S4). Seven candidate mutations (Supplementary Table S5) were selected from the putative Chinese-derived set of alleles that potentially have a high functional impact (see methods). These variants are likely to have a strong impact on the phenotype as derived from the pCADD model, with the strongest located within the three prime UTR region of PRKG1.

We analyzed genotype and phenotype data of 11,255 animals from a commercial Duroc population to assess the potential phenotypic impact of the introgressed haplotypes. We screened the (Illumina) Geneseek custom 50 K SNP array for SNPs in highest LD with the introgressed haplotypes, and a SNP (INRA0045978) was selected as a proxy for the introgressed segment due to its high LD (r2 range from 0.65 to 0.73) with the seven candidate alleles in the Duroc population (Supplementary Table S6, Supplementary Figure S9). Next, we used the genotypes for this selected SNP from 11,255 Duroc animals from the same commercial breed to test the association of INRA0045978 with a set of production traits (See methods).

We found a significant association with backfat (genotype “0/0” versus “1/1; t-test p-value 0.016; Figure 8C) with a and with loin depth (genotype “0/1” versus “1/1; t-test p-value 0.028; Figure 8E and Supplementary Table S7). The INRA0045978 SNP has a low Duroc reference allele frequency in Western wild boar (0.0125) but higher in Chinese pigs (0.5517) and Duroc (0.8782). These results suggest that the PRKG1-haplotype may decrease backfat (mean difference of 2.3 mm) and increase loin depth (mean difference of 6.1 mm) in Duroc pigs.