The animal experiments used in this study strictly complied with the guidelines of the Animal Welfare Committee of Guizhou University (EAE-GZU-2023-E015). A total of 142 healthy Xiangsu breeding sows under uniform feeding conditions were selected from the Guizhou University pig farm. The litter size traits, including the total number born (TNB), number born alive (NBA), and weaned piglet number were recorded per sow for the first and second litters (Supplementary material). Subsequently, blood samples were collected from the anterior vena cava of 142 pigs using a 5 mL EDTA anticoagulant tube and stored at −20°C. After the data collection of farrowing traits, three sows were randomly selected for slaughter, and their heart, liver, spleen, lung, kidney, longissimus dorsi muscle, and ovary tissue samples were collected and stored at −80°C in a refrigerator.

Genomic DNA was extracted from the anterior vena cava blood samples using a DNA extraction kit (DP348; Tiangen, Beijing, China) following the manufacturer’s guidelines. In addition, the total RNA was extracted from the heart, liver, spleen, lung, kidney, longissimus dorsi muscle, and ovary tissue samples using the TRIzol extraction kit (9109, TaKaRa, Dalian, China). The first strand of cDNA was synthesized via reverse transcription using a reverse transcription kit (A230; Kangrun, Beijing, China) according to the manufacturer’s interactions.

The primers for DNA and cDNA amplification were designed according to the pig NOBOX DNA (accession number: NC_010451.4) and RNA (accession number: NM_001195116.1) sequences published on NCBI using Primer Premier 5.0 software (PREMIER Biosoft International, Palo Alto, CA, United States). Primers were synthesized by Qingdao Biotechnology Co., Ltd. (Chongqing, China), and the primer information is shown in Supplementary Table S1.

The PCR amplification of the genomic DNA was performed in a total volume of 20 μL, including 10 μL of 2× Taq PCR Master Mix (GeneStar, Beijing, China), 1 μL of forward primer, 1 μL of reverse primer, 1 μL of genomic DNA, and 7 μL of double-distilled water (ddH2O). The amplification conditions consisted of 35 cycles of 3 min of pre-denaturation at 94°C, 30 s of denaturation at 94°C, 30 s of denaturation at 63°C, 72°C annealing for 1 min, 72°C extension for 5 min, and preservation at 4°C. Next, 5 μL of amplification product was aspirated and subjected to 1% agarose gel electrophoresis for 25 min. The target bands were verified using a gel imaging system, and the amplicons were sent to Qingdao Biotech (Chongqing, China) for sequencing. The cDNA was amplified in a 20 μL q-PCR reaction system consisting of 1 μL of cDNA template, 0.5 μL of forward primer, 0.5 μL of reverse primer, 10 μL of 2× RealStar Fast SYBR qPCR Mix (GeneStar, Beijing, China), and 8 μL of ddH₂O. The cDNA amplification was replicated four times. The reaction conditions were pre-denaturation at 95°C for 2 min, denaturation at 95°C for 15 s, annealing at 60°C for 30 s, and extension at 72°C for 30 s for 40 cycles. GAPDH was the fluorescence quantitative reference gene.

NOBOX amino acid sequences from eight species, including Homo sapiens (XP_016867231.1), Mus musculus (NP_570939.1), Equus caballus (XP_ 023494578.1), Sus scrofa (NP_001182045.1), Ovis aries (XP_042104903.1), Bos taurus (XP_024846980.1), Cervus elaphus (XP_043727747.1), and Gallus gallus (XP_040516952.1) were obtained from NCBI,¹ and a phylogenetic tree was constructed using MEGA7² (17). Subsequently, the common motifs in the super secondary structure were predicted using the MEME tool³ to reveal the structural characteristics and functions of NOBOX proteins in the selected eight species (18). Besides, the multiple sequence comparison of the amino acid sequences of NOBOX was performed using Clustal Omega⁴ (19), and the results were uploaded to the online server Con Surf⁵ for sequence conservation evaluation. For the functional and stability studies of NOBOX proteins, PhD-SNP,⁶ SNAP2,⁷ I-Mutant2.0,⁸ and MuPro⁹ were used for prediction analysis (20–23). Additionally, Sopma¹⁰ and AlphaFold2 were used to predict the secondary structure and construct the 3D model of the tertiary structure of the protein (24, 25).

The presence of SNPs in NOBOX sequence was determined via peak plotting against the PCR sequencing reads using the SeqMan software (26). Wild-type and mutant sequences were aligned and compared using MegAlign and ClustalW software in the DNA Star package. Genotype and gene frequencies at each mutation locus were calculated directly. Hardy–Weinberg equilibrium (HWE) was evaluated using the chi-squared (χ²) test, and the gene polymorphism parameters included homozygosity (Ho), heterozygosity (He), number of effective alleles (Ne), and polypeptide information content (PIC) (27). Linkage disequilibrium (LD) and haplotype analyses among SNPs were performed using the SHEsis Main¹¹ software (28), and the degree of chain imbalance was evaluated using the r² value, where r² > 0.33 indicated a strong chain imbalance state (29). On the other hand, diplotypes were evaluated based on haplotypes.

The difference in the number of litters between the different genotype groups was compared using one-way analysis of variance in the average drop-down option in IBM SPSS Statistics 25 software. The following general linear model formula was used: Y_ijk = μ + G_i + S_k + A_j + e_ijk, where Y_ijk is the litter size and number of pigs weaned, μ is the mean, G_i is the fixed effect of genotype, S_k is the random effect of sire, A_j is fixed effect of age, and e_ijk is the residual effect. The results are presented as the mean ± standard error (30).

NOBOX expression levels at different mutation sites were calculated using the 2^−ΔΔCtd method (31), and its expression patterns were mapped using the GraphPad Prism 8 software. Data are expressed as mean ± standard deviation at two decimal places.

The tissue expression profiles of NOBOX are shown in Figure 1A NOBOX was expressed in the heart, liver, kidney, and ovary, with the highest expression level in the ovarian tissues, which was significantly higher than that in the heart, liver, and kidney (p < 0.01). On the contrary, there were no significant differences in NOBOX expression levels among the heart, liver, and kidney (p > 0.05). In addition, NOBOX was not expressed in the spleen, lung, and longest dorsal muscle.

The gel electrophoresis imaging is shown in Figure 1B. The PCR amplification products were consistent with the target fragment size, with clear and single bands, non-specific amplification, and no obvious trailing phenomenon, implying the primers had good specificity.

Sequence alignment between the sequencing results and the reference sequence of porcine NOBOX (Accession number: NC_010451.4) revealed six SNPs in NOBOX, which were labelled g.1624 T>C, g.1858 G>A, g.2770 G>A, g.2821 A>G, g.5659 A>G. Notably, two alleles and three genotypes were present in all the six SNPs (Figure 2). The wild-type and mutant sequence alignment revealed that the base G at the g.1858 G>A locus was mutated to A, altering the codon-GUG- to -AUG-. Consequently, methionine (M) replaced valine (V), and g.1858 G>A was a missense SNP. The mutation at g.5659 A>G changed the codon-CCA- to -CCG-, resulting in a synonymous mutation because-CCA- and-CCG- are simple codons and proline (P) was not replaced (Supplementary Figure S1).

The phylogenetic tree of NOBOX sequences from the eight species is presented in Figure 3. Among the NOBOX species affinities, pig (S. scrofa) was most closely related to human (H. sapiens), followed by house mouse (M. musculus) and horse (E. caballus), and was furthest removed from the chicken (G. gallus). Fifteen significant amino acid sequences were detected in the eight species, indicating functional similarity at the super-secondary structure (Supplementary Figure S2). In addition, NOBOX was poorly conserved across species, with p.V82M as the nonconserved mutation (Figure 4).

The predicted effect of missense SNPs on protein function based on the online prediction tool PhD-SNP yielded a score of 8, indicating a neutral effect. However, the prediction tool SNAP2 yielded a score of 37, suggesting altered protein function. Additionally, the prediction of protein stability using I-Mutant 2.0 and MuPro revealed the free energy changes of −1.21 and −0.53, respectively, with the P.V82M mutation reducing protein stability. A comparison of the secondary structure prediction results between the wild type and mutant revealed that p.V82M mutation increased the α helix and β turn and decreased the random coil (Supplementary Tables S2–S4).

The 3D model of the NOBOX protein at p.V82M constructed using AlphaFold2 is presented in Figure 5. The p.V82M mutation replaced nonpolar, uncharged valine with a large, nonpolar, uncharged methionine, which altered the polar interactions with surrounding amino acids and affected the protein structure and function after the mutation.

According to NOBOX population genetic analyses (Table 1), the purity of each SNP locus was greater than the heterozygosity and the effective allele number ranged from 0.12 to 0.37. Besides, the polypeptide of SNP loci g.1624 T>C, g.1858 G>A, g.2821 A>G, g.5659 A>G, and g.6025 T>A was at the moderate polymorphism level (0.25 < PIC < 0.50). These five loci were in the Shangsu hybridization annotation, with strong selection potential and rich genetic diversity. The χ² test revealed that all five SNP loci were in HWE (p > 0.05). However, the polypeptide content at locus g.2770 G>A was 0.12, a lower polymorphism level (PIC <0.25), contrary to the HWE based on the χ² test (p < 0.05).

The LD analysis of NOBOX SNPs using D′ and r² tests is shown in Figure 6 (32). LD analysis revealed that the D′ values ranged from 0.08 to 1.00 and r² values from 0 to 0.874. SNP loci between g.1624 T>C and g.5659 A>G, g.1624 T>C and g.6025 T>A, g.2821 A>G and g.5659 A>G, and g.5659 A>G, and g.6025 T>A, with the r² of 0.44, 0.52, 0.37, and 0.87, respectively, were n strong chain imbalance, with the strongest degree of chain imbalance between g.5659 A>G and g.6025 T>A (Table 2).

Haplotype analysis identified three dominant haplotypes with frequencies greater than 5% from the Xiangsu pig population. They included Hap1 (TAGAAT), Hap2 (TGGAAT), and Hap3 (CGGGGA), with the frequencies of 36.20, 18.30, and 18.10%, respectively (Table 3).

The correlation between the SNPs of pig NOBOX and litter size traits is shown in Table 4. The GG genotypes were significantly higher than AA genotypes at the g.1858 G>A locus and GG genotypes were significantly higher than AG genotypes at the g.5659 A>G locus in first-born TNB and NBA (p < 0.05). For g.6025 T>A, the TT and TA genotypes were significantly higher than AA genotype in the TNB and NBA of second-born sows (p < 0.05).

The association analysis between diplotype and litter size traits is summarized in Table 5. The TNB and number of weaned piglets in the first-born sows were significantly higher in diploid Hap1/1 (TTAAGGAAAATT) and Hap1/3 (TCGAGGAGAGTA) than in Hap1/2 (TTGAGGAAAATT) (p < 0.05). Diplotypes with <5.0% frequency were not considered.