A total of 3,072 pigs were used in the present study, and they were from two farms in Liaoning and Shanxi Provinces, China. The GL, breed, farm, parity, and total litter birth weight were collected for these sows, and 5,662 records were shown in Supplementary Table S1. The pedigrees of these sows were traced to four generations with 5,786 individuals to be used to estimate the genetic parameters of GL.

By the dmuai of DMU Version 6.5.1, we estimated the variance and covariance components of GL based on a repeatability model, which was described as follows: y = X b + Z a a + Z p p + e in which, y was the vector of GL phenotype; b was the vector of fixed effects including farm (two farms), breed (three variables: Yorkshire, Landrace, and Duroc pigs), pairty (four variables: 1, 2, 3, and 4–8 parity), and total litter birth weight (31 variables: 1–29, 32, and 33 kg); X was the design matrix associated b with y ; p was the vector of random permanent environmental effects; a was the vector of additive genetic effects; Z p and Z a referred to the corresponding residual effects; and e was the vector of random residual effects.

Further, we calculated the heritability and repeatability of GL by the following formulas: h 2 = σ a 2 σ a 2 + σ p 2 + σ e 2   a n d   r e 2 = σ a 2 + σ p 2 σ a 2 + σ p 2 + σ e 2 ; where h 2 and r e 2 were the heritability and repeatability, respectively. σ a 2 was the additive genetic variance, σ p 2 was the permanent environmental effect variance, and σ e 2 was the residual effect variance.

Of all 3,072 sows with GL phenotypes, 906 individuals were used for GWAS. The ear samples of these 906 sows were collected and the DNA was isolated from them by a commercially available kit (Q1Aamp DNA Mini Kit, QIAGEN, Germany). Then these pigs were genotyped with the GeneSeek Genomic Profiler (GGP) Porcine 50K Chip (Illumina, San Diego, CA, United States), which contained 50,697 SNPs. With the PLINK, we performed the quality control by removing the SNPs with Minor Allele Frequencies (MAF) < 0.05, call rate <95%, and a deviation from Hardy-Weinberg equilibrium (HWE, p < 0.000001). Finally, we obtained 34,609 SNPs and 906 animals for GWAS, in which, all the individuals had a high SNP detection rate (>90%). All SNP positions were annotated according to the pig genome assembly Sscrofa 11.1.

To improve the identification of QTL, the whole-genome sequence data is expected to use (van den Berg et al., 2019), while sequencing thousands of individuals is expensive. In this study, we conducted the imputation for the GGP Porcine 50K Chip after quality control according to the online software PHARP v1 (http://alphaindex.zju.edu.cn/PHARP/index.php), which included 2,012 haplotypes, 34 million SNPs of autosomes, and 71 pig breeds (Wang et al., 2022). After imputation, we carried out the quality control (R ² ≥ 0.8 and MAF = 0.05), and obtained 9,212,179 SNPs, in which, the average R ² was 0.9248. Further, we conducted the quality control and GWAS according to the filter and model used in the front GGP Porcine 50K Chip, where 9,043,409 SNPs and 906 individuals were utilized.

We conducted the GWAS by GCTA software with the model as follows: y = μ + b X + G + F + P + m × M + e where y is a vector of phenotype values for all individuals; b is the average effects of gene substitution of a particular SNP; X is a vector of the SNP genotype; G is the genomic relationship matrix created by all SNPs; F and P were the farm (two factors) and parity (four factors) fixed effects, respectively; M is the fixed effects of total litter birth weight; m is the regression coefficient of M and e is the random residual.

The results of GWAS were shown by Manhattan and quantile-quantile (Q-Q) plots. A significant SNP was determined if the raw p-value <0.05/n, in which, n was the number of SNPs. Further, the threshold of suggestive significant association was calculated according to 1/n.

In population genetics, Fst is an index to measure whether the actual frequency of genotypes in a population deviates from the theoretical proportion of genetic balance. In recent years, the number of litter sizes has increased, and the GL has been obvious for a long. To evaluate the degree of genetic differentiation of individuals with the different GL, Fst was carried out through the vcftools software. We calculated the Fst value of each SNP between GL 114 days and 115/116/117/118 days based on the following formulas: F S T = H T − H S H T , H T = 2 × 1 − P E + P F 2 × P E + P F 2 , and H S = 2 × P E × 1 − P E + 2 × P F × 1 − P F 2 ;

In which, P E and P F are the allele frequencies of individuals with GL 114 days and 115/116/117/118 days, respectively.

The SNPs involved in the top 1% (Fst values) were selected as the core SNPs associated with the GL trait.

By BioMart of Ensembl database, the candidate genes within 500 kb of significant SNPs were screened based on the pig reference genome (Sscrofa11.1). Further, we investigated the functions of these genes through the KOBAS (http://kobas.cbi.pku.edu.cn/kobas3/genelist/), and the significance threshold was corrected p-value < 0.05.

GL phenotype distribution and descriptive statistics were shown in Figure 1. The average GL was 115.97 days in this population, and the maximum and minimum values were 120 and 112, respectively. By estimation, the additive genetic variance V(a), permanent environmental effect variance V(p), and residual effect variance V(e) were detected, namely, 0.3014, 0.1594, and 1.4299, respectively. Correspondingly, the heritability and repeatability were calculated: h 2 = 0.1594 and r e 2 = 0.2437.

Through the GWAS based on the GGP Porcine 50K Chip, four genome-wide significant SNPs were identified (p < 1.44E-06 (0.05/34,609), Table 1 and Figures 2A,B), including WU_10.2_8_32052217 and WU_10.2_3_124292040 located on (Sus scrofa Chromosome (SSC) 3, and ALGA0073634 and WU_10.2_13_204954425 located on SSC13). According to the suggestive significant threshold p < 2.89E-05 (1/34,609), additional five SNPs were detected (Table 1 and Figures 2A,B), namely, WU_10.2_3_124371187 and WU_10.2_3_124478097 located on SSC3, H3GA0016351 and WU_10.2_5_111419379 located on SSC5, and CASI0009756 located on SSC13.

Further, the genotype was imputed by PHARP v1 to be used for a GWAS, and a total of 53 genome-wide significant SNPs [p < 5.52E-09 (0.05/9,043,409)] on pig chromosomes SSC13 were identified (Figures 2C,D; Table 2), including SNP ALGA0073634. Based on the suggestive significant threshold p < 1.11E-07 (1/9,043,409), 112 SNPs were detected to be associated with GL, containing 111 located on SSC13 and one located on SSC7 (Figures 2C,D; Table 2).

GLs of 114, 115, 116, 117, and 118 days were involved in 309, 263, 116, 117, and 49 individuals, respectively. Based on the GGP Porcine 50K Chip data, we calculated the Fst values of GL between 114 days and 115/116/117/118 days for each SNP, and the FST density distribution maps were constructed (Figure 3). The average Fst values of the four groups were 0.0041, 0.0279, 0.0341, and 0.0376, indicating that there was a differentiation between individuals with GL of 114 days and 115/116/117/118 days. In addition, the group of 114 and 118 days had the highest degree of differentiation than other groups.

For each group, 320 SNPs in the top 1% (FST > 0.0317, 0.1489, 0.1748, and 0.2108 for 114 and 115/116/117/118 days, respectively) were the core SNPs in the region where the signal occurred (Figure 3), a total of 885 core SNPs were identified (Supplementary Table S2). Out of them, SSC13 had the largest number of core SNPs with 119, followed by SSC7 (116 core SNPs) and SSC9 (89 core SNPs). The highest FST value (FST = 0.4035) was located on SSC1. According to the chromosomal location and FST value of the selection signal, the distribution of 885 core SNPs for each group was marked in Figure 4, which showed the differences in the distribution of selection signals related to GL across the genome. Combined with GWAS results, four SNPs, WU_10.2_8_32052217, WU_10.2_3_124292040, WU_10.2_3_124371187 and WU_10.2_3_124478097 located on SSC3, were co-identified by GWAS and Fst analysis.

The genes located within a 500 kb distance on either side of the significant SNPs were considered candidates for GL. In total, 4,588 candidate genes were detected, of which, 27 were identified by GWAS of GGP Porcine 50K Chip data (Table 1), 12 were identified by GWAS of imputation SNPs (Table 2), and 4,574 were screened by core SNPs of Fst analysis (Supplementary Table S2). There were 20 genes simultaneously confirmed by GWAS and Fst analysis, namely HUNK, ARHGDIB, ERP27, RERG, NEDD9, TMEM170B, SCAF4, SOD1, TIAM1, ENSSSCG00000048838, ENSSSCG00000047227, EDN1, HIVEP1, ENSSSCG00000043944, LRATD1, ENSSSCG00000048577, ENSSSCG00000042932, ENSSSCG00000041405, ENSSSCG00000045589, and ADTRP.

To research the functions of candidate genes, we performed functional annotation by KOBAS. 12 genes from GWAS of GGP Porcine 50K Chip data were significantly involved in four KEGG pathways and 84 GO terms, including protein and lipid metabolism, regulation of multicellular organism growth, embryo implantation, and ovarian follicle development (Supplementary Table S3). Seven candidates from GWAS of imputation SNPs were strongly involved in 11 KEGG and 109 GO enrichments, which were mainly associated with body and organ signaling pathways, cellular metabolism, regulation of neuron apoptotic process, and heart development (Supplementary Table S3). 2,518 genes from Fst results heavily participated in 24 KEGG and 128 GO enrichments, containing cell differentiation, early endosome, and carbohydrate metabolic process etc. (Supplementary Table S3).