Data were obtained from two Large White populations with different genetic backgrounds in one Chinese commercial pig company in Shanghai City, which were originated from Canadian and French lines. Feeding and performance testing of animals from these two lines were carried out at two different farms. When the average live weight per batch was approximate 100 kg, individual tests were performed. Body weight and ultrasonic backfat thickness between 11^th to 12^th ribs were measured. The initial and ending dates were recorded. To uniform the data, the measured age was adjusted to 100 kg live weight using the equation (Wang, 2007, 18–19): A G E 100 = m e a s u r e d   a g e + 100   k g − m e a s u r e d   a g e C F , where C F m a l e = m e a s u r e d   w e i g h t m e a s u r e d   a g e × 1.826 and C F f e m a l e = m e a s u r e d   w e i g h t m e a s u r e d   a g e × 1.715 , and the measured backfat thickness was adjusted to 100 kg live weight using the formula (Chen et al., 2019): B F 100 m a l e = m e a s u r e d   b a c k f a t × 12.402 12.402 + 0.106 × ( m e a s u r e d   w e i g h t − 100 ) and B F 100 f e m a l e = m e a s u r e d   b a c k f a t × 13.706 13.706 + 0.119 × ( m e a s u r e d   w e i g h t − 100 ) .

Table 1 shows a summary of the descriptive statistics of AGE100 and BF100 traits. Totally, there were 3,727 observations available for both AGE100 and BF100. Out of them, 2,138 were from the Canadian lines and 1,589 were from French lines. These Large White pigs were born between 2015 and 2020. According to the pedigree information, there were no genetic connectedness between two populations. As seen in Table1, the heritability of AGE100 in Canadian and France lines are 0.15 and 0.30, while the heritability of BF100 in Canadian and France lines are 0.21 and 0.44, respectively.

Genotyping all the individuals with phenotypes was carried out using GeneSeek GGP Porcine HD array. Since the SNP chip is composed of 50,915 probes according to the Sus Scrofa 10.2 version, the autosomal SNPs were further liftovered to the latest version of the pig genome Sus Scrofa 11.1. Thus, 46,258 autosomal SNPs were kept for the subsequent analysis.

Quality control was executed by PLINK (v1.90) (Purcell et al., 2007). Pigs with call rate < 0.9 were excluded. SNPs with a minor allele frequency (MAF) below 0.05 and call rate < 0.9 were excluded in each population. Finally, the remaining autosomal SNPs were 41,172 and 40,506 with the average distances of 53.87 Kb and 54.85 Kb between adjacent SNPs in Canadian and French lines, respectively.

To investigate the population stratification, principal component analysis (PCA) was conducted using the remaining SNPs to obtain eigenvalues and eigenvectors by PLINK (v1.90) (Purcell et al., 2007). In addition, we performed ancestry estimation using ADMIXTURE (v1.3.0) (Alexander et al., 2009). The number of level of genetic structure were estimated from K = 1 to four for all individuals jointly, followed by the cross-validation error (CV) procedure. The linkage disequilibrium (LD, expressed as r ² ) was calculated using PLINK (v1.90) (Purcell et al., 2007) within each line. In this study, the population effect size (Ne) was computed by SNeP software (v1.1) (Barbato et al., 2015).

The mixed model was executed for analyzing the traits under study as following: y = μ + X b + W g + Z u + e (1) where y is the vector of target phenotypes of individuals; μ is overall mean; b is the vector of fixed effects: sex (two levels) and year-season in which seasons were comprised of four levels (Spring:  March to May; Summer:  June to August; Autumn: September to November; Winter:  December to February); g is the vector of the SNP effects, X is the matrix of incidence associating each observations to the pertinent level of fixed effects, W is the incidence matrix relating observations to SNPs effects with elements coded as 0, one and two for genotype A₁A₁, A₁A₂, and A₂A₂, respectively, u is the random additive genetic effect of the individual and is assumed to be distributed as N ( 0 , G σ u 2 ) , G is the genomic relationship matrix and σ u 2 is the polygenic additive genetic variance, Z is incidence matrix for u, e is the random residual and is assumed to be distributed as N ( 0 , I σ e 2 ) , where I is the identity matrix and σ e 2 is the residual variance. Associations between the target traits and the SNPs were analyzed using single-SNP association tests in each population, which were implemented by the mlma option of the software GCTA (Version 1.93.3beta) (Yang et al., 2011).

Furthermore, we integrated two Large White pig lines to identify the candidate genes using combined-population GWAS and cross-population GWAS. In the combined-population GWAS, line effect was taken into account the fixed effect, and the analysis procedure was followed the single-population GWAS mentioned above. The cross-population GWAS also utilized the summary data of single-population GWAS to implement meta-analysis by METAL (version 2011–03–25) (Willer et al., 2010). In meta-analysis, the weighted Z-score model took account of the p-values, direction of SNP effects and the number of individuals. In each case, threshold p-values were set to -log₁₀ (1/SNPs) and -log₁₀ (0.05/SNPs) for suggestive and Bonferroni-adjusted genome-wide significance, respectively. Quantile-quantile (QQ) plot of–log (p-values) was examined to determine how well GCTA accounted for population structure and family relatedness.

In response to intensive artificial selection pressures, the porcine genome has been sculpted signals at the underlying genomic regions harboring functional genetic variants, which are termed as selection signatures (Bertolini et al., 2018). Population differentiation-based methods were performed, including F_ST, hapFLK and runs of homozygosity (ROH). The VCFtools software (Danecek et al., 2011) was used to compute the Weir and Cockerham’s F_ST estimator (Weir and Cockerham, 1984) per site between the two pig lines. The hapFLK statistic (Fariello et al., 2013) was estimated using hapFLK module in python. It should be mentioned that before calculation of hapFLK values, fastPHASE (Scheet and Matthew, 2006) was used to determine the optimum number of haplotype clusters. ROH analysis was performed using PLINK (v1.90) (Purcell et al., 2007), and the parameters were assigned following our previous study (Shi et al., 2020).

In addition, we further investigated the impact of SNPs undergoing selection signatures on the traits under study. First, according to the sizes of selection signals, we sorted the whole SNPs and split them into five sets [0–20%, 20–40%, 40–60%, 60–80% and 80–100%]. To quantify the relative importance of SNPs sets, we calculated one of these SNP sets and the remaining four SNP sets explaining the proportion of phenotypic variance. The statistical model was employed as below: y = μ + X b + u 1 + u 2 + e (2) where u ₁ is the vector of the first random additive genetic effect, which is distributed as u 1 ∼ N ( 0 , G 1 σ u 1 2 ) , where G 1 is the first genomic relationship matrix that was constructed using one of quantile SNPs, u ₂ is the vector of second random additive genetic effect, which is distributed as u 2 ∼ N ( 0 , G 2 σ u 2 2 ) , where G 2 is the second genomic relationship matrix which was calculated using the remaining four SNP sets. Other notations are same as Eq. 1. These two random effect variances can be computed by Qgg R package (Rohde et al., 2020). The proportion of phenotypic variance contributed by the selected SNP set ( h u 1 2 ) was computed using the following equation: h u 1 2 = σ u 1 2 σ u 1 2 + σ u 2 2 + σ e 2 .

To identify positional candidate genes, the BioMart database (http://www.ensembl.org/) was implemented. The candidate genes resided within the genomic regions of up- and downstream 500 kb around the significant SNPs were taken into account in our study. Functional annotation of the genes located in the regions of interest was performed with R package WebGestaltR (Wang J et al., 2017).

Table 1 summarized the descriptive statistics of AGE100 and BF100 in the two Large White pig lines. Phenotypes of both traits followed a normal distribution (Supplementary Figure S1). PCA showed a substantial genetic diversity between these two populations. The total genetic variance in these animals was explained 14.92%, 1.25% and 1.16% by the first three principal components, respectively (Figure 1A). As seen in Figure 1A, PC1 distinctly divided them into Canadian and French line, suggesting that two pig populations shared less similarity in the genetic background. Moreover, in PC2 Canadian line had a widespread. This finding agreed with the results of model-based analysis of population admixture which was showed in Figure 1B (K = 2, 3 and 4). Figure 1C displayed average LD (r ² ) at various physical distances between two loci on all the autosomes. The average r ² at pair-wise SNP distance of less than 5 Mb on autosomes ranged from 0.432 to 0.582 in Canadian Large White population. The LD decay pattern in French population was similar to the Canadian line with average r ² ranging from 0.431 to 0.588. The average r ² decreased much more slowly with the increase of pair-wise SNP physical distance and remained constant beyond 1 Mb in two lines (Figure 1C). Ne estimated at 99 generations ago were 110 for Canadian line and 107 for French line, respectively (Figure 1D). It should be noted that MAFs in Canadian and French lines had no significant differences (Supplementary Figure S2).

To identify the selection signatures between these two populations, F_ST, hapFLK and ROH were used to detect selection signatures across the whole genome in two pig lines (Supplementary Figure S3). In each selection signature method, the top 20% of genomic regions were selected. All the genes in the genomic regions under selection were further analyzed for functional annotation using the WebGestaltR package (Wang J et al., 2017) (Supplementary Figure S4). To find the proportions of phenotypic variances explained by SNPs subjected to selection, the entire genome was divided into five groups according to selection signatures. The heritability of AGE100 and BF100 explained by these groups were jointly estimated using a multi-components (a GMR for each group) linear mixed model. However, we did not observe heritability of AGE100 and BF100 tended to enrich in the regions under selection (Supplementary Figure S5). The possible reason for this result was the lack of statistical power due to the small population size and the traits are complex (Ma et al., 2019).

For both lines, no genome-wide significant SNPs were detected associated with AGE100 trait (Figure 2A, Figure 2B), while only four SNPs at the suggestive significant level were observed (Figure 2C; Table 2). Only one SNP located in Sus scrofa chromosome 7 (SSC7) was identified by cross-population GWAS to be significantly associated with AGE100 at the suggestive significant level, which was also detected in the combined-population (Figure 2D; Table 2).

For Canadian line, two SNPs (SSC2:236157, SSC2:3285386) were significantly associated with BF100 at genome-wide significant level. These two SNPs were also detected by the cross-population GWAS analysis, which resided in the genic regions of ANO9 and ANO1 (Figure 3A; Table 3). There was no significant SNPs identified to be associated with BF100 trait in the French line (Figure 3B).

Three SNPs (SSC1:159697691, SSC2:3285386 and SSC18:48235741) were identified to be significantly associated with BF100 using the combined-population GWAS at the suggestive significant level (Figure 3C; Table 3). Cross-population GWAS analysis of association results from the two pig populations revealed two additional SNPs for BF100, and the most significantly associated SNP from the cross-population GWAS was located on chromosome 2 (Figure 3D; Table 3).