Whole‐genome SNP markers reveal conservation status, signatures of selection, and introgression in Chinese Laiwu pigs

Abstract Laiwu pigs are a Chinese indigenous breed that is renowned for its exceptionally high intramuscular fat content (average greater than 6%), providing an excellent genetic resource for the genetic improvement of meat quality of modern commercial pigs. To uncover genetic diversity, population structure, signature of selection, and potential exotic introgression in this breed, we sampled 238 Laiwu pigs from a state‐supported conservation population and genotyped these individuals using GeneSeek 80K SNP BeadChip. We then conducted in‐depth population genetics analyses for the Laiwu pig in a context of 1,116 pigs from 42 Eurasian diverse breeds. First, we show that the current Laiwu population has more abundant genetic diversity than the population of 18 years ago likely due to gene flow from European commercial breeds. Both neighbor‐joining (NJ) and principal component analyses indicate the introgression of European haplotypes into Laiwu pigs. The admixture analysis reveals that an average 26.66% of Laiwu genetic components are of European origin. Then, we assigned the tested individuals to different families according to their clustering patterns in the NJ tree and proposed a family‐based conservation strategy to reduce the risk of inbreeding depression in Laiwu pigs. Next, we explored three statistics (ROH and iHS and EigenGWAS) to identify a list of candidate genes for fat deposition, reproduction, and growth in Laiwu pigs. Last, we detected a strong signature of introgression from European pigs into Laiwu pigs at the GPC6 locus that regulates the growth of developing long bones. Further association analyses indicate that the introgressed GPC6 haplotype likely contributed to the improvement of growth performance in Laiwu pigs. Altogether, this study not only benefits the better conservation of the Laiwu pig, but also advances our knowledge of the poorly understood effect of human‐mediated introgression on phenotypic traits in Chinese indigenous pigs.


| INTRODUC TI ON
Laiwu pigs ( Figure 1) were originally distributed in Laiwu City, Shandong Province of China. This breed is renowned for its desirable meat quality and exceptionally high intramuscular fat content in pork with an average value of ~6% as compared to less than 2% in European breeds (Chen, Fang, Wang, Wang, & Zeng, 2017). Hence, Laiwu pigs provide excellent genetic materials for the improvement of meat quality in the present-day pig industry. For this reason, Laiwu pigs have been included in the conservation list of China's livestock and poultry genetic resource by the Ministry of Agriculture of China (Wang et al., 2011).
Archaeological evidence shows that modern and ancient (5,300-6,500 years ago) Laiwu pigs had nearly identical head bones, indicating that Laiwu pigs have a long breeding history of more than 5,000 years. In the early 20th century, Yantai and Qingdao cities of Shandong province were colonized by Britain and Germany, respectively. Western pig breeds such as Berkshire, Yorkshire, and Duroc were then introduced to the Shangdong province to cross with local breeds (Wang et al., 2011). It hence raises the possibility that Shangdong indigenous pigs including the Laiwu pig could have genetic components introgressed from Western breeds. The possibility was supported by our recent study based on whole-genome SNP markers. We showed that approximately one-fourth of genetic components in the genome of Laiwu pig are of European origin (Huang et al., 2019).
However, the genomic distribution and phenotypic effect of the introgressed European haplotypes in Laiwu pigs remain elusive.
In addition, the conservation population of Laiwu pigs is restricted to a single state-supported farm (Wang et al., 2011). The conservation status including genetic diversity, potential risk of inbreeding depression, and family structure is largely unknown in Laiwu pigs.
In this study, we explored whole-genome SNP markers to investigate the conservation status, genomic signatures of selection and introgression and their effect on phenotypic traits of all Laiwu pigs (n = 238) from the current conservation population in a context of 1,116 individuals from 42 Eurasian diverse breeds.
The findings enable us to propose a reliable and sustainable conservation strategy for the Laiwu pig, and provide novel insight into historical contribution of European lineages to Chinese indigenous pigs.

| Ethics statement
All experimental procedures were approved by the Animal Care and Use Committee of South China Agricultural University, Guangzhou, China.

| Sample and genotyping
Ear tissues of 238 Laiwu pigs were sampled from the national conservation farm for Laiwu pigs in Shandong province of China in 2018. These individuals represent the most comprehensive genetic diversity of the Laiwu pig. Genomic DNA was extracted from the ear tissue using a traditional phenol/chloroform method. DNA samples were genotyped for 68,516 SNPs using GeneSeek Genomic Profiler Porcine HD BeadChip (Neogen Corporation, USA) according to the manufacturer's instructions.
To investigate the genetic diversity and population structure of Laiwu pigs in a global perspective, we integrated the 80K SNP data with the Illumina 60K (61,565 SNPs) SNP data of 878 pigs from 35 Chinese pig breeds, one hybrid breed and six European commercial pig breeds (Table S1). The 878 pigs included 18 Laiwu pigs sampled in 1999. Of the 878 pigs, 299 individuals were tested in a previous study (Ai, Huang, & Ren, 2013), 107 individuals were reported in (Wang et al., 2018), and the other 472 individuals were studied in (Xu, Wang, et al., 2019;Xu, Sun, et al., 2019;Yang et al., 2017). We converted the raw data to PLINK v1.9 (Purcell et al., 2007) input files and obtained a common set of 42,464 SNPs from 1,116 individuals. Then, we conducted the following quality control procedures for the merged SNP data using PLINK v1.9 (Purcell et al., 2007): (1) We randomly selected one individual from one pair of highly related animals with an identity-by-state score of greater than 0.99 in Laiwu pigs, (2) retained SNPs with minor allele frequencies (MAF) of no less than 0.01, (3) removed SNPs and individuals with call rates of lower than 90%, and (4) discarded all unmapped SNPs and those on sex chromosomes. A final set of 35,027 SNPs from 1,111 pigs were used for the following analyses.

| Genetic diversity analyses
Expected heterozygosity (He), observed heterozygosity (Ho), effective population size (Ne), and linkage disequilibrium (LD) decay were calculated to evaluate the genetic diversity of the tested populations, including the current (LWH, n = 233) and early-day conservation population (LWH1, n = 18) of Laiwu pigs. He and Ho were calculated using PLINK v1.9 (Purcell et al., 2007) with default setting. Pair-wise LD was evaluated by the correlation coefficient (r 2 ) between alleles at two separate SNP loci using PLINK v1.9 (Purcell et al., 2007) under the default setting. According to a previous study (Sved, 1971), Ne was estimated based on the equation Ne t = (1/4c) (1/r 2 − 1), where Ne t is the effective population size of t generations age and is calculated as t = 1/2c, r 2 is the LD between pairwise SNPs, and c is the genetic distance in Morgan between a pair of SNPs (Tortereau et al., 2012).

| Inbreeding coefficient
Two measures of genomic inbreeding coefficient were calculated for each population using PLINK v1.9 (Purcell et al., 2007).
(1) SNPbased inbreeding coefficient (F) that was calculated using 25,839 SNPs with pair-wise LD values of less than 0.5, and the command was set as "--indep-pairwise 50 10 0.5." (2) Runs of homozygosity (ROH)-based inbreeding coefficient (F ROH ) that was measured by the ratio of the total length of ROH to the length of autosomes (2.45 Gb in this study; Mcquillan et al., 2008). ROHs were identified for each individual using PLINK v 1.9 (Purcell et al., 2007) with the following parameters (Shi et al., 2020;Xie et al., 2019): (a) The minimum number of SNPs in a sliding window was 50; (b) the minimum ROH length was set to 1 Mb to eliminate the impact of strong linkage disequilibrium; (c) each ROH need contain a minimum of 80 consecutive SNPs, which was calculated by following equation (Lencz et al., 2007),
The two matrices were then explored to construct neighbor-joining (NJ) trees of individuals and populations using PHYLIP v3.69 (Felsenstein, 1989 (Yang, Lee, Goddard, & Visscher, 2011), and the first two PCs were plotted via in-house R scripts. The ADMIXTURE v1.30 software (Alexander, Novembre, & Lange, 2009) was employed to infer the proportion of introgressed ancestry in the tested populations. To reduce the effect of ascertainment bias, we used the genotype data of the pruned 25,839 SNPs with LD (r 2 ) values of less than 0.5. To avoid sampling bias, we randomly selected 10 pigs from each population to calculate the ancestral lineage compositions. Considering the large sample size and heterogenous genetic background of current Laiwu pigs, we classified the 233 Laiwu pigs to 23 groups (each group contained ten pigs, except the last group with 13 pigs) with other random sampling breeds to conduct the Admixture analyses and obtained well-consistent results (data not shown). The optimal number was determined by cross-validation error. Finally, the inferred population structure for one of group was visualized using in-house R scripts.

| Detection of introgression
The f3 test ( Malinsky (Malinsky et al., 2015) was calculated, using a 10 SNPs sliding window with 2 SNPs stepping . Positive values of the f dM stand for introgression between P3 and P2 and negative values for introgression between P3 and P1. P values then estimated by Z-transformed f dM values assuming the standard normal distribution, and windows with p < .01 (f dM = 0.611) were identified as significantly introgressed genomic regions. Moreover, absolute divergence (dxy) and nucleotide diversity (π) in genome-wide were calculated with the python script from https://github.com/giber t-Fab/ABBA-BABA (Martin, Davey, & Jiggins, 2015), using a 10 SNPs sliding window with 2 SNPs stepping. We also used fastPHASE (Scheet & Stephens, 2006) to construct haplotypes using 21 SNPs at the significant selection and introgression gene GPC6 (11:62,436,560,908) in tested populations. Finally, we downloaded the resequencing data of 259 Eurasian pigs (including six Laiwu pigs) from previous studies  and retrieved the GPC6 region (13,899 qualified SNPs) to perform NJ tree analysis.

| Association analyses
To examine the putative effect of the introgression region, we genotyped 365 pigs from the Sujiang breed that was derived from a cross between Chinese Jiangquhai pigs and European Duroc pigs using the GeneSeek Genomic Profiler Porcine HD BeadChip (Xu et al., 2020). A total of 21 qualified SNPs extracted from GPC6 region on pig chromosome 11. The linear mixed model was fit to this data using GEMMA software (Zhou & Stephens, 2012) to test for the significant effect of SNPs on body weight, body length, and chest circumference. The statistical model was described as following: where y is the phenotype; W is a matrix of covariates (i.e., fixed effects that contain batch, age and a column of 1s); α is a vector of corresponding coefficients that includes the intercept; x is a vector of SNP genotypes; β is the effect size of SNPs; u is a vector of random polygenic effects with a covariance structure that follows a normal distribution u ∼ N(0, K × Vg), where K is a genomic relationship matrix derived from independent SNPs and Vg is the polygenic additive variance, and ε is a vector of random errors.

| Genetic diversity of Laiwu pigs
The current conservation population of Laiwu pigs had the greatest value of He (0.284) and the third greatest value of Ho (0.281) among 35 Chinese indigenous pig breeds tested in this study. In contrast, the early-day Laiwu population (LWH1) had the third least values of Ho (0.15) and He (0.15) among these 35 breeds (Table S1). The SNPbased (F) and ROH-based (F ROH ) inbreeding coefficients showed positive correlation (r = 0.513, p < .001). We subdivided F ROH values into three categories, that is, those derived from ROHs of less than 5 Mb (F ROH 0-5Mb ), from 5 to 10 Mb (F ROH 5-10Mb ), and greater than 10 Mb (F ROH >10Mb ;  (Table S2).
LD extent in each population was estimated as the physical genomic distance at which the genotypic association (r 2 ) was 0.3 (Ai et al., 2013). LWH1 pigs had the longer of LD extent (r 2 0.3 = 195.11 kb) among Chinese pigs. The LD value was even greater than those of European pigs from Large White, Landrace, Duroc, and Hampshire breeds, and was more than three and two times than that of LWH pigs (r 2 0.3 = 62.14 kb) and LWH2 population (r 2 0.3 = 79.31 kb). The average effective population size (Ne) of each population was estimated using the method as previously described (Herrero-Medrano et al., 2013). To reduce potential sampling bias, we calculated Ne for the past five generations. Ne ranged from 84 in Ganxi pigs to 241 in LWH pigs. The Ne of LWH was twice as much as that of LWH1  Table S1 (Ne = 120) and was roughly equal to the actual number of the LWH population (No = 238). Surprisingly, the Ne of LWH2 (Ne = 271) greater than whole LWH pigs, we speculated some individuals with large LD were removed in LWH2 population (Table S1).

| Principal component analysis
European modern breeds, European-Chinese hybrid breed (Sutai), Laiwu pigs, and other Chinese breeds formed four groupings in the PCA plots ( Figure S1). PC1 accounted for 19.9% variance between Chinese and European pigs, and PC2 explained 4.8% variance between Laiwu and other Eurasia pigs. The two Laiwu populations (LWH1 and LWH) clustered together, and a number of LWH individuals showed close relationships with European modern pigs ( Figure   S1), which is consistent with the NJ clustering pattern.

| Population structure and subfamilies in Laiwu pigs
We made a close examination on the ancestral lineage compositions of all 233 tested Laiwu pigs using the ADMIXTURE analysis (Alexander et al., 2009). When K = 2, a proportion of Chinese lineage ranging from 29.10% to 86.64% was observed in these 233 Laiwu pigs with an average value of 73.25% and a standard deviation (SD) of 12.22% (Table S3) Figure S3). We highlighted nine and two individuals as the outliers, because these individuals had the ROH and F values of greater than the mean value plus two standard deviations. Hence, special attention should be paid to these highly inbred individuals to avoid inbreeding depression. As expected, the LWH3 and LWH4 individuals had lower F ROH (0.076) and F (0.045) values than LWH2 pigs (F ROH = 0.164; F = 0.270), which is attributable to a high proportion of European lineage in the LWH3 and LWH4 pigs (Table S3).

| Runs of homozygosity
We found that 0.39% of SNPs not occur in any ROHs of Laiwu pigs, and 2 chromosomes had long non-ROH cold-spot regions: SS2  (Table S4).

| iHS
The greatest value |iHS| values on SSC1, and the top 1% SNPs (452 SNPs; Table S4) were considered candidate loci under selection (Figure 5b). We observed a significantly positive correlation coefficient between the ROH and iHS statistics (r = 0.071, p < 10 -16 ). This supports our assumption that the significant ROH regions are not only due to demography but also positive selection (Xu, Sun, et al., 2019;Zhang, Guldbrandtsen, Bosse, Lund, & Sahana, 2015).

| EigenGWAS
In the EigenGWAS analysis, we identified 110 SNPs surpassing the genome-wide significant threshold (0.05/35027; Figure 6a, Table S4). The F ST and EigenGWAS statistics showed significant positive correlation (r = 0.91, p < 2.2 × 10 -16 ). Of note, the top SNP located in a large region (47.87-64.05 Mb) on SSC11.  (Table 1); fatness-related candidate genes including ESR1, TP53INP1, and GPC6; and reproduction-related candidate genes including VIP, ESR1, AKAP12, and TGFBR3 (Table 1). A total of 333 quantitative trait loci (QTL) have been reported around these regions according to the pig QTL database (https://www.anima lgeno me.org/ cgi-bin/QTLdb/ SS/index), and 66.97% are related to meat and carcass traits, 13.81% are associated with health traits, 13.21% are related to production traits, and 4.2% are associated with reproduction (Table S6). It is worth mentioning that only two significant SNPs on SSC4 were detected by the three methods ( Figure S4a). One SNP locates at 5.8 kb downstream of the ZNF326 gene, and the other SNP is in an intron of the GDF6 gene.

| Detection of introgression
There were 40 extreme Z-scores of the f3 test indicated that LWH pigs were significantly admixed by LWH1, eleven Chinese pig breeds (CNref), six European commercial pig breeds (EUref; Table S7). The pigs, 3 BMEI pigs, 3 LW pigs, 3 USBK pigs, and 1 LR pig (Table S9). Most of the individuals in LWH1 (22/36) and LWH2 (218/300) carried this main haplotype, and this haplotype was not found in any other Chinese pig breeds, excepted for BMEI pigs which had been introgressed by European pigs reported in previous study (Huang et al., 2019). Then, we retrieved genomic DNA sequence of GPC6 (13,899 qualified SNPs) from whole-genome sequence data of 259 Eurasian pigs  Figure 7, Table S10). The IBS-based NJ tree showed that three Laiwu pigs clustered with Large White pigs and European wild boars, supporting the introgression event at the GPC6 loci.

| Associations of introgressed haplotypes with body size of Laiwu pigs
We identified three significant (p < .05) SNPs of GPC6 by linear mixed model association analysis (Table 2, Table S11). We focused on the locus at 62755707 bp (SNP ASGA0051220) with significant selection  Figure 8a). The association analysis showed that ASGA0051220 significantly associated with body weight (p = .005) and chest circumference (p = .024), and individuals with the TT genotype had larger body size than those with genotypes TC and CC (Figure 8b-d). Our findings suggest that the introgressed European haplotypes have been preferentially selected for possibly improving growth performance, leading to their distribution at high frequencies in Laiwu pigs.

| Genomic signature of introgression in Laiwu pigs
In this study, multiple evidences support that European haplotypes had been introgressed into Laiwu pigs. First, Laiwu pigs did not cluster with the Chinese major clade in the NJ tree and multiple Laiwu individuals displayed close genetic relationships with European F I G U R E 7 Neighbor-joining (NJ) tree based on an identity-by-state matrix among 259 Eurasian pigs. NJ tree based on 13,899 qualified SNPs in the region of GPC6. CN, Chinese local pig breeds; AWB, Asian wild boar; EU, European modern pig breeds; and EWB, European wild boar modern pigs in the PCA plot. Moreover, the analyses of Admixture, f3, and D-statistic clearly showed gene flow from European modern pigs to Laiwu pigs, and 13 Laiwu pigs had more than half genetic components of European origin. These findings were in agreement with the historical documents that European breeds were introduced to cross with local pigs in Shangdong province of China at three periods during the past century, including from 1920s to 1940s, from 1950s to 1970s, and from 1980s to 1990s (Wang et al., 2011). It is most likely that indiscriminate crossbreeding between Laiwu pigs and introduced European modern pigs occurred during these periods, leading to an average of 26.66% (Admixture, qpAdm, and f dM ) genetic components from European modern pigs in the present conservation population of Laiwu pigs.

| Effect of introgression events on the genetic diversity of Laiwu pigs
We analyzed the genetic diversity of two Laiwu populations, the present-day and early-day, in a world panel of pigs. Six statistics (Ho, He, Ne, F, F ROH, and LD) collectively show that the current population has more abundant genetic diversity than the early-day population. There findings were not expected, because the conservation population of Laiwu pigs was raised in a single farm without introduction of Laiwu pigs from other locations for the past decades. It is thus inevitable that inbreeding occurred in this population to a certain degree, resulting in a decreased level of genetic diversity.
One reasonable explanation for this inconsistence is that recent admixture with exotic lineage increased the genetic diversity level in the current Laiwu population. This assumption is supported by the fact that Laiwu pigs have genomic signatures of admixture with European modern pigs, and the fact that the nucleus population of Lulai pig, a synthetic line derived from a cross between Laiwu and Large White, was also raised in this conservation farm for Laiwu pigs.
It is likely that the two populations were not well managed, and occasional admixture between individuals from the two populations led to gene flow from Lulai pigs to Laiwu pigs and consequently increased genomic variability of the current conservation population of Laiwu pigs. This speculation was also confirmed by the result of nucleotide diversity (π) promotes 25.11% in the introgressed regions of current Laiwu pigs.

| Sustainable conservation strategy for Laiwu pigs
Of the 233 Laiwu pigs tested in this study, 150 individuals (LWH2) had more than the average proportion (73.25%) of Chinese lineage, and the average proportion of LWH2 was 81%. We suggest that these 150 individuals should be included as the nucleus conservation population of Laiwu pigs, and the considerable level of genetic diversity has been retained from the current population. As revealed by the NJ tree, the 150 individuals pertain to ten families. The previous studies demonstrated that family rotational mating (F: R) could maintain 90% of genetic diversity in a livestock for more than 100 years (Lu, 2013) and effectively reduce inbreeding in populations (Honda, Nomura, & Mukai, 2004;Windig & Kaal, 2008

| Genomic signatures of selection in Laiwu pigs
We performed ROH, iHS, and EigenGWAS to uncover genome-wide footprints caused by natural and artificial selection in the 150 LWH2 pigs. We found few overlaps between the candidate regions identified by the three approaches. This is not surprising as there are differences in the statistics underlying each approach allowing to reveal the signatures of different types of selection across different timescales (Mastrangelo et al., 2020). Runs of homozygosity mainly results from population phenomena such as genetic drift, population bottleneck and inbreeding, the intensive artificial selection of superior animals as one factor has also reshaped the ROH patterns in various regions of the genome (Peripolli et al., 2017). The same as previous studies reported, some short and shared ROH islands were observed to be overlapped with regions under selection based the iHS testing, and showed significant positive correlation (Xu, Sun, et al., 2019;Zhang et al., 2015). The iHS test is especially powerful in detection of recent selection that has swept the selected allele to moderate frequencies but the selected allele has not yet been fixed (Tang, Thornton, & Stoneking, 2007 (Lacombe et al., 2007). ESR1 variants have been associated with litter sizes in pigs (Muñoz et al., 2007), and ESR1 knockout mice develop obesity after sexual maturation (Ohlsson et al., 2000).
We also detected several candidate genes relating to growth and fat deposition traits. The genes MYCT1, BARHL2, ZNF326, ABCC4, and GPC5 are identified as candidate genes for human body height in previous GWAS studies (Fox et al., 2007;Allen et al., 2010;Winkler et al., 2015;Wood et al., 2014). SYNE1 knockout mice show growth retardation and decreased survival rates (Jianlin et al., 2010). GDF6 plays an important role in formation of a diverse subset of skeletal joints (Clendenning & Mortlock, 2012). Mice lacking SOX21 was with reduced growth, and SOX21 variants may be a cause of nonendocrine short stature in human (Cheung, Okano, & Camper, 2017).
Gain-of-function mutations in the phosphatidylserine synthase 1 (PTDSS1) gene cause Lenz-Majewski hyperostotic dwarfism and developed a severe skeletal dysplasia (Piard et al., 2018). Mice lacking out TP53INP1 were prone to obesity (Seillier et al., 2015). It is reported that the body weight, body length, reproductive performance, and lean meat percentage of Laiwu pigs have been improved from 1988 to 2007 (Wang et al., 2011). These improvements likely benefit by the positive selection and made the Laiwu pigs more competitive. These findings advance our understanding of the molecular mechanisms underlying the germplasm characteristics in Laiwu pigs.

| Effect of introgression events on phenotypic traits in Laiwu pigs
Pigs were domesticated independently from European and Chinese wild boars nearly 10,000 years ago (Frantz et al., 2013), and different phenotypic characteristics harbor in Eurasian pigs. Over past decades, European modern breeds have been imported into China and admixed with many Chinese indigenous pig breeds (Wang et al., 2011). Our results confirmed the introgression has occurred from European modern pigs to Laiwu pigs genome. For introgressed regions, it is possible to identify several genes involved the growthrelated GO terms, considering that European modern pigs have been intensively selected for large body size, fast growth and lean pork yield.
We are particularly interested in the most significant selected Loss-of-function mutations in the GPC6 gene cause autosomal-recessive omodysplasia in mice that is characterized by short stature, shortened limbs, and facial dysmorphism. Moreover, GPC6 stimulates Hedgehog (Hh) signaling by binding to Hh and Patched 1 at the cilium and increasing the interaction of the receptor and ligand, and then promotes the growth of developing long bones (Capurro et al., 2017). In addition, the previous GWAS study identified that GPC6 was a candidate gene correlated with the intramuscular fat of Duroc pigs (Ding et al., 2019). We explored genomic sequence data to confirm the introgression event at the GPC6 locus, and further showed that the introgressed European GPC6 haplotype was likely favorable for larger body size and weight and has been preferentially selected in Laiwu pigs, given that three SNPs at the GPC6 locus were significantly associated with large body size in Sujiang pigs. Considering that the relatively small sample sizes of association analysis, and the tested population was not Laiwu pigs, additional studies are needed to definitively determine the role of PGC6. To our best knowledge, this study provides the first example of European introgressed haplotypes at a specific locus affecting growth performance in Chinese indigenous pigs, which advances our knowledge of the poorly understood effect of human medicated introgression on phenotypic traits in Chinese indigenous pigs.

| COMPE TING INTERE S TS
None of the authors have any competing interests in the manuscript.

ACK N OWLED G M ENTS
We are gratefully to Dr. Wanbo Li, Dr. Wen Huang, and three anon-