Abstract
Extreme environmental conditions are a major challenge for livestock production. Changes in climate conditions, especially those that lead to extreme weather, can reduce livestock production. The screening of genes and molecular markers is of great significance to explore the genetic mechanism of sheep prolificacy traits in Taklimakan Desert environment. We selected healthy adult Pishan Red Sheep (PRS) and Qira Black Sheep (QR) which live in Taklimakan Desert environment, collected blood from jugular vein, extracted DNA, and prepared Illumina Ovine SNP50 chip. For PRS, linkage disequilibrium (LD) was calculated using the ovine SNP50 Beadchip and the effective population size (Ne) was estimated using SMC++. The genetic characteristics of PRS were analyzed by integrated haplotype score (iHS) and fixation index (FST). The result showed that r2 of PRS was 0.233 ± 0.280 in the range of 0–10 Kb and decreased with increasing distances. SMC++ tested that the Ne of PRS remained at 236.99 in recent generations. 184 genes were screened out under iHS 1% threshold, and 1148 genes were screened out with FST under the 5% threshold, and 29 genes were obtained from the intersection of the two gene sets. In this study, the genetic characteristics of PRS and QR were compared by ovine genome chip, and the related excellent genes were searched, providing reference for the protection of sheep germplasm resources and molecular breeding in a desert environment.
1 Introduction
Sheep is one of the earliest domesticated animals in the world, and also one of the most successful animals domesticated by human beings in the Neolithic age. After long-term domestication and different environments, sheep have great changes in morphology, physiology and behavior. Pishan Red sheep (PRS) living in the Taklimakan Desert is characterized by perennial estrus and multiple fetuses after long-term selection by nature and people. In addition to its high-quality production traits, PRS is also a rare breed of sheep because its origin is on the southern edge of the Taklimakan Desert and north of the Karakoram Mountains (). Pishan red sheep is a local sheep breed formed under the local cultural and geographical conditions. The origin and formation history of PRS is not fully understood. Therefore, to dig the genetic structure and molecular markers of PRS can better protect PRS.
Linkage disequilibrium (LD) can improve the accuracy of genomic association analysis and predict marker regions (). LD decay patterns also provide information about the evolutionary history of the population and can be used to estimate ancestral effective population size (Ne) (). Ne and other genetic events can also influence the extent of LD in the population (). Therefore, LD helps in understanding the selection patterns experienced by individual breeds. Currently, LD estimates have been reported in several studies for a variety of livestock species, such as cattle (), pigs (), horses (), chickens () and sheep ().
The selected regions were searched through different chromosomes to provide a molecular genetic basis for sheep protection. Compared with traditional selection methods, genomics can be evaluated early with higher accuracy (). Voight and Kudaravalli proposed an integrated haplotype score (iHS) test based on extended haplotype homozygosity (). The incomparability of test statistics caused by differences in recombination rates between different chromosome segments was corrected by calculating the Extended Haplotype Homozygosity (EHH) statistics and genetic distance integration. Fixation index (FST) is used to measure the degree of population differentiation, indicating that there are obvious allele frequency differences between populations (). Genetic drift and selection process can usually cause the genetic differentiation between populations. This method is suitable for selection signal detection of multiple populations. Now, it is necessary to selectively intervene in the breeding work of PRS through genomic selection technology to guarantee a better inheritance of the excellent production traits of PRS to future generations.
In order to analyze the genetic mechanism of prolific traits in PRS under desert environment, we selected prolific PRS and singleton pregnancy Qira Black Sheep (QR) under relevant survival background as research objects. Genetic basis of PRS was analyzed based on LD and Ne and the genetic mechanism of prolific traits in Taklimakan Desert was explored by using genomic selection method.
2 Materials and methods
2.1 Animal care
This work was conducted in accordance with the specifications of the Ethics Committee of Tarim University of Science and Technology (SYXK 2020-009).
2.2 Animal collection
33 PRS (polyembryony) and 40 QR (singleton pregnancy) were randomly selected from Pishan County and Cele County in Hetian region and all of them were healthy adult sheep with no genetic relationship. Blood samples were collected from the jugular vein and DNA was extracted with a DNA kit (Tiangen Biotech Co. Ltd., Beijing, China). The samples were sent to Beijing Compass Agritechnology Co., Ltd. to prepare the Illumina Ovine SNP50 Beadchip (The number of SNPs in this chip is 50 k). The validation samples were 130 PRS (first lambing) from Pishan Farm, with the same feeding conditions, 1.5–2.5 years old, no relationship.
2.3 Genotyping and data quality control
Genome Studio software was used to process the preliminary data results and obtain the VCF files. Plink () software was used for quality control. Unqualified SNP sites were eliminated. The quality control criteria of this study were as follows: 1) individual detection rate >0.95, 2) SNP detection rate >0.95, and 3) Hardy-Weinberg equilibrium (HWE) with p values .
2.4 LD calculation method
LD is the basis of association analysis, and the analysis of LD between loci helps to understand the LD level of the PRS genome. Since r2 is more capable of objectively reflecting the LD between different loci, it was adopted in this study as the LD measurement standard. The LD values ranged from 0 to 1. As the LD level increased with the value of r2, the linkage degree increased. The calculation formula of r2 is as follows ():Where PA1 and PB1 are the frequency of the first allele at the two marker loci and the haplotype frequency formed between alleles. The correlation coefficient (r2 mean) of alleles was calculated to measure the level of linkage disequilibrium (LD) using PopLDdecay V3.41 (), and perl scripts were used to visualize the results.
2.5 Estimation of effective population size
The SMC++ () method was used to estimate Ne. The population size history and splitting time of the PRS can be predicted with SMC++. A new spline regularization scheme was adopted in this method, significantly reducing estimation errors. The conversion of each VCF file into an input file in SMC++ format was made using the vcf2smc script distributed by SMC++. All the simulations were performed under the initial condition of a mutation rate of 1.25 × 10−8.
2.6 Genetic diversity and population structure
The genotypic data after quality control was subjected to Principal Component Analysis (PCA) using MingPCACluster (https://github.com/hewm2008/MingPCACluster). The VCF2Dis v1.09 (https://github.com/BGI-shenzhen/VCF2Dis) was used to calculate the P distance matrix, and then the NJ-tree was constructed by ATGC:FastME (http://www.atgc-montpellier.fr/fastme) program. Genetic admixture calculations were performed using Admixture ().
2.7 Fixation index
FST is used to measure the degree of population differentiation and can reflect the level of species population differentiation. This method is suitable for selective signal detection of multiple populations and as follows:Where, MSG is the mean square of error within the population, MSP is the mean square of error between the populations, and nc is the average sample size between the populations after correction. By using a sliding window with a window size of 50 Kb and a sliding step size of 25 Kb, the FST value of each sliding window SNP is calculated. Vcftools was used to calculate the FST value for each window, and then CMplot was used to plot Manhattan .
2.8 Integrated haplotype score
An intrapopulation selective genomic sweep analysis was performed on all individuals using iHS. iHS is an alternative EHH statistic using a single marker loci haploid type. It is defined as the core in the site, the expansion of the ancestors in the core loci alleles in the haploid type, and new mutant alleles in the extension of the haploid type EHH statistics for the integral genetic distance. It is possible to calculate the ratio between the previously mentioned genetic metrics to select a signal detection statistic using this method, as expressed by :uniHS is:Where IHH (Integrated EHH) refers to integrating genetic distance with EHH; A is the ancestral allele, and D is the newly derived allele.
2.9 Enrichment analysis of candidate genes
The iHS results and FST results were selected for intersection analysis, with annotations with the sheep genome Ovis Oar_v4.0. Gene functional annotation was performed referencing the NCBI databases (http://www.ncbi.nlm.nih.gov/gene) and OMIM database (http://www.ncbi.nlm.nih.gov/omim). The g:Profiler (https://biit.cs.ut.ee/gprofiler/gost) was used for autosomal enrichment of candidate genes for GO and Reactome/KEGG pathway analysis.
2.10 PCR amplification of BMPR1B
The FecB locus of the BMPR1B gene in 130 PRS was amplified by PCR according to the primers shown in Supplementary Material. After the PCR products were detected by 1.5% agarose gel electrophoresis, all qualified PCR products were sent to Beijing Compass Agritechnology Co., Ltd. for DNA sequencing and genotype identification.
3 Results
3.1 Descriptive statistics
Genotypic quality control was conducted on the SNPs of the 33 PRS and 40 QR used in the experiment. After the unqualified SNPs were removed, there were 49,219 informative SNPs in the PRS and QR population.
3.2 The extent of genome-wide LD and effective population size of PRS
When LD was calculated, the distance between markers was set in the 0–1 Mb (0–10, 10–25, 25–50, 50–100, 100–500 Kb, 0.5–1 Mb) autosomal range (Table 1). The average r2 decreased with increasing physical distance (Figure 1A). When the distance between markers was 10–25 Kb, the average r2 value was 0.151 ± 0.200. Compared with the average r2 value between SNPs within 0–10 Kb (0.233 ± 0.280), the difference was 0.082, much larger than the average r2 value between other adjacent distance regions. When the distance increased to 50–100 Kb, the r2 value was lower than 0.1, indicating that the LD was weak.
TABLE 1
| Distance | Average r2 | Number of SNP pairs | Proportion of r2 > 0.2 | Proportion of r2 > 0.3 |
|---|---|---|---|---|
| 0-10 Kb | 0.233 ± 0.280 | 8561 | 0.364 | 0.288 |
| 10-25 Kb | 0.151 ± 0.200 | 10,948 | 0.233 | 0.150 |
| 25-50 Kb | 0.120 ± 0.158 | 18,123 | 0.177 | 0.099 |
| 50-100 Kb | 0.096 ± 0.115 | 36,123 | 0.122 | 0.050 |
| 100-500 Kb | 0.078 ± 0.080 | 185,353 | 0.075 | 0.019 |
| 0.5‐1 Mb | 0.074 ± 0.073 | 368,857 | 0.062 | 0.014 |
LD statistical analysis between different distances (0–1 Mb).
r2: denotes the extent of LD.
FIGURE 1
Based on LD, six generations Ne of PRS were estimated (Figure 1B). The PRS remained relatively stable at about 631.97 during the 2,500–2,000 generations. It rose to 1,119.99 with 1,500 generations, 1,328.81 with 1,000 generations, and fell to 898.12 with 500 generations. It declined even faster until recent generations stood at 236.99. During the first 200-100 years, PRS populations remained relatively stable. It was only in the last 10 years that the population increased slightly due to the introduction of measures to protect the endemic species.
3.3 Genetic diversity and population structure
PCA analysis can separate PRS from QR, and part of PRS extended outward (Figure 2A). In Figure 2A, k-means clustering was carried out according to the predefined subsets. When k = 2, the distinction was obvious. Neighbor-Joining (N-J) Tree showed that PRS and QR were divided into two varieties (Figure 2B), which was consistent with PCA analysis. Admixture analysis showed that PRS and QR had similar genetic backgrounds (Figure 2C).
FIGURE 2
3.4 Selective gene sweep
The FST values of PRS and QR were calculated, and the FST values were arranged in descending order. The first 5% was regarded as the significant window (Figure 3A), and a total of 1148 genes were obtained. Under the 1% threshold, 184 candidate genes were screened out with iHS (Figure 3B). A total of 29 genes were obtained from the intersection of the two gene sets (Figure 3C), scuh as BMPR1B, 3BHSD, STPG2, ATRN, GANS, etc (Tabel A2. xlsx).
FIGURE 3
3.5 Functional annotation of genes
A total of 184 genes were screened by iHS and GO and Reactome pathway analysis were performed on the 184 genes. The enrichment results of the Reactome pathway of the candidate genes were shown in Figure 4. Some genes were related to animal reproduction (Table 2). Genes ABCA4, ARL4C, and PARP14 can affect the function of ion channels. Genes SOX2, DAB1, and COL5A2 regulate cell differentiation. 29 genes were obtained by the intersection of the results of FST and IHS, and after pathway enrichment, these genes were found to be related to GnRH signalingpathway, Ovarian steroidogenesis and Estrogen signaling pathway (Figure 5), indicateing that these genes may control PRS prolific trait.
FIGURE 4
TABLE 2
| Gene symbol | NCBI gene ID | Go term name | Go term ID | OAR | Coordinates (bp) |
|---|---|---|---|---|---|
| PAX5 | 101108719 | developmental process involved in reproduction | GO:0003006 | 2 | 51663045–51849640 |
| reproductive process | GO:0022414 | ||||
| sexual reproduction | GO:0019953 | ||||
| HERC2 | 101102534 | developmental process involved in reproduction | GO:0003006 | 2 | 113570525–113814058 |
| reproductive process | GO:0022414 | ||||
| Reproduction | GO:0000003 | ||||
| multicellular organism reproduction | GO:0032504 | ||||
| multi-organism reproductive process | GO:0044703 | ||||
| HSPA5 | 780447 | reproductive structure development | GO:0048608 | 3 | 10783159–10787445 |
| developmental process involved in reproduction | GO:0003006 | ||||
| reproductive process | GO:0022414 | ||||
| multicellular organismal reproductive process | GO:0048609 | ||||
| INSR | 443431 | reproductive structure development | GO:0048608 | 5 | 13788125–13936266 |
| developmental process involved in reproduction | GO:0003006 | ||||
| CAST | 443364 | reproductive process | GO:0022414 | 5 | 93695915–93785491 |
| sexual reproduction | GO:0019953 | ||||
| BMPR1B | 443454 | reproductive structure development | GO:0048608 | 6 | 30030664–30482585 |
| developmental process involved in reproduction | GO:0003006 | ||||
| reproductive process | GO:0022414 | ||||
| cellular process involved in reproduction in multicellular organism | GO:0022412 | ||||
| Reproduction | GO:0000003 | ||||
| reproductive system development | GO:0061458 | ||||
| multicellular organism reproduction | GO:0032504 | ||||
| multi-organism reproductive process | GO:0044703 | ||||
| GABRB1 | 101122358 | Reproduction | GO:0000003 | 6 | 66336302–66776044 |
| multicellular organism reproduction | GO:0032504 | ||||
| HAS2 | 101110341 | multicellular organism reproduction | GO:0032504 | 9 | 31108300–31140164 |
| TLR5 | 554256 | reproductive structure development | GO:0048608 | 12 | 25835725–25861758 |
| reproductive system development | GO:0061458 | ||||
| C11ORF74 | 101107893 | developmental process involved in reproduction | GO:0003006 | 15 | 65840095–65907755 |
| reproductive process | GO:0022414 | ||||
| Reproduction | GO:0000003 | ||||
| multi-organism reproductive process | GO:0044703 | ||||
| FGF10 | 443074 | reproductive structure development | GO:0048608 | 16 | 30604216–30700637 |
| developmental process involved in reproduction | GO:0003006 | ||||
| reproductive process | GO:0022414 | ||||
| reproduction | GO:0000003 | ||||
| ABHD2 | 101103014 | developmental process involved in reproduction | GO:0003006 | 18 | 19910266–20020987 |
| reproductive process | GO:0022414 | ||||
| cellular process involved in reproduction in multicellular organism | GO:0022412 | ||||
| Reproduction | GO:0000003 | ||||
| sexual reproduction | GO:0019953 | ||||
| multicellular organism reproduction | GO:0048609 | ||||
| multicellular organismal reproductive process | GO:0044703 | ||||
| multi-organism reproductive process | |||||
| ZFP42 | 101110027 | reproductive structure development | GO:0048608 | 26 | 16740249–16747149 |
| developmental process involved in reproduction | GO:0003006 | ||||
| reproductive process | GO:0022414 |
Candidate genes related to reproduction.
FIGURE 5
3.6 PCR product sequencing and genotyping
After 1.5% agarose gel electrophoresis, the PCR products were found to conform to the expected size, indicating that the target fragment was successfully amplified (Figure 6A). The sequencing results are shown in Figure 6B. SNP locus detection revealed that BMPR1B gene G.431965A > G, and three genotypes, B+, ++, and BB, were detected (Table 3). The number of lambs of the BB genotype was significantly higher than that of the B+ and ++ genotypes (p < 0.05). The number of lambs of the B+ genotype was significantly higher than that of the ++ genotype (p < 0.05). It was shown in the χ2 test that the PRS reached Hardy-Weinberg equilibrium at the FecB locus.
FIGURE 6
TABLE 3
| Genotype | Number of samples | Genotype frequency | Allele frequency | Average number of lambs | χ2 | |
|---|---|---|---|---|---|---|
| B | + | |||||
| B+ | 64 | 0.492 | 0.5077 | 0.4923 | 1.453 ± 0.063a | 1.1151 |
| BB | 34 | 0.262 | 1.735 ± 0.077b | |||
| ++ | 32 | 0.246 | 1.125 ± 0.059c | |||
The genotype frequency and litter size of different BMPR1B genotypes of Pishan Red sheep.
Note: Different letters of shoulder labels of data in the same column indicate significant difference (p < 0.05).
4 Discussion
The Pishan Red sheep is a newly discovered local sheep group known for its stress resistance, perennial estrus, and high fecundity. Understanding the breed characteristics and evolutionary history of PRS will help to further strengthen the conservation and utilization of its genetic resources. PCA, N-J tree and Admixture suggest that PRS and QR could be subdivided into two genetic clusters, and have similar ancestral components.
4.1 LD and Ne population size
In PRS, LD was only moderate at 0–10 Kb and rapidly decreased to 0.120 ± 0.158 at 25–50 Kb. The r2 we calculated for PRS was close to that of other sheep breeds and species. In Iranian Zandi sheep, the average r2 for the pairwise space of 0–10 Kb was 0.26 (). In the Barbaresca sheep, the average r2 for the intermarker distance of 0.5–1.0 Mb was 0.12 (). However, in Border Leicester and Poll Dorset, the average r2 for the intermarker distance of 0–10 Kb was 0.34 and 0.33, respectively (). The average r2 for the intermarker distances of 0–10 Kb and 10–25 Kb was 0.43 and 0.26, respectively, in Vrindavani cattle (). In Gir cattle, the average r2 value of 0.5–1 Mb was 0.032 (). Thus, the variation in the reported r2 in different breeds and species suggests that LD are highly specific in sheep breeds. The PRS is located at the edge of the Taklimakan Desert. Its unique geographical environment may be the main reason for its high genetic diversity.
Mean r2 values varied on chromosomes (from 0.182 ± 0.269 in OAR25 to 0.285 ± 0.367 in OAR19, distance <10 Kb), consistent with previous reports in sheep (), beef cattle (), and dairy cows (). This phenomenon may be due to differences in recombination rates in different chromosomes, natural or artificial selection, and genetic drift (). In the Vrindavani cattle, the highest mean r2 value of chromosome 28 was 0.643, and the lowest r2 value of chromosome 18 was 0.172 (0–10 Kb) (). Moreover, the variation in r2 estimated for different chromosomes was higher in short SNP pair distances, which is in line with the results reported for Chinese Merino sheep (Xinjiang type) (). As seen from the attenuation of LD of 26 chromosomes of PRS, the LD of each chromosome was weak. The r2 values were higher where the distance between the marked sites was close, but there were also higher r2 values between the two sites, indicating a certain pattern of LD between the distant sites. The weak LD levels of OAR12, OAR18, OAR20, and OAR25 indicated that the degree of purification of these chromosomes was not high, leading to higher genetic diversity than other chromosomes. The chromosomes with higher r2 values were OAR19, OAR15, OAR13, and OAR14, indicating that these four chromosomes may be more strongly selected than the other chromosomes.
The Ne of 236.99 of PRS 50 generations ago was similar to that of Chinese Merino sheep (). We observed that Ne decreased more strongly from about 380 generations ago, consistent with the results of Chinese Merino sheep (). The low level of r2, even at relatively short distances, showed that the Ne in PRS was large in recent past generations compared with other species. For example, the r2 for SNP pairs within 0.9–1.0 Mb and the Ne in recent generations in Duroc pigs were reported to be 0.2 and 75, respectively (). However, given the sharp drop in Ne in recent generations, we should be careful to maintain Ne larger than 100 individuals.
4.2 Adaptive mechanisms of the desert environment
The environmental adaptability differs between the Taklimakan Desert sheep breeds and those of other areas. PRS can adapt to extreme conditions, such as high salinity, drought, and ultraviolet rays. Under the iHS 1% threshold, 184 genes were screened. It was revealed that the genetic evidence and physiological mechanism of PRS adaptation to the desert environment. In terms of linoleic acid (LA) metabolism (R-MMU-2046105), FADS1 and FADS2 control polyunsaturated fatty acids (; ; ; ), which accumulate fat to cope with extreme weather. The activation of phospholipase C enzymes results in the generation of second messengers of the phosphatidylinositol pathway in terms of adaptability. The events resulting from this pathway increase intracellular calcium and protein kinase C (PKC) activation. Phospholipase C cleaves the phosphodiester bond in PIP2 to form 1,2 diacylglycerol (DAG) and 1,4,5-inositol trisphosphate (IP3). IP3 opens Ca2+ channels in the platelet-dense tubular system, raising intracellular Ca2+ levels (R-MMU-76005). In terms of immunity, HERC2 and USP10 can promote the activation of the ATR-CHK1 pathway, triggering cell cycle checkpoints (). PAALD is involved in the process of phagocytosis (). KIF2A, TRIM66, BCL11B, and CLEC14A play essential roles in repairing damaged cell DNA (; ; ; ). GPRIN3 and HERC3 control cell senescence and apoptosis (; ). In terms of growth and development, ZFP42 can control the differentiation of embryonic stem cells (). EVC can promote chondrogenesis (). GABRB1 correlates with hypothalamic volume and regulates intelligence ().
4.3 Genetic mechanisms of perennial estrus and reproduction
The estrus cycle refers to the time between the previous and next ovulation. The ovary undergoes follicular growth, maturation, ovulation, luteal formation, and degeneration during the estrus cycle. The vast majority of sheep are singleton and seasonally estrus, leading to the failure of a balanced supply of lamb meat in the four seasons, seriously restricting the production efficiency of the meat sheep industry. PRS live in desert environments, and after natural and artificial selection, it forms the characteristics of perennial estrus and the early onset of puberty. TLR5 can effectively alleviate the stimulation caused by radiation (), and SOX10 independently regulates the expression of IRF1 in melanoma through the JAK-STAT signaling pathway (). ATP6V0A can affect fetal brain development (). AUH can cause the early onset of puberty (). During the luteal phase of the estrus cycle, IFNE is highly expressed in the uterus and has a protective effect against uterine infection (). As a member of the BMP system, BMPR1B plays a significant role in the sheep ovary. The BMP system can control the proliferation and differentiation of ovarian granulosa cells and the development of oocytes, among which BMPR1B plays an essential role in the regulation of ovarian function. The SNP of BMPR1B C.746 A > G, the 249th amino acid change, partially inactivates BMPR1B protein. This change affects the reaction of the ligands GDF5 and BMP4 recognized by BMPR1B to steroid production, making follicles mature earlier and increasing the ovine ovulation number (). BMPR1B is the primary gene affecting the trait of lambing abundance in sheep (; ).
Fibroblast growth factor (FGF) is a large family of paracrine cells that can regulate follicular development and oocyte maturation. FGF10 can interact with BMP15 to increase cumulus cell diffusion and improve glucose utilization (). FGF10 was expressed in oocytes and membrane cells in cattle follicles, acting on granulosa cells to inhibit steroid production. FGF10 regulates the cumulus-oocyte complex, improving expansion and development ability (; ). FGF10 can also reduce the proportion of apoptotic oocytes and increase the number of cells developing to the blastocyst stage (). In embryonic development, FGF10 can activate the MAPK pathway, increase the phosphorylation level of MAPK, and mediate the migration of sheep trophoblast cells (). FGF10 can reduce the level of FSHR mRNA in granulosa cells in vitro. FGF10 can also inhibit the secretion of estrogen. The injection of FGF10 at the initial stage of follicular selection can inhibit follicular growth, which may be achieved by reducing FSHR mRNA levels and thus inhibiting estrogen expression. At the same time, the FGF10 mRNA concentration in the ovaries of healthy and growing bovine follicles was higher than that in atresia follicles (), which indicated that FGF10 played an inhibitory role in early follicular development and promoted follicular maturation in late follicular development. FGF10 can also regulate the expression of genes CD9, CD81, DNMT1, and DNMT3B to improve embryo quality ().
Hyaluronic acid (HA), also called hyaluronate or hyaluronan, is an essential extracellular matrix component widely available in various mammalian tissues (). HA is synthesized by three HA synthase families (HAS1, HAS2, and HAS3). HAS2 is so critical to development that the HAS2-HA system is considered an essential HAS-HA system. HAS2 is responsible for the rapid hyaluronic acid synthesis in cumulus-oocyte complexes and granulosa cells. In the dominant follicles of mammals, most of the hyaluronic acid is secreted by cumulus cells and is also present in the follicular fluid. HAS2 plays a vital physiological role in oocyte maturation, ovulation, in vivo fertilization, and early embryonic development (). C11orf74 enhances sperm motility and improves fertilization (). Progesterone can upregulate ABHD2 to activate the camp-PKA signaling pathway (). These genes were strongly selected, suggesting that perennial estrus may result from a combination of these genes.
The number of lambs is a complex character affected by many factors, among which heredity is the main factor. The FecB gene is the most studied gene and can significantly affect the ovulation number and lambing number of sheep breeds. The B+ genotype exists among multiple sheep breeds in Xinjiang, China, including Duolang and Hotan sheep. The number of lambs of the B+ genotype was significantly higher than that of the ++ genotype. We detected FecB in 130 leather-red sheep and found three genotypes: BB, B+, and ++. The average number of lambs born in the BB genotype was the highest, followed by the B+ and ++ genotypes, indicating that the FecB gene is an effective gene for improving the fecality of PRS.
5 Conclusion
In the extreme environment of the desert, PRS has formed the characteristics of perennial estrus, multiple pregnancies, and adequate resistance to stress. The LD decay rate of the PRS was faster, and LD and the Ne were at a relatively low level, suggesting that we should take reasonable protection measures to increase the PRS population. Genomic selection signal analysis and population validation showed that FecB could be used as a molecular breeding marker for multiple fetal lines of PRS in desert environments.
Statements
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: Our data has been uploaded to this website https://figshare.com/articles/dataset/e_vcf/21700937.
Author contributions
C-lZ and JZ led the bioinformatic and statistical analyses of data and helped to draft the first version of the manuscript. SL generated and contributed 50 K SNP data for sheep breeds in the Xinjiang biota. JZ, SL, MT, and C-lZ wrote and/or revised the manuscript. All authors read and approved the manuscript.
Funding
This study was funded by grants form National Natural Science Foundation of China (32060743) supported by Bintuan Science and Technology Program (2022CB001-09).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2023.1092066/full#supplementary-material
References
1
AkeyJ. M.ZhangG.ZhangK.JinL.ShriverM. D. (2002). Interrogating a high-density SNP map for signatures of natural selection. Genome Res.12, 1805–1814. 10.1101/gr.631202
2
Al-MamunH. A.ClarkS. A.KwanP.GondroC. (2015). Genome-wide linkage disequilibrium and genetic diversity in five populations of Australian domestic sheep. Genet. Sel. Evol.47, 90. 10.1186/s12711-015-0169-6
3
Al-SamerriaS.Al-AliI.McFarlaneJ. R.AlmahbobiG. (2015). The impact of passive immunisation against BMPRIB and BMP4 on follicle development and ovulation in mice. Reproduction149, 403–411. 10.1530/rep-14-0451
4
AotoK.KatoM.AkitaT.NakashimaM.MutohH.AkasakaN.et al (2021). ATP6v0a1 encoding the a1-subunit of the v0 domain of vacuolar h+-ATPases is essential for brain development in humans and mice. Nat. Commun.12, 2107. 10.1038/s41467-021-22389-5
5
BizjakN.TansekM. Z.StefanijaM. A.LampretB. R.MezekA.TorkarA. D.et al (2020). Precocious puberty in a girl with 3-methylglutaconic aciduria type 1 (3-MGA-i) due to a novel AUH gene mutation. Mol. Genet. Metabolism Rep.25, 100691. 10.1016/j.ymgmr.2020.100691
6
BrackettC. M.GreeneK. F.AldrichA. R.TrageserN. H.PalS.MolodtsovI.et al (2021). Signaling through TLR5 mitigates lethal radiation damage by neutrophil-dependent release of MMP-9. Cell. Death Discov.7, 266. 10.1038/s41420-021-00642-6
7
BuratiniJ.PintoM.CastilhoA.AmorimR.GiomettiI.PortelaV.et al (2007). Expression and function of fibroblast growth factor 10 and its receptor, fibroblast growth factor receptor 2B, in bovine follicles. Biol. Reproduction77, 743–750. 10.1095/biolreprod.107.062273
8
CaixetaE. S.Sutton-McDowallM. L.GilchristR. B.ThompsonJ. G.PriceC. A.MachadoM. F.et al (2013). Bone morphogenetic protein 15 and fibroblast growth factor 10 enhance cumulus expansion, glucose uptake, and expression of genes in the ovulatory cascade during in vitro maturation of bovine cumulus–oocyte complexes. Reproduction146, 27–35. 10.1530/rep-13-0079
9
CalderP. C. (2015). Functional roles of fatty acids and their effects on human health. J. Parenter. Enter. Nutr.39, 18S–32S. 10.1177/0148607115595980
10
ChenJ.WangZ.GuoX.LiF.WeiQ.ChenX.et al (2019). TRIM66 reads unmodified h3r2k4 and h3k56ac to respond to DNA damage in embryonic stem cells. Nat. Commun.10, 4273. 10.1038/s41467-019-12126-4
11
ChenY.LiY.PengY.ZhengX.FanS.YiY.et al (2018). δnp63α down-regulates c-myc modulator mm1 via e3 ligase herc3 in the regulation of cell senescence. Cell. Death Differ.25, 2118–2129. 10.1038/s41418-018-0132-5
12
CorbinL. J.BlottS. C.SwinburneJ. E.VaudinM.BishopS. C.WoolliamsJ. A. (2010). Linkage disequilibrium and historical effective population size in the thoroughbred horse. Anim. Genet.41, 8–15. 10.1111/j.1365-2052.2010.02092.x
13
DasU. (2006). “Essential fatty acids,” in Encyclopedia of biophysics (Heidelberg: Springer Berlin Heidelberg), 706–714. 10.1007/978-3-642-16712-6_533
14
DingF.LuoX.TuY.DuanX.LiuJ.JiaL.et al (2021). Alpk1 sensitizes pancreatic beta cells to cytokine-induced apoptosis via upregulating TNF-α signaling pathway. Front. Immunol.12, 705751. 10.3389/fimmu.2021.705751
15
EdeaZ.DadiH.KimS.-W.ParkJ.-H.ShinG.-H.DessieT.et al (2014). Linkage disequilibrium and genomic scan to detect selective loci in cattle populations adapted to different ecological conditions in Ethiopia. J. Animal Breed. Genet.131, 358–366. 10.1111/jbg.12083
16
FischerC. D.Wachoski-DarkG. L.GrantD. M.BramerS. A.KleinC. (2018). Interferon epsilon is constitutively expressed in equine endometrium and up-regulated during the luteal phase. Animal Reproduction Sci.195, 38–43. 10.1016/j.anireprosci.2018.05.003
17
GhoreishifarS. M.Moradi-ShahrbabakH.ParnaN.DavoudiP.KhansefidM. (2019). Linkage disequilibrium and within-breed genetic diversity in iranian zandi sheep. Arch. Anim. Breed.62, 143–151. 10.5194/aab-62-143-2019
18
GrossiD. A.JafarikiaM.BritoL. F.BuzanskasM. E.SargolzaeiM.SchenkelF. S. (2017). Genetic diversity, extent of linkage disequilibrium and persistence of gametic phase in canadian pigs. BMC Genet.18, 6. 10.1186/s12863-017-0473-y
19
HansonL. Å.KorotkovaM. (2002). The role of breastfeeding in prevention of neonatal infection. Seminars Neonatol.7, 275–281. 10.1016/s1084-2756(02)90124-7
20
HayesB. J.LewinH. A.GoddardM. E. (2013). The future of livestock breeding: Genomic selection for efficiency, reduced emissions intensity, and adaptation. Trends Genet.29, 206–214. 10.1016/j.tig.2012.11.009
21
HillW. G. (1974). Estimation of linkage disequilibrium in randomly mating populations. Heredity33, 229–239. 10.1038/hdy.1974.89
22
Irving-RodgersH.RodgersR. (2006). Extracellular matrix of the developing ovarian follicle. Seminars Reproductive Med.24, 195–203. 10.1055/s-2006-948549
23
JiangF.ZhuY.ChenY.TangX.LiuL.ChenG.et al (2021). Progesterone activates the cyclic AMP-protein kinase a signalling pathway by upregulating ABHD2 in fertile men. J. Int. Med. Res.49, 300060521999527. 10.1177/0300060521999527
24
KijasJ. W.Porto-NetoL.DominikS.ReverterA.BunchR.McCullochR.et al (2014). Linkage disequilibrium over short physical distances measured in sheep using a high-density SNP chip. Anim. Genet.45, 754–757. 10.1111/age.12197
25
KumarS.RajputP. K.BahireS. V.JyotsanaB.KumarV.KumarD. (2020). Differential expression of BMP/SMAD signaling and ovarian-associated genes in the granulosa cells of FecB introgressed GMM sheep. Syst. Biol. Reproductive Med.66, 185–201. 10.1080/19396368.2019.1695977
26
LamuedraA.GratalP.CalatravaL.Ruiz-PerezV. L.Palencia-CamposA.Portal-NúñezS.et al (2022). Blocking chondrocyte hypertrophy in conditional evc knockout mice does not modify cartilage damage in osteoarthritis. FASEB J.36, e22258. 10.1096/fj.202101791rr
27
LawrenceR. M.LawrenceR. A. (2004). Breast milk and infection. Clin. Perinatology31, 501–528. 10.1016/j.clp.2004.03.019
28
LiuS.HeS.ChenL.LiW.DiJ.LiuM. (2017). Estimates of linkage disequilibrium and effective population sizes in Chinese merino (xinjiang type) sheep by genome-wide SNPs. Genes. and Genomics39, 733–745. 10.1007/s13258-017-0539-2
29
LvX.ChenL.HeS.LiuC.HanB.LiuZ.et al (2020). Effect of nutritional restriction on the hair follicles development and skin transcriptome of Chinese merino sheep. Animals10, 1058. 10.3390/ani10061058
30
MajkowskiM.LaszkiewiczA.SniezewskiL.GrzmilP.PawlickaB.TomczykI.et al (2018). Lack of NWC protein (c11orf74 homolog) in murine spermatogenesis results in reduced sperm competitiveness and impaired ability to fertilize egg cells in vitro. PLOS ONE13, e0208649. 10.1371/journal.pone.0208649
31
MastrangeloS.PortolanoB.GerlandoR. D.CiampoliniR.ToloneM.SardinaM.et al (2017). Genome-wide analysis in endangered populations: A case study in barbaresca sheep. Animal11, 1107–1116. 10.1017/s1751731116002780
32
MulsantP.LecerfF.FabreS.SchiblerL.MongetP.LannelucI.et al (2001). Mutation in bone morphogenetic protein receptor-IB is associated with increased ovulation rate in booroola mérino ewes. Proc. Natl. Acad. Sci.98, 5104–5109. 10.1073/pnas.091577598
33
OspinaA. M. T.MaioranoA. M.CuriR. A.PereiraG. L.Zerlotti-MercadanteM. E.CyrilloJ. N. S. G.et al (2019). Linkage disequilibrium and effective population size in gir cattle selected for yearling weight. Reproduction Domest. Animals54, 1524–1531. 10.1111/rda.13559
34
PanY.WangM.BalochA. R.ZhangQ.WangJ.MaR.et al (2019). FGF10 enhances yak oocyte fertilization competence and subsequent blastocyst quality and regulates the levels of CD9, CD81, DNMT1, and DNMT3b. J. Cell. Physiology234, 17677–17689. 10.1002/jcp.28394
35
PattersonN.MoorjaniP.LuoY.MallickS.RohlandN.ZhanY.et al (2012). Ancient admixture in human history. Genetics192, 1065–1093. 10.1534/genetics.112.145037
36
PintoR. P.FontesP.LoureiroB.CastilhoA. S.TicianelliJ. S.RazzaE. M.et al (2014). Effects of FGF10 on bovine oocyte meiosis progression, apoptosis, embryo development and relative abundance of developmentally important genes in vitro. Reproduction Domest. Animals50, 84–90. 10.1111/rda.12452
37
Porto-NetoL. R.KijasJ. W.ReverterA. (2014). The extent of linkage disequilibrium in beef cattle breeds using high-density SNP genotypes. Genet. Sel. Evol.46, 22. 10.1186/1297-9686-46-22
38
PurcellS.NealeB.Todd-BrownK.ThomasL.FerreiraM. A.BenderD.et al (2007). Plink: A tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet.81, 559–575. 10.1086/519795
39
QanbariS.PimentelE. C. G.TetensJ.ThallerG.LichtnerP.SharifiA. R.et al (2009). The pattern of linkage disequilibrium in German holstein cattle. Anim. Genet.41, 346–356. 10.1111/j.1365-2052.2009.02011.x
40
RaoY. S.LiangY.XiaM. N.ShenX.DuY. J.LuoC. G.et al (2008). Extent of linkage disequilibrium in wild and domestic chicken populations. Hereditas145, 251–257. 10.1111/j.1601-5223.2008.02043.x
41
RodgersR. J.Irving-RodgersH. F. (2010). Formation of the ovarian follicular antrum and follicular fluid. Biol. Reproduction82, 1021–1029. 10.1095/biolreprod.109.082941
42
ScotlandK. B.ChenS.SylvesterR.GudasL. J. (2009). Analysis of rex1 (zfp42) function in embryonic stem cell differentiation. Dev. Dyn.238, 1863–1877. 10.1002/dvdy.22037
43
SeiraO.LiuJ.AssinckP.RamerM.TetzlaffW. (2019). KIF2a characterization after spinal cord injury. Cell. Mol. Life Sci.76, 4355–4368. 10.1007/s00018-019-03116-2
44
SinghA.KumarA.MehrotraA.KarthikeyanA.PandeyA. K.MishraB. P.et al (2021). Estimation of linkage disequilibrium levels and allele frequency distribution in crossbred vrindavani cattle using 50k SNP data. PLOS ONE16, e0259572. 10.1371/journal.pone.0259572
45
SuZ.LiY.LvH.CuiX.LiuM.WangZ.et al (2021). CLEC14a protects against podocyte injury in mice with adriamycin nephropathy. FASEB J.35, e21711. 10.1096/fj.202100283r
46
SunH.-M.ChenX.-L.ChenX.-J.LiuJ.MaL.WuH.-Y.et al (2017). PALLD regulates phagocytosis by enabling timely actin polymerization and depolymerization. J. Immunol.199, 1817–1826. 10.4049/jimmunol.1602018
47
TenesaA.NavarroP.HayesB. J.DuffyD. L.ClarkeG. M.GoddardM. E.et al (2007). Recent human effective population size estimated from linkage disequilibrium. Genome Res.17, 520–526. 10.1101/gr.6023607
48
TerhorstJ.KammJ. A.SongY. S. (2016). Robust and scalable inference of population history from hundreds of unphased whole genomes. Nat. Genet.49, 303–309. 10.1038/ng.3748
49
UimariP.TapioM. (2011). Extent of linkage disequilibrium and effective population size in Finnish Landrace and Finnish Yorkshire pig breeds. J. Animal Sci.89, 609–614. 10.2527/jas.2010-3249
50
ValisnoJ. A. C.MayJ.SinghK.HelmE. Y.VenegasL.BudbazarE.et al (2021). BCL11b regulates arterial stiffness and related target organ damage. Circulation Res.128, 755–768. 10.1161/circresaha.120.316666
51
VoightB. F.KudaravalliS.WenX.PritchardJ. K. (2006). A map of recent positive selection in the human genome. PLoS Biol.4, e72. 10.1371/journal.pbio.0040072
52
WangJ. (2005). Estimation of effective population sizes from data on genetic markers. Philosophical Trans. R. Soc. B Biol. Sci.360, 1395–1409. 10.1098/rstb.2005.1682
53
YangQ. E.GiassettiM. I.EalyA. D. (2011). Fibroblast growth factors activate mitogen-activated protein kinase pathways to promote migration in ovine trophoblast cells. Reproduction141, 707–714. 10.1530/rep-10-0541
54
YinL.ZhangH.TangZ.XuJ.YinD.ZhangZ.et al (2021). rMVP: A memory-efficient, visualization-enhanced, and parallel-accelerated tool for genome-wide association study. Genomics, Proteomics Bioinforma.19, 619–628. 10.1016/j.gpb.2020.10.007
55
YokoyamaS.TakahashiA.KikuchiR.NishibuS.LoJ. A.HejnaM.et al (2021). SOX10 regulates melanoma immunogenicity through an IRF4–IRF1 axis. Cancer Res.81, 6131–6141. 10.1158/0008-5472.can-21-2078
56
YuanJ.LuoK.DengM.LiY.YinP.GaoB.et al (2014). HERC2-USP20 axis regulates DNA damage checkpoint through claspin. Nucleic Acids Res.42, 13110–13121. 10.1093/nar/gku1034
57
ZhangC.DongS.-S.XuJ.-Y.HeW.-M.YangT.-L. (2018). PopLDdecay: A fast and effective tool for linkage disequilibrium decay analysis based on variant call format files. Bioinformatics35, 1786–1788. 10.1093/bioinformatics/bty875
58
ZhangK.HansenP. J.EalyA. D. (2010). Fibroblast growth factor 10 enhances bovine oocyte maturation and developmental competence in vitro. Reproduction140, 815–826. 10.1530/rep-10-0190
59
ZhangX.ZhangL.SunW.LangX.WuJ.ZhuC.et al (2020). Study on the correlation between BMPR1b protein in sheep blood and reproductive performance. J. Animal Sci.98, skaa100. 10.1093/jas/skaa100
60
ZhuB.ChenC.XueG.LeiX.LiJ.MoyzisR. K.et al (2014). The GABRB1 gene is associated with thalamus volume and modulates the association between thalamus volume and intelligence. NeuroImage102, 756–763. 10.1016/j.neuroimage.2014.08.048
Summary
Keywords
desert environment, genomic selection, linkage disequilibrium, perennial estrus, litter size
Citation
Zhang C, Zhang J, Tuersuntuoheti M, Chang Q and Liu S (2023) Population structure, genetic diversity and prolificacy in pishan red sheep under an extreme desert environment. Front. Genet. 14:1092066. doi: 10.3389/fgene.2023.1092066
Received
07 November 2022
Accepted
28 March 2023
Published
11 April 2023
Volume
14 - 2023
Edited by
Ibrar Muhammad Khan, Fuyang Normal University, China
Reviewed by
Hongyu Liu, Anhui Agricultural University, China
Samiullah Khan, Guizhou University, China
Updates
Copyright
© 2023 Zhang, Zhang, Tuersuntuoheti, Chang and Liu.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Shudong Liu, liushudong63@126.com
† These authors have contributed equally to this work
This article was submitted to Livestock Genomics, a section of the journal Frontiers in Genetics
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.