Screening of Litter-Size-Associated SNPs in NOX4, PDE11A and GHR Genes of Sheep

Simple Summary NOX4, PDE11A and GHR have previously been screened as candidate genes for litter size in sheep. Therefore, it is necessary to validate these genes in the sheep population and identify loci that are associated with litter size. In this study, the results of association analysis showed that c.1057-4C > T in NOX4 and c.1983C > T in PDE11A were significantly associated with litter size in Small Tail Han sheep, and their effects were independent of each other. In summary, this study provided two new genetic markers for improving litter size in sheep. Abstract In previous studies, NOX4, PDE11A and GHR genes have been screened as important candidate genes for litter size in sheep by using the GWAS method; however, neither their effects on litter size nor the loci associated with litter size have been identified. In this study, three candidate loci (c.1057-4C > T in NOX4, c.1983C > T in PDE11A and c.1618C > T in GHR) were first screened based on our previous resequencing data of 10 sheep breeds. After the three loci were genotyped using Sequenom MassARRAY technology, we carried out population genetics analysis on the three loci and performed association analysis between the polymorphism of the three loci and the litter size of sheep. The results of population genetics analysis suggested that c.1057-4C > T in NOX4 and c.1983C > T in PDE11A may be subject to natural or artificial selection. The results of association analysis indicated that litter size was significantly associated with c.1057-4C > T in NOX4 and c.1983C > T in PDE11A (p < 0.05) in Small Tail Han sheep, and there was no significant interaction effect between the two loci on the litter size. In summary, c.1057-4C > T in NOX4 and c.1983 C > T in PDE11A can be considered candidate molecular markers for improving litter size in sheep.


Introduction
Increasing litter size is one of the most effective ways to improve the economic benefits of sheep farms.However, as litter size is a trait with low heritability, it is difficult to improve in a short time period through conventional breeding methods.The marker selection breeding method is more effective for litter size improvement.Nevertheless, very few of the major genes and markers for litter size are currently known, so it is essential to explore new genetic markers for accelerating the genetic progress of litter size in sheep.
Through resequencing the genome of 248 sheep from 36 local breeds and 6 improved breeds around the world, Li et al. [1] screened reproductive candidate genes such as the reduced nicotinamide adenine dinucleotide phosphate oxidase 4 gene (NOX4) and the 11th member of the phosphodiesterase family gene (PDE11A) by simultaneously using genome-wide association study (GWAS) analysis and selection signal analysis.It was found that NOX4 is highly expressed in the follicles of rats [2].Activated NOX4 can produce excess reactive oxygen species (ROS) and lead to sustained oxidative stress, impaired oocyte quality or follicular atresia in the ovary [3,4].Furthermore, ROS was involved in the transduction of growth factors and played a key role in the ovulatory signal cascade of rodents [5], suggesting that ROS is an important regulatory factor during the ovulation process of animals.The PDE11A gene can catalyze the hydrolysis of cAMP [6].In the animal ovary, cAMP acts as the second messenger for follicle-stimulating hormone (FSH) and luteinizing hormone (LH) receptors [7,8], and plays a central role in all stages of follicular development and ovulation [9].Similarly, through GWAS analysis of five prolific and one low-fertility sheep breed, Xu et al. [10] screened the growth hormone receptor (GHR) as a candidate gene for litter size in sheep.GHR was also highly expressed in reproductiverelated tissues such as the ovary and uterus.The GHR is located on the surface of early follicles, which can directly regulate the formation of late follicles by binding to the growth hormone (GH) [11].Bachelot et al. found that the GHR gene is involved in the regulation of estrogen secretion, follicular growth and development, ovulation number and litter size in mice [12].In addition, there were fewer follicles in GHR-knockout mice than in wild mice [11].These results suggested that the above three candidate genes were closely related to follicular development and ovulation; however, the effect of these three genes on litter size in sheep remains unclear.

Animal Preparation and Sample Collection
All the sheep breeds that were analyzed in this experiment were divided into a polytocous group and a monotocous group.The polytocous group included Hu sheep (n = 96, Xuzhou City, Jiangsu Province, China), Cele Black sheep (n = 96, Qira County, Hotan Prefecture, Xinjiang Uygur Autonomous Region, China) and Small Tail Han sheep (n = 384, Heze City, Shandong Province, China).The monotocous group included Sunit sheep (n = 96, Urad Front Banner, Bayannaoer City, Inner Mongolia Autonomous Region, China) and Bamei mutton sheep (n = 96, Linhe District, Bayannaoer, Inner Mongolia Autonomous Region, China).The litter size information was recorded in detail in three parities of Small Tail Han sheep.However, it was difficult to obtain litter size records of all three parities in other sheep populations.Jugular vein blood was collected from all these ewes.

DNA Extraction
DNA was extracted from blood samples of the above sheep using a DNA extraction Kit (Beijing, China, Tiangen Biochemical Technology Co., Ltd., DP304-03).The concentration of DNA samples was determined by using a NanoDrop 2000 (Thermo Scientific, Madison, WI, USA), and the quality of DNA was detected through 1.5% agarose gel electrophoresis.

Selection of Candidate Loci in Three Genes
These three genes (NOX4, PDE11A and GHR) were selected as candidate genes to analyze their association with sheep litter size.First, based on our previous resequencing data of 10 sheep breeds [13], the polymorphic loci in these three genes were screened between the high-fecundity group (Hu sheep, Cele Black sheep, Small Tail Han sheep and Australian Merino sheep) and the low-fecundity group (Sunit sheep, Bamei mutton sheep, Tan sheep, Prairie Tibetan sheep, Valley Tibetan sheep, Oula sheep and Bayinbuluke sheep).Then, the allele frequency of each polymorphic locus in these 10 breeds and the population differentiation coefficient (F st ) of each locus between the two groups were calculated.Next, the three candidate loci (c.1057-4C > T in NOX4, c.1983C > T in PDE11A and c.1618C > T in GHR) were screened based on their potential effects on the sequence and structure of protein and their F st values (their F st values need to meet any one of the following conditions: 1 ⃝ non-synonymous mutations or SNPs located in a splice region or regulatory region with F st > 0.05; 2  ⃝ mutations in intron with F st > 0.15).

Statistical Analysis
The allele and genotype frequency, polymorphism information content (PIC), heterozygosity (He) and number of effective alleles (Ne) were calculated using the following formulae: where n is the number of alleles, p i is the allele frequency of the ith allele and p j is the allele frequency of the jth allele.The chi-square test was used to detect whether the genotype distribution of each locus deviated from Hardy-Weinberg equilibrium.The association of litter size with the genotypes of three loci was analyzed using the linear model: the least-squares means for multiple comparisons of litter size among the different genotypes in Small Tail Han sheep.y is the phenotypic value of litter size, µ is the population mean, P is the fixed parity effect, G1 is the fixed effect for candidate SNP in NOX4, G2 is the fixed effect for candidate SNP in PDE11A, G3 is the fixed effect for candidate SNP in GHR, and e is the random error effect of each observation.Meanwhile, we also examined the interaction effect between each two loci, in order to check whether the interaction between candidate loci needs to be supplemented in the model.However, we did not find a significant interaction between the candidate loci.Analysis was performed using R software (avo, Version 4.0.3).

Protein Functional Domain and Interaction Networks Analysis
For SNP loci that were significantly associated with litter size, we further analyzed the location of the SNPs in the gene and their potential influence on nearby protein functional domains.In addition, in order to understand what kind of interaction network, including the candidate proteins, regulates follicular development, the interaction networks of the candidate proteins were predicted using the STRING database v.12.0 (https://cn.string-db.org/, accessed on 20 November 2023) [16].

Genotyping and Population Genetic Analysis of Candidate SNPs in NOX4, PDE11A and GHR Genes
The candidate SNPs (c.1057-4C > T in NOX4, c.1983C > T in PDE11A and c.1618C > T in GHR) were genotyped in all ewe samples, and the results are shown in Figure 1.It can be seen that there were three genotypes in each SNP locus, and three genotypes could be well distinguished by the MassARRAY system.
the interaction effect between each two loci, in order to check whether the interaction between candidate loci needs to be supplemented in the model.However, we did not find a significant interaction between the candidate loci.Analysis was performed using R software (avo, Version 4.0.3).

Protein Functional Domain and Interaction Networks Analysis
For SNP loci that were significantly associated with litter size, we further analyzed the location of the SNPs in the gene and their potential influence on nearby protein functional domains.In addition, in order to understand what kind of interaction network, including the candidate proteins, regulates follicular development, the interaction networks of the candidate proteins were predicted using the STRING database v.12.0 (https://cn.string-db.org/,accessed on 20 November 2023) [16].

Genotyping and Population Genetic Analysis of Candidate SNPs in NOX4, PDE11A and GHR Genes
The candidate SNPs (c.1057-4C > T in NOX4, c.1983C > T in PDE11A and c.1618C > T in GHR) were genotyped in all ewe samples, and the results are shown in Figure 1.It can be seen that there were three genotypes in each SNP locus, and three genotypes could be well distinguished by the MassARRAY system.The results of population genetic analysis are summarized in Table 2  The results of population genetic analysis are summarized in Table 2  The association between the three SNP loci in NOX4, PDE11A and GHR genes with the litter size was analyzed in Small Tail Han sheep (Table 3).It is worth noting that there was a significant association between the polymorphism of the NOX4 c.1057-4C > T locus and the litter size in the three parities of Small Tail Han sheep (p < 0.05).The litter size of ewes with the TT-genotype was significantly higher than that of ewes with the CT-genotype and CC-genotype (p < 0.05) in all three parities, and the litter size of CT-genotype ewes was also significantly higher than that of CC-genotype ewes (p < 0.05) in the first parity.For the PDE11A c.1983C > T locus, there was also a significant association between its polymorphism and the litter size of Small Tail Han sheep (p < 0.05).The litter size of TT-genotype ewes was significantly higher than those of CT-and CC-genotype ewes in all three parities (p < 0.05), and the litter size of CT-genotype ewes was also significantly higher than that of CC-genotype ewes in the second and third parities (p < 0.05).However, there was no significant association between GHR c.1618C > T locus and litter size in the Small Tail Han population (p > 0.05).Additionally, parity had a significant effect on the litter size of Small Tail Han sheep (p = 8.75 × 10 −12 < 0.01).NOTE: The number in brackets denotes the number of ewes with the kind of genotype.Different lower-case letters next to the litter size mean that there are significant differences between genotypes (p < 0.05).

Protein Structure and Interaction Network Analysis of NOX4 and PDE11A
For the two SNPs (c.1057-4C > T and c.1983C > T) that were significantly associated with litter size, we further investigated the potential effect of the mutations on protein structure (Figure 2).Both SNPs belong to the splice polypyrimidine tract variant, which may lead to selective splicing.The SNP c.1057-4C > T is located at the posterior boundary of intron 10 in NOX4 (Figure 2B), which might affect the splicing, composition and function of the domain including exon 11.Exon 11 of the NOX4 gene is located in the FAD-binding_6 Domain (Figure 2A), which is the core catalytic part of ROS induction [17] and an important regulatory factor in the ovulation process of animals [18,19].The SNP c.1983C > T is located at the last position of exon 12 in PDE11A (Figure 2D), which may affect the composition and function of the domain including exon 12. Exon 12 of the PDE11A gene is located in the PDEase_1 Domain, and this domain participates in catalyzing the hydrolysis of cAMP [20] and plays a central role in all stages of follicular development and ovulation [21].In order to further understand the biological functions of NOX4 and PDE11A genes in the reproduction process of sheep, the protein networks that closely interact with them were constructed.Among the 10 strongest interacting proteins with NOX4 (Figure 3A,B), 3 (TLR4, NCF2 and CYBB) are crucial for follicular development as well as ovulation [22][23][24][25].Among the 10 proteins that closely interact with PDE11A, 4 (AK3, APRT, ENPP3 and In order to further understand the biological functions of NOX4 and PDE11A genes in the reproduction process of sheep, the protein networks that closely interact with them were constructed.Among the 10 strongest interacting proteins with NOX4 (Figure 3A,B), 3 (TLR4, NCF2 and CYBB) are crucial for follicular development as well as ovulation [22][23][24][25].Among the 10 proteins that closely interact with PDE11A, 4 (AK3, APRT, ENPP3 and ENTPD1) are closely associated with oocyte development and ovulation [26][27][28].

Discussion
Through GWAS analysis and selection signal sweeps, Li et al. [1] and Xu et al. [10] screened key reproductive candidate genes, including NOX4, PDE11A and GHR.Therefore, these three genes are notable candidate genes, and their association with litter size should be studied further.In this study, based on our previous resequencing data of 10 sheep breeds, we screened three SNP loci in these three genes and further analyzed their association with litter size.
Population genetic analysis showed that c.1983 C > T in the PDE11A gene deviated from Hardy-Weinberg equilibrium in Small Tail Han sheep and Sunit sheep, and c.1057-4C > T in the NOX4 gene deviated from Hardy-Weinberg equilibrium in Hu sheep, Cele Black sheep and Sunit sheep.The results suggest that both SNP loci might be under the influence of natural or artificial selection within these sheep breeds.NOX4 gene expression regulates the production of ROS [29], whose excess generation can increase the incidence of diseases such as preeclampsia and lead to placental abruption or stillbirths [30].Coincidentally, mutations in the PDE11A gene have also been found to be associated with preeclampsia, premature birth and stillbirths [31].Thus, it is plausible that during the domestication and breeding processes of these sheep populations, these two SNP loci undergo strong artificial selection, aiming to minimize the occurrences of stillbirths and preterm births.Conversely, the SNP locus in the GHR gene maintained Hardy-Weinberg equilibrium across the five sheep breeds, suggesting that this particular locus may not have been influenced by either natural or artificial selection.
The association analysis demonstrated a significant relationship between the c.1057-4 C > T locus in the NOX4 gene and litter size in Small Tail Han sheep.NOX4, a primary source of ROS production, is crucial for follicular maturation.Briefly, NOX4 plays an important role in the survival of granulosa cells and follicular development.However, the association between litter size and mutations in the gene has not been studied in sheep.The NOX4 protein consists of 594 amino acids that are encoded by 18 exons and contains three functional domains: Ferric reduct Family Domain (amino acids 73 to 220), FAD_binding_8 Domain (amino acids 321-432) and Ferric reductase NAD binding Domain (amino acids 438-577).Exon 11 is located in the FAD_binding_8 Domain, which is essential for electron transfer [32] and a central catalytic part of ROS production induction

Discussion
Through GWAS analysis and selection signal sweeps, Li et al. [1] and Xu et al. [10] screened key reproductive candidate genes, including NOX4, PDE11A and GHR.Therefore, these three genes are notable candidate genes, and their association with litter size should be studied further.In this study, based on our previous resequencing data of 10 sheep breeds, we screened three SNP loci in these three genes and further analyzed their association with litter size.
Population genetic analysis showed that c.1983 C > T in the PDE11A gene deviated from Hardy-Weinberg equilibrium in Small Tail Han sheep and Sunit sheep, and c.1057-4C > T in the NOX4 gene deviated from Hardy-Weinberg equilibrium in Hu sheep, Cele Black sheep and Sunit sheep.The results suggest that both SNP loci might be under the influence of natural or artificial selection within these sheep breeds.NOX4 gene expression regulates the production of ROS [29], whose excess generation can increase the incidence of diseases such as preeclampsia and lead to placental abruption or stillbirths [30].Coincidentally, mutations in the PDE11A gene have also been found to be associated with preeclampsia, premature birth and stillbirths [31].Thus, it is plausible that during the domestication and breeding processes of these sheep populations, these two SNP loci undergo strong artificial selection, aiming to minimize the occurrences of stillbirths and preterm births.Conversely, the SNP locus in the GHR gene maintained Hardy-Weinberg equilibrium across the five sheep breeds, suggesting that this particular locus may not have been influenced by either natural or artificial selection.
The association analysis demonstrated a significant relationship between the c.1057-4 C > T locus in the NOX4 gene and litter size in Small Tail Han sheep.NOX4, a primary source of ROS production, is crucial for follicular maturation.Briefly, NOX4 plays an important role in the survival of granulosa cells and follicular development.However, the association between litter size and mutations in the gene has not been studied in sheep.The NOX4 protein consists of 594 amino acids that are encoded by 18 exons and contains three functional domains: Ferric reduct Family Domain (amino acids 73 to 220), FAD_binding_8 Domain (amino acids 321-432) and Ferric reductase NAD binding Domain (amino acids 438-577).Exon 11 is located in the FAD_binding_8 Domain, which is essential for electron transfer [32] and a central catalytic part of ROS production induction [17].Unbalanced ROS levels may lead to oxidative stress and a significant decline in the number and quality of oocytes [18].In addition, Gnainsky et al. found that reducing FAD levels can inhibit mitochondrial activity and oogenesis in oocytes [19].It is worth noting that the c.1057-4 C > T mutation in NOX4 is located at the anterior boundary adjacent to exon 11 and belongs to the splice polypyrimidine tract variant, which may affect the function of the FAD_binding_8 Domain through selective splicing and may change the concentration of ROS.Thus, follicular development, oocyte activity and oogenesis may ultimately be affected due to changes in ROS.
Similarly, according to the association analysis, c.1983C > T in PDE11A was significantly associated with litter size in Small Tail Han sheep.In mammalian oocytes, the meiotic cycle is influenced by the levels of cAMP and cGMP [33][34][35], which enhances the survival ability of granulosa cells [21].The activity of cAMP-PDE is mainly provided by PDE11A and PDE8A [36].In addition, mutations in PDE11A will affect the number of meiosis occurrences and oocyte ovulation [21].However, we also did not find any study on the association between PDE11A gene mutations and litter size in animals.The PDE11A protein is encoded by 914 amino acids across 20 exons and consists of three functional domains: two GAF Domains (amino acids 195-348 and amino acids 380-536) and one PDEase_1 Domain (amino acids 643-878).Exon 12 is located in the PDEase_1 Domain, which is associated with the ability to catalyze cAMP hydrolysis [20].The c.1983 C > T mutation in PDE11A is located at the last base of exon 12, which belongs to the splice region variant.Thus, it may affect the structure of the PDEase_1 Domain by alternative splicing and change the ability to catalyze the hydrolysis of cAMP.Ultimately, the above changes might affect oocyte meiosis, the survival ability of granulosa cells and the number of ovulations.
Based on the protein interaction network provided by the STRING database, we identified the 10 proteins that closely interact with NOX4.Among them, three proteins (TLR4, NCF2 and CYBB) were crucial for follicular development and ovulation.TLR4 is highly expressed in mouse cumulus cells and is involved in the proliferation of cumulus-oocyte complexes and ovulation.Moreover, its deficiency can affect ovulation and pregnancy rates [37].Additionally, its level will increase in the induced PCOS mice [38].The activity of NCF2 may lead to increased ROS accumulation, thereby reducing the development ability of oocytes [39,40].Furthermore, Kuokkanen et al. found that elevated levels of CYBB and NCF2 mRNA may potentially harm luteal cells [25].Among the 10 proteins that are closely related to PDE11A, 4 (AK3, APRT, ENPP3, ENTPD1) play important roles in different periods from oocyte development to embryonic development.The expression of AK3 in oocytes indirectly regulates ATP levels and thus influences the activity of oocytes [41].The activity of APRT varies across different periods of ovulation and embryonic development [26].ENPP3 was identified as a new biological marker for tubal metaplasia [27] and is very important for morphological changes and inflammatory responses during ovulation and luteinization [42].Kauffenstein et al. [43] found that the mice carrying a mutation in ENTPD1 had a smaller litter size compared to the wild type.In summary, these networks further confirm the relationship between candidate genes and follicular development and litter size.However, the interplay among these genes and how they synergistically regulate follicular development and ovulation require further exploration.

Conclusions
The results of association analysis indicated that litter size was significantly associated with c.1057-4C > T in NOX4 and c.1983C > T in PDE11A (p < 0.05) in Small Tail Han sheep, and there was no significant interaction effect between the two loci on the litter size.Therefore, these two loci can be considered candidate molecular markers for improving litter size in sheep.
. The c.1057-4C > T locus of the NOX4 gene displayed a moderate polymorphism in Bamei mutton sheep (0.25 ≤ PIC < 0.5) and a low polymorphism in Small Tail Han sheep, Hu sheep, Cele Black sheep and Sunit sheep (PIC < 0.25).The c.1983 C > T locus of the PDE11A gene showed a moderate polymorphism (0.25 ≤ PIC < 0.5) in Hu sheep, Cele Black sheep and Bamei mutton sheep, and a low polymorphism (PIC < 0.25) in Small Tail Han sheep and Sunit sheep.The c.1618C > T locus of the GHR gene presented a moderate polymorphism (0.25 ≤ PIC < 0.5) in Small Tail Han sheep, Hu sheep, Sunit sheep and Bamei mutton sheep, and a low polymorphism (PIC < 0.25) in Cele Black sheep.According to the results of the chi-square test, the c.1057-4C > T locus of the NOX4 gene was in Hardy-Weinberg equilibrium (p > 0.05) in Small Tail Han sheep and Bamei mutton sheep only.The PDE11A gene c.1983 C > T locus was in Hardy-Weinberg equilibrium in Hu sheep, Cele Black sheep and Bamei mutton sheep (p > 0.05).The c.1618 C > T locus of the GHR gene was in Hardy-Weinberg equilibrium in all five sheep breeds (p > 0.05).
. The c.1057-4C > T locus of the NOX4 gene displayed a moderate polymorphism in Bamei mutton sheep (0.25 ≤ PIC < 0.5) and a low polymorphism in Small Tail Han sheep, Hu sheep, Cele Black sheep and Sunit sheep (PIC < 0.25).The c.1983 C > T locus of the PDE11A gene showed a moderate polymorphism (0.25 ≤ PIC < 0.5) in Hu sheep, Cele Black sheep and Bamei mutton sheep, and a low polymorphism (PIC < 0.25) in Small Tail Han sheep and Sunit sheep.The c.1618C > T locus of the GHR gene presented a moderate polymorphism (0.25 ≤ PIC < 0.5) in Small Tail Han sheep, Hu sheep, Sunit sheep and Bamei mutton sheep, and a low polymorphism (PIC < 0.25) in Cele Black sheep.According to the results of the chi-square test, the c.1057-4C > T locus of the NOX4 gene was in Hardy-Weinberg equilibrium (p > 0.05) in Small Tail Han sheep and Bamei mutton sheep only.The PDE11A gene c.1983 C > T locus was in Hardy-Weinberg equilibrium in Hu sheep, Cele Black sheep and Bamei mutton sheep (p > 0.05).The c.1618 C > T locus of the GHR gene was in Hardy-Weinberg equilibrium in all five sheep breeds (p > 0.05).

Animals 2024 , 12 Figure 2 .
Figure 2. The location and functional domain of NOX4 c.1057-4C > T and PDE11A c.1983 C > T. (A) Functional domain of NOX4 and the location of exon 11. (B) Location of c.1057-4C > T in the NOX4 gene of sheep.(C) Functional domain of PDE11A and the location of exon 12. (D) Location of c.1983C > T in the PDE11A gene of sheep.

Figure 2 .
Figure 2. The location and functional domain of NOX4 c.1057-4C > T and PDE11A c.1983 C > T. (A) Functional domain of NOX4 and the location of exon 11. (B) Location of c.1057-4C > T in the NOX4 gene of sheep.(C) Functional domain of PDE11A and the location of exon 12. (D) Location of c.1983C > T in the PDE11A gene of sheep.

Animals 2024 , 12 Figure 3 .
Figure 3. (A) Protein networks that closely interact with NOX4 in sheep, provided by the STRING database (version 12).(B) Protein networks that closely interact with PDE11A in sheep, provided by the STRING database (version 12).The proteins in the red box are closely related to ovulation and follicular development.

Figure 3 .
Figure 3. (A) Protein networks that closely interact with NOX4 in sheep, provided by the STRING database (version 12).(B) Protein networks that closely interact with PDE11A in sheep, provided by the STRING database (version 12).The proteins in the red box are closely related to ovulation and follicular development.

Table 1 .
The primer sequences for MassARRAY analysis.

Table 2 .
Population genetic analysis of candidate loci of NOX4, PDE11A and GHR genes in five sheep breeds.

Table 3 .
Least-squares means and standard deviations of litter size for different genotypes in Small Tail Han sheep.