Multiomics Analyses Provide New Insight into Genetic Variation of Reproductive Adaptability in Tibetan Sheep

Abstract

Domestication and artificial selection during production-oriented breeding have greatly shaped the level of genomic variability in sheep. However, the genetic variation associated with increased reproduction remains elusive. Here, two groups of samples from consecutively monotocous and polytocous sheep were collected for genome-wide association, transcriptomic, proteomic, and metabolomic analyses to explore the genetic variation in fecundity in Tibetan sheep. Genome-wide association study revealed strong associations between BMPR1B (p.Q249R) and litter size, as well as between PAPPA and lambing interval; these findings were validated in 1,130 individuals. Furthermore, we constructed the first single-cell atlas of Tibetan sheep ovary tissues and identified a specific mural granulosa cell subtype with PAPPA-specific expression and differential expression of BMPR1B between the two groups. Bulk RNA-seq indicated that BMPR1B and PAPPA expressions were similar between the two groups of sheep. 3D protein structure prediction and coimmunoprecipitation analysis indicated that mutation and mutually exclusive exons of BMPR1B are the main mechanisms for prolific Tibetan sheep. We propose that PAPPA is a key gene for stimulating ovarian follicular growth and development, and steroidogenesis. Our work reveals the genetic variation in reproductive performance in Tibetan sheep, providing insights and valuable genetic resources for the discovery of genes and regulatory mechanisms that improve reproductive success.

Keywords: GWAS; Tibetan sheep; adaptation; multiomics; reproduction.

Fig. 1.

Whole-genome sequencing found that BMPR1B and PAPPA were associated with reproductive phenotypes in Tibetan sheep. a) Manhattan plot showing the significant association between the mutation in BMPR1B with litter size differences between HF and LF groups. The significance threshold P–value [−log₁₀(P) = 8.87] is denoted by the line. b) Statistics of litter size in the two groups of Tibetan sheep. HF and LF indicate high-fertility and low-fertility groups. The bars represent mean ± SD. A two-tailed t test was used (***P < 0.001). c) Manhattan plot showing the significant association between a SNP in PAPPA on chromosome 2 with the lambing interval phenotype. The red horizontal line corresponds to the genome-wide significance threshold [−log₁₀(P) = 8.87]. d) Statistics of lambing interval in the two groups of Tibetan sheep (*P < 0.05; **P < 0.01; ***P < 0.001). e) Statistics on the proportion of BMPR1B genotypes in a large population with different reproductive phenotypes. The orange and blue histograms represent sheep with high and low fertility in the large population (n = 1,130). f) Statistics on the proportion of PAPPA genotypes in a large population (n = 1,130).

Fig. 2.

AS events and Co-IP analysis revealed that BMPR1B was associated with high production in Tibetan sheep. a) Schematic representation of Tibetan sheep ovarian tissue and follicular fluid preparation for multiomics analysis. b) Statistics of dominant follicle/total follicle numbers in HF and LF groups of Tibetan sheep (**P < 0.01). c) Statistics of the dominant follicle volume in the two groups of Tibetan sheep (***P < 0.001). d) GSEA results of the KEGG pathway “Ovarian steroidogenesis” is shown. e) Differences in MXE between the HF and LF groups. f) The 3D structure, predicted by AlphaFold2, of proteins B1, B2, and B3 encoded by variable splicing, with cysteine residues in the extracellular domain, the GS motif, transmembrane domain, loop 45, and kinase domain. Protein B3 included a mutation from Q to R on exon 8 of protein B1. g) Co-IP assays show that FKBP1A binds to B1 in Tibetan sheep. h) Co-IP assays show that sheep BMP4 complexes with B1, B2, and B3.

Fig. 3.

Clustering analysis of various cells in the ovaries of Tibetan sheep. a) UMAP cluster map revealing 21 (CL0 to CL20) specific clusters representing the major ovarian cell types. b) UMAP cluster map showing the seven major ovarian cell types in Tibetan sheep. c) Dot plot showing distinct expression patterns of the selected features genes for each cell type. The expression level of each gene is indicated by a color gradient from low to high. The percentage of cells expressing a specific gene is indicated by the size of a dot. d) Heatmap and hierarchical clustering based on the expression of the top 50 most variable genes.

Fig. 4.

Analysis of DEGs in the ovaries of Tibetan sheep. a) Up- and downregulated genes across all clusters. The x axis represents CL0 to CL20, and y axis indicates the average log₂FC. b) GSEA results of the KEGG pathway “TGF-β signaling pathway” is shown. c) Expression heatmap showing DEGs between HF and LF groups. d) UMAP cluster map showing expression of genes characteristic of the major ovarian GC types. Blue dashed lines indicate the boundaries of the main clusters of interest.

Fig. 5.

Pseudotime and clustering analyses of mural GCs and oocytes. a) Pseudotime trajectory analysis of mural GC (CL14). b) Heatmap pseudotime-dependent marker genes for mural GC subtypes. The expression level of each gene is indicated by a color gradient from low to high . c) KEGG enrichment of DEGs in ovarian mural GC subtypes of sheep with different lambing numbers. d) Expression trends of signature genes in mural GC subtypes arranged along pseudotime. e) Analysis of cell trajectories of oocytes (CL8) by Monocle. f) Pseudotime-ordered heatmap demonstrating the pseudotime order of selected marker genes for oocyte subtypes. The expression level of each gene is indicated from low to high by a color gradient from blue to red, respectively. g) KEGG enrichment of DEGs in ovarian oocyte subtypes of sheep with different lambing numbers. h) Expression plots demonstrating expression trends of the signature genes in oocyte subtypes arranged along the pseudotime.

Fig. 6.

Proteomic and metabolomic analyses of ovaries and follicular fluid. a) A volcano plot shows the results of a Student's t test for DEPs (P < 0.05, fold change > 1) between the LF and HF groups. b) Sankey plot showing DMs were involved in each of the enriched pathways obtained via KEGG. Dot plot shows the ratio between DMs and the total number of metabolites in each enriched pathway (false detection rate–adjusted P ≤ 0.05). c to f) Determination of serum hormone levels in both groups by ELISA. E₂, estradiol; P₄, progesterone; T, testosterone; LH, luteinizing hormone. The bars display mean ± SD. A two-tailed t test was used (**P < 0.01; ***P < 0.001).

Fig. 7.

The main pathways of reproductive candidate genes for follicular growth and ovulation in Tibetan sheep. Upregulated and downregulated reproductive candidate genes during follicle growth and ovulation are highlighted in red and blue, respectively.