Integration analysis of pituitary proteome and transcriptome reveals fertility–related biomarkers in FecB mutant Small Tail Han sheep

The Booroola fecundity mutation (FecB) in Small Tail Han sheep has been shown to enhance ovulation rates and litter sizes by affecting the hypothalamic–pituitary–gonadal (HPG) axis. Despite the pituitary’s role in reproductive regulation, its involvement in FecB-induced ovulation remains understudied. Our study aimed to fill this gap by analyzing pituitary tissues from FecB homozygous (BB) and wild-type (WW) ewes during luteal and follicular phases using tandem mass tag–based protein quantification and the DIABLO framework for proteomic and transcriptomic data integration. Significant differences in 277 proteins were observed across estrus periods, with network analysis highlighting the voltage-dependent calcium channel L-type alpha-1C as a key convergence point in oxytocin signaling and GnRH secretion pathways. The DIABLO method revealed a strong correlation (0.98) between proteomic and transcriptomic datasets, indicating a coordinated response in FecB ewes. Notably, higher expression levels of Follicle Stimulating Hormone Subunit Beta (FSHB) and Luteinizing Hormone Subunit Beta (LHB) were found in BB ewes during the follicular phase, potentially due to elevated E2 concentrations. Furthermore, our analysis identified genes related to the Gamma–aminobutyric acid type A receptor family (GABRA2, GABRG1, GABRB1) in the pituitary, with GABRB1 showing higher expression in BB ewes. This suggests a role for GABA in modulating GnRH and gonadotropin feedback loops, potentially contributing to the FecB mutation’s effect on ovulation. This study provides novel insights into the pituitary’s role in fertility among FecB sheep, identifying GABA as a potential regulatory factor within the HPG axis. The findings also open avenues for discovering new biomarkers in pituitary endocrinology for sheep breeding purposes.

1 Introduction

The Bone morphogenetic mechanism type 1 B (BMPR1B) gene was the first major gene reported to increase ovulation rate, which found in the Booroola prolific sheep strain. BMPR1B^Q249R (FecB) gene mutation significantly increases the number of ovulations per estrus cycle in ewes (1–3). Recently, Zhou et al. introduced the FecB mutation into low–fertility sheep using gene editing technology and observed that this mutation is a crucial causal mutation influencing high sheep productivity (4, 5). Our laboratory results have found that Small Tail Han sheep (STH) carrying the FecB mutation exhibit earlier estrus onset and ovulation time compared to wild–type (WW) ewes. Additionally, these STHs have smaller mature follicle diameters but significantly higher rates of ovulation and litter sizes compared to their WW counterparts (6, 7).

Both in vivo pituitary cell and sheep studies at the in vitro level showed that follicle stimulating hormone levels were significantly higher in BB ewes than in WW ewes (P < 0.05) (8, 9). In STH sheep, at the 3rd hour after the first estrus, follicle stimulating hormone (FSH) concentrations were significantly higher in BB sheep than in WW sheep (P < 0.05), while luteinising hormone (LH) concentrations were significantly lower than in WW sheep (P < 0.05) (6). The FecB mutation have been observed to potentially lead to a partial attenuation of the transforming growth factor beta (TGFβ) pathway, consequently impacting the follicular response to FSH and LH. Furthermore, this mutation may change the hormonal feedback loop between the follicle and the pituitary gland, thereby influencing the release of FSH and LH (10).

In mammals that ovulate spontaneous like sheep, gonadotropin (FSH and LH) secretion by the pituitary gland drives follicle formation and a massive discharge or surge of gonadotropin that triggers ovulation during a typical ovarian cycle (11). The pituitary, the “master gland of the body”, secretes hormones that are in essence peptides or glycoproteins. Thus, transcriptomics and quantitative proteomics have expanded our understanding of the pituitary endocrine and pituitary regulation of reproduction (12–18). However, these studies have mainly focused on pituitary regulation of estrus, and the molecular mechanisms of pituitary involvement in follicular development to regulate litter size remain poorly understood with the FecB mutations. Integration of systems biology and innovative data can unravel complex biology with new insights (19).

In this experiment, we collected pituitary tissues from luteal and follicular phases of STH ewes with different FecB genotypes, and detected the changes in protein abundance in these tissues by using tandem mass tagging (TMT) quantitative proteomics and parallel reaction monitoring (PRM) methods. Previously, our group obtained transcriptomic data from these tissues. We performed a combination of these two histological data to screen for core biomarkers, signaling pathways, and biological processes associated with FecB gene mutations, in the pituitary gland, that are involved in ovulation.

2 Materials and methods

2.1 Study design and tissues collection

TaqMan MGB method was used to detect FecB mutant genotypes FecB (6) in 2–4 years old, similar weight and non–pregnant ewes in the core flock of STH in Yuncheng County, Shandong Province, China. The study was conducted in accordance with the Declaration of Helsinki, and approved by the Animal Ethics Committee of the Institute of Animal Sciences, Chinese Academy of Agriculture Science (No. IAS 2019–49).

Six WW genotyped ewes and six BB genotyped ewes were then selected based on the FecB genotype and reared under natural light conditions with ad libitum feeding. Vaginal sponges (progesterone 300 mg) (InterAg Co., Ltd., New Zealand) were used for 12 days of simultaneous estrus in the above sheep. We set the time of vaginal sponge removal at 0 h (around 10:00 am), and then euthanized 3 WW ewes and 3 BB ewes at 45 h (follicular phase, F) and 216 h (luteal phase, L), respectively, and obtained pituitary tissues, and then transferred those 12 tissues to –80°C for storage after snap–freezing in liquid nitrogen. Jugular venous blood samples were collected from sheep at the time of vaginal sponge removal (0 hours) and at the time of euthanasia (45 and 216 hours). The concentration of estradiol (E₂) was measured using an E₂ RIA kit (BNIBT, China) with a sensitivity of <2 pg/mL, while progesterone (P₄) levels were quantified using RIA kits (BNIBT, China) with a sensitivity of <0.1 ng/mL. To validate the precision of hormone detection, triplicate assays were conducted for each sample (Figure 1). The data for P₄ and E₂ were presented as mean ± standard error, and differences in hormone concentrations between genotypes at the same time were assessed using the T-TEST method, with statistical significance defined as P ≤ 0.05.

Figure 1

Figure 1 Proteomic and transcriptomic experiment sample group and process of data acquisition and analysis.

2.2 Peptide preparation and TMT labeling

The proteins from the pituitary gland tissues were extracted using the SDT lysis method (20), and protein concentrations were quantified using the Pierce™ BCA Protein Assay Kit (Thermo Fisher Scientific, USA) (21). Each sample was digested with trypsin using the filter aided sample preparation method, with approximately 200 µg of protein per sample (22). The peptide content was quantitated using ultraviolet spectrophotometry (optical density: 280) for subsequent proteomic experiments. Peptides (100 μg per sample) were labeled using the six plex tandem mass tag of TMT Mass Tagging kit (Thermo Fisher Scientific, USA) according to the manufacturer’s instructions. The labeling was performed on WW ewes at the luteal phase (WW_L), BB ewes at the luteal phase (BB_L), BB ewes at the follicular phase (BB_F), and WW ewes at the follicular phase (WW_F) using TMT–128, TMT–129, TMT–130, and TMT–131, respectively. Each group consisted of three biological replicates (Figure 1).

2.3 Peptide fractionation and Liquid chromatography–mass spectrometry/mass spectrometry analysis

The labeled peptides from 12 samples were mixed in equal volumes and subsequently lyophilized. A quantity of 100 μg of mixed peptides was dissolved in 300 µL of 0.1% trifluoroacetic acid. The solubilized peptide solution was then transferred to a High pH Reversed–Phase spin column (Thermo Fisher Scientific, USA), which had been pre–equilibrated with 0.1% trifluoroacetic acid and acetonitrile aqueous buffer. By employing gradient elution with high pH acetonitrile solutions of increasing concentration, ten fractions of column–bound peptides were obtained. Each eluted peptide sample was subjected to vacuum drying, followed by reconstitution with 12 μL of 0.1% formic acid, and the concentration were determined by measuring the absorbance at an optical density of 280.

Each sample was subjected to separation using an Easy nLC nanoflow liquid chromatography instrument (Thermo Fisher Scientific, USA). The column was initially equilibrated with 95% 0.1% formic acid in water (buffer A) employing a Thermo Scientific Acclaim PepMap100 (100 μm*2 cm, nanoViper C18, Thermo Fisher Scientific, USA). Subsequently, the samples were introduced into the column and separated utilizing a Thermo scientific EASY column (10 cm, inner diameter 75 μm, resin 3 μm, C18–A2, Thermo Fisher Scientific, USA). A 1–hour liquid gradient strategy was employed, wherein the concentration of buffer B, a 0.1% formic acid acetonitrile aqueous solution (comprising 84% acetonitrile), was gradually increased from 0% to 35% over the initial 50 minutes. Subsequently, the concentration of buffer B reached 100% within the time frame of 50 to 55 minutes. Finally, the concentration of buffer B was maintained at 100% for the remaining 5 minutes. The flow rate utilized was 300 nL/min.

The fractions of samples conducted through chromatography were analyzed using a Q–Exactive mass spectrometer (Thermo Fisher Scientific, USA). Detection was achieved through positive ions, and the acquisition of peptides and peptide fragments was performed by collecting Top20 mass spectrometry (MS) 2 scan fragments subsequent to each full scan. The parent ion scanning range for the full scan encompassed a range of 300 to 1800 m/z, with a resolution of 70000 at 200 m/z. The Automatic Gain Control (AGC) target was set at 1e6, while the Maximum injection time was limited to 50 ms and the Dynamic exclusion time was set to 60.0 s. The MS2 scan activation type employed is higher energy dissociation (HCD), with an isolation window of 2 m/z. The resolution achieved is 17500 at 200 m/z, while the Normalized Collision Energy utilized is 30eV. Additionally, the Underfill ratio is observed to be 0.1%.

2.4 The identification and relative quantification of proteomic and transcriptomic data

The MS/MS spectra data were subjected to analysis using Proteome Discoverer v.1.4 (Thermo Fisher Scientific, USA) and MASCOT engine v.2.2 (Matrix Science, UK). The protein reference database was constructed by translating transcriptome data from sheep used in the proteomic analysis (PRJNA782215) (23). The library searching parameters are detailed in Supplementary Table S1. Protein abundance differences were calculated after applying a filtering step, where proteins lacking abundance measurements in two or more samples from each group were excluded. The protein abundances were log2–transformed and normalized after filtration, and differences in protein abundance were determined using linear models for microarray data (LIMMA v 3.52.1) package (24). A factorial design was employed for the calculations, incorporating the factors of genotype (two levels: WW and BB) and estrus (two levels: F and L), with consideration given to interactions between the two factors. Proteins exhibiting significant differences in abundance were selected based on a P–value of less than 0.05 and a fold change in abundance comparison between the two groups greater than 1.2 or less than 0.83 (25).

The transcriptome data (PRJNA782215) was aligned and assembled with the sheep reference genome (Oar4.0: GCF_000298735.2) using HISAT2 (v.2.0.5, http://daehwankimlab.github.io/hisat2/) and String Tie (v.1.3.2d, https://ccb.jhu.edu/software/stringtie/). For the analysis of differentially expressed genes (DEGs), read counts were calculated using HTSeq (v0.6.1, https://lira.no-ip.org:8443/doc/python-htseq-doc/html/index.html). Subsequently, DESeq2 (v1.28.1) (26) was employed to normalize the counts and identify genes with differential expression. The model utilized was consistent with the one employed for identifying differentially abundant proteins (DAPs). The criteria for identifying differentially expressed genes included a P value less than 0.05, a corrected P value of less than 0.05, and a fold change greater than 2.0 between the two groups.

2.5 The integration of proteomic and transcriptomic data

We analyzed the repeatability and reliability of the test samples using Pearson’s correlation and Principal component analysis (PCA). The Data Integration Analysis for Biomarker discovery using Latent variable approaches for Omics studies (DIABLO) framework of mixOmics (v.6.20.0) was utilized to integrate transcriptomic and proteomic data in order to identify correlations between the two datasets. The DIABLO core algorithm is an extension of the partial least–squares method (27). Prior to data integration, transcript count data were normalized to log₂ (counts per million) (logCPM) using the edgeR package (v.3.38.1) (28). Transcripts with a logCPM value below 0 in more than 75% of the sample set were excluded. Similar filtering procedures were applied to the proteomic data during the identification of DAPs. The DIABLO design matrix incorporated estrus stage and genotype factors, and block linkage was 0.1. The determination of the ideal number of components and variables to incorporate in the final model was accomplished through the utilization of the “perf” and “tune.block.splsda” functions, and subsequently validated via 3 × 5–fold cross–validation (times). The model’s performance was evaluated by calculating the balanced error rate and overall misclassification error rate, while the model’s discrimination accuracy was assessed by calculating the area under the curve of the receiver operating characteristic curve (ROC).

2.6 Functional annotation and enrichment analysis of biomarkers

To explore the potential role of proteins in the luteal–follicular phase transition and the role of these biomarkers in FecB mutations affecting fertility, Gene Ontology (GO) categories (29) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (30) over–representation analyses were ran by clusterProfiler (v4.4.4) (31). The parameter of pAdjustMethod is false discovery rate (FDR). Top enriched KEGG pathways and GO terms were visualized using the two clusterProfiler package functions dotplot and cnetplot.

2.7 Quantitative analysis of selected proteins with PRM

To validate the protein abundance obtained by TMT quantification method, the abundance of randomly selected proteins in Pituitary were quantified using PRM analysis (32). The peptides used for PRM were shown in Supplementary Table S2. PRM assay was performed using peptides prepared for TMT analysis, approximately 1 μg of peptide from each sample, mixed with 20 fmol of the standard peptide (SAAGAFGPELSR, Thermo Fisher Scientific, USA). Reversed–phase chromatographic analysis was performed with the Easy nLC–1200 system (Thermo Fisher Scientific, USA) For the chromatographic analysis, the C18 Trap Column was 100 μm * 50 mm, with a resin of 5 μm, and the C18 Analytical Column was 75 μm * 200 mm, with a resin of 3 μm. The chromatographic separation exhibited a flow rate of 300 nl/min, with a gradient characterized by the following changes in concentration: buffer B experienced an increase from 5% to 10% within the initial 2 minutes, followed by a subsequent increase from 10% to 30% between 2 minutes and 45 minutes. From 45 minutes to 55 minutes, there was a further increase from 30% to 100%, and during the final 5 minutes, the concentration of buffer B remained constant at 100%. The fractions of samples conducted through chromatography were analyzed using a Q–Exactive mass spectrometer (Thermo Fisher Scientific, USA). Detection was achieved through positive ions, and the acquisition of peptides and peptide fragments was performed by collecting Top20 MS2 scan fragments subsequent to each full scan. The parent ion scanning range for the full scan encompassed a range of 300 to 1800 m/z, with a resolution of 60000 at 200 m/z. The AGC target was set at 3e6, while the Maximum injection time was limited to 200 ms. The MS2 scan activation type employed is higher energy dissociation (HCD), with an isolation window of 1.6Th. The resolution achieved is 30000 at 200 m/z, while the Normalized Collision Energy utilized is 27eV. The AGC target was set at 3e6, while the Maximum injection time was limited to 120 ms (33). The peptide were extracted and calculated by Skyline v.3.5.0 (34)

3 Results

3.1 Overview of pituitary proteomic data

Twelve pituitary samples were randomly divided into three experimental groups, with four samples in each group. Care was taken to ensure that each group included one sample from BB_F, BB_L, WW_F, and WW_L, thus rendering them appropriate for the six–plex TMT labeling strategy. The triple TMT experiments yielded 4896, 4860, and 5176 peptides, respectively. The subsequent integration of these three sets of results and their annotation ultimately led to the identification of 5836 distinct proteins, as depicted in Supplementary Figure S1A, for subsequent analysis. The analysis of protein molecular weight distribution showed that the molecular weights of 99% of proteins were between 6 and 312 kD (Supplementary Figure S1B). Approximately, the protein sequence coverage of 1545 proteins were > 30% (Supplementary Figure S1C). More than 1541 proteins included less than two unique peptides. The proteins XP_004018043.1, XP_014949338.1, XP_012026869.1, and XP_014952969.1 were identified to have over 100 unique peptides (Supplementary Table S3).

3.2 Results of PRM quantification

Albumin (ALB), calcium voltage–gated channel auxiliary subunit alpha 2 delta 2 (CACNA2D2), calcium/calmodulin dependent protein kinase I (CAMK1), glycoprotein hormones, alpha polypeptide (CGA), G protein subunit alpha o 1 (GNAO1), luteinizing hormone beta polypeptide (LHB), glutathione S–transferase Mu 1 (LOC101107401), solute carrier family 27 member 1 (SLC27A1), and alcohol dehydrogenase class–2 isozyme 2–like (LOC101110134) proteins were selected for validation of the differences in protein abundance using the PRM in the following subgroups: follicular versus luteal phase in BB genotype ewes (BB_F/L) group, follicular phase in BB ewes compared to WW ewes (BW_F) group, and luteal phase BB ewes compared with WW ewes (BW_L) groups, and follicular versus luteal phase in WW genotype ewes (WW_F/L) group. The differences in protein abundance obtained from TMT versus PRM in the above groups were then expressed as log2 (ratio) values. Comparison of these two methods showed a consistent trend in the differences in protein abundance between the different groups they detected (Figure 2, Supplementary Tables S2, S3).

Figure 2

Figure 2 Quantitative validation of the abundance of important proteins obtained by tandem mass tag (TMT) by parallel reaction monitoring (PRM). albumin (ALB), calcium voltage–gated channel auxiliary subunit alpha 2 delta 2 (CACNA2D2), calcium/calmodulin dependent protein kinase I (CAMK1), glycoprotein hormones, alpha polypeptide (CGA), G protein subunit alpha o 1 (GNAO1), luteinizing hormone beta polypeptide (LHB), glutathione S–transferase Mu 1 (LOC101107401), solute carrier family 27 member 1 (SLC27A1), and alcohol dehydrogenase class–2 isozyme 2–like (LOC101110134) were selected for validation of the differences in protein abundance using the PRM in the following subgroups: follicular versus luteal phase in BB genotype ewes (BB_F/L) group (A), follicular phase in BB ewes compared to WW ewes (BW_F) group (B), and luteal phase BB ewes compared with WW ewes (BW_L) groups (C), and follicular versus luteal phase in WW genotype ewes (WW_F/L) group (D). The log2(ratio) values were computed between the aforementioned groups. The protein abundance, as determined by the PRM method, exhibited similar trends to that of the TMT method for the selected proteins across all four groups.

3.3 The profile of the DAPs in proteomics data

The effects of genotype and estrus, as well as the interactions between them, were evaluated using the limFit and eBayes functions of the LIMMA package. Then, the contrast function was used to further analyze which proteins changed during the transition from the luteal to follicular phase in WWgenotype (FvsLWW) and BB genotype (FvsLBB) ewes, respectively. Additionally, the contrast function was used to determine which proteins responded differently (interactions) in BB vs. WW Following the thresholds described in the Material Methods, 170 DAPs were detected in the FvsLBB group. Out of these, 55 proteins were up–regulated during the follicular phase, while 115 proteins were down–regulated (Figure 3A, Supplementary Table S4). In the FvsLWW group during the follicular phase, 144 DAPs were screened, with 54 being up–regulated and 90 being down–regulated (Figure 3B, Supplementary Table S5). In the interaction group, 145 DAPs were detected, with 73 being up–regulated and 72 being down–regulated for BB genotypes relative to the wild type (Figure 3C, Supplementary Table S6).

Figure 3

Figure 3 Volcano plots were demonstrating the differentially abundant proteins (DAPs) in FvsLBB group (A), FvsLWW group (B) and Interaction group (C). The threshold for determining upregulated and downregulated DAPs was set at a foldchange higher than 1.2 or lower than 0.83, with a P value lower than 0.05.

We obtained the intersection of DAPs between the FvsLWW and FvsLBB groups and plotted the Venn diagram (35). The intersection of the FvsLWW and FvsLBB groups had 37 DAPs (Figure 4A). The abundance of a total of 277 unique DAPs in these two groups was then analyzed by hierarchical clustering heatmap, and the 37 DAPs in both groups showed the same trend (Figure 4B). We annotated the unique DAPs in these two groups with GO and KEGG enrichment (Figure 4C, Supplementary Tables S7, S8). Network clustering of the top 30 KEGG–enriched terms using cnetplot (Figure 4D). Pathways of the Oxytocin signaling pathway and GnRH secretion were linked together with CACNA1C, and this linkage contained 7 proteins, including LHB. The cnetplot also describes the corresponding network of three metabolically relevant KEGG–enriched terms (Glutathione metabolism, Fructose and mannose metabolism, and Propanoate metabolism) (Figure 4D), which are related to pituitary endocrine function.

Figure 4

Figure 4 The functional annotation and pathway enrichment of differentially abundant proteins (DAPs) during the luteal–follicular phase transition. (A) The Venn plot that illustrates the overlap of DAPs between the FvsLWW and FvsLBB groups. (B) Additionally, a hierarchical cluster analysis is performed on the DAPs in the FvsLWW and FvsLBB groups, with all protein abundance values normalized and scaled. (C) Furthermore, a dotplot is presented, showcasing the top 30 KEGG enrichment terms for the concatenated DAPs in the FvsLWW and FvsLBB groups. The pAdjust parameter used for statistical analysis is the false discovery rate (FDR). (D) Finally, the networks corresponding to the top 30 KEGG enrichment terms are visualized using cnetplot.

3.4 DIABLO multi–omics analysis to screen potential biomarkers associated with reproduction in STH sheep

The correlation coefficient of the gene expression levels between samples exceeded 0.85 in one group, suggesting a high degree of consistency and reliability in the sample selection (Supplementary Figure S3). Proteomic and transcriptomic data from the same ewes were integrated using the DIABLO framework. This supervised learning approach was used to screen for the best classification feature variables from multifactorial data containing follicular stage and genotype. In the follicular phase, there was significant variation in mRNA expression and protein abundance patterns between the two genotypes of type BB and WW sheep. However, during the luteal phase, the data patterns of the two omcis of type BB and WW sheep were found to be similar. This observation is illustrated in Supplementary Figure S2A. Furthermore, the proteomic and transcriptomic datasets exhibited a strong positive correlation, as evidenced by a Pearson’s correlation score of 0.98 on two components, as shown in Supplementary Figure S2B. Furthermore, the scatter plots (Supplementary Figure S2B) show that the biomarkers selected by integrating the two omics data can effectively discriminate between data from different estrus stages and different genotypes in the follicular phase. These results suggest that the FecB mutation primarily affects ovulation during the follicular phase. The “centroids.dist,” “mahalanobis.dist,” and “max.dist” are used to estimate the classification error rate. The findings derived from a 3 × 5–fold cross–validation demonstrate that the sPLS–DA model exhibits superior performance when ncomp = 3 (Supplementary Figure S2C). The circos plot (Figure 5) was employed to visualize the positive or negative correlation (Supplementary Table S9) between protein abundance and mRNA expression, utilizing ncomp = 1–2 and a correlation cutoff of 0.9. Among these protein–mRNA interaction pairs, 24 exhibited a positive correlation while 10 displayed a negative correlation. These proteins and mRNAs exhibited significant interactions, with three specific mRNAs (COL28A1, LOC106991550, and ST14) demonstrating connectivity. Notably, the LOC106991550 mRNA displayed interactions with a noteworthy total of 20 proteins.

Figure 5

Figure 5 A circos plot displays the positive (red lines) and negative (blue lines) correlations (r > 0.9) between features variables in the quadrants.

The plotLoadings function in the DIABLO package was used to visualize the degree of contribution of the featured variables in the three components obtained from the combined analysis. Among the variables in component 1 (Figure 6A), 90 were from proteomic data and 5 from transcriptomic data; among the variables in component 2 (Figure 6B), 20 were from proteomic data and 40 from transcriptomic data; and among the variables in component 3 (Figure 6C), 60 were from proteomic data and 5 from transcriptomic data (Supplementary Table S10). The auroc() function is utilized to calculate the Receiver Operating Characteristic (ROC) curve and Area Under the Curve (AUC) for comparing one group against the others. The AUC values, exceeding 0.8, obtained from two blocks of 2 components indicate that the sPLS–DA analysis conducted in this study is effective (Figure 6D). The 319 differential expressed genes (DEGs) (fold change (FC) > 2, pvalue < 0.05, Supplementary Table S11) and 121 DAPs obtained interaction group were then compared with the featured mRNAs (DIABLO_mRNA) and proteins (DIABLO_protein) selected by DIABLO (Supplementary Table S10). The four featured datasets obtained by these two methods differed significantly and did not intersect with each other. CALY and SPON1 were the three datasets intersected by DAP, DEG, and DIABLO_mrna. NUDT16, a DAP, also appeared in the DIABLO_mRNA and DIABLO_protein datasets, respectively. These results were visualized by the online tool ChiPlot (https://www.chiplot.online/) in Figure 6E.

Figure 6

Figure 6 The figure showcases the featured biomarkers and the accuracy of models in proteomic and transcriptomic data through the utilization of DIABLO. The loading vector weights of the prominent biomarkers in different omics are presented individually for component 1–3 (A–C). (D) Additionally, the Receiver operating characteristic area under the curve result higher than 0.85 demonstrates the high effectiveness of the DIABLO multi–omics joint analysis in identifying characteristic biomarkers. (E) Furthermore, the intersection between differentially abundant proteins (DAPs), differentially expressed genes (DEGs), and the featured biomarkers identified by the DIABLO method is illustrated using the upset plot. BB_F, BB genotype ewes at the follicular phase; BB_L, BB genotype ewes at the luteal phase; WW_F, WW genotype ewes at the follicular phase; WW_L, WW genotype ewes at the luteal phase.

The proteins and mRNAs from the four datasets were annotated according to KEGG enrichment analysis (Supplementary Tables S12–15). The top 30 KEGG pathways are presented in Figure 7. It is interesting to note that some pathways were related to energy and amino acid metabolic pathways related to reproduction processes, such as Cortisol synthesis and secretion, Arginine biosynthesis, Fructose and mannose metabolism, Pyruvate metabolism, and Tryptophan metabolism. These genes are also enriched for pituitary signaling molecules and interactions and related pathways of signal transduction, including the calcium signaling pathway and the cAMP signaling pathway. There were also several hormone–related pathways such as Growth hormone synthesis, secretion and action, Steroid hormone biosynthesis, and estrogen signaling pathway.

Figure 7

Figure 7 The top 30 KEGG pathways that exhibit enrichment in terms of differentially abundant proteins (DAPs) (A) and differentially expressed genes (DEGs) (B) within the interaction group, as well as the top 30 KEGG pathways enriched based on the biomarkers selected by DIABLO of protein (DIABLO_protein) (C) and mRNA (DIABLO_mRNA) (D) block.

4 Discussion

The pituitary is an important regulator involved in the hypothalamic–pituitary–gonad axis (HPG) axis and it orchestrate the complex neuroendocrine regulation of reproduction (36). FecB mutation may be involved in hormonal regulation and signal transduction during follicular development and ovulation (37). The combined analysis of proteome and transcriptome data provides an effective method to explain complex biological processes and screen valuable biomarkers for breeding (38, 39). Due to RNA–mediated transcription regulation and protein post–translational modification, the correlation between proteome and transcriptome data calculated by using Pearson correlation alone is very low, and it is difficult to effectively use the information of the two omics at the same time (40, 41). We choose DIABLO method, which mainly uses the principle of machine learning to solve the problem of multi–omics data integration mining with small sample size (42). Using this method to reduce the complexity of the two platforms of pituitary transcriptome and proteome data, we calculated the correlation of the two omics to be 0.98, and also screened important signature pathways.

The pituitary gland is an important gland for controlling the physiologic functions of reproduction, and it consists of the adenohypophysis and the neurohypophysis (43). In our proteomic results, we identified reproduction–related hormones secreted by the pituitary, FSH, LH, growth hormone, prolactin and thyroid–stimulating hormone (TSH). We also identified oxytocin synthesized by the hypothalamus and stored in the posterior pituitary (44). The process of follicular growth and ovulation is intricately regulated by interactions between the hypothalamus, pituitary gland, ovary, and uterus (37). The hypothalamus secretes gonadotropin-releasing hormone (GnRH) to stimulate the synthesis and release of FSH and LH from the pituitary gland (45). These hormones then act on the cells of the follicle to support its growth and maturation. As the follicle develops, the granulosa cells utilize androstenedione produced by the endometrium to synthesize E₂, which in turn inhibits the synthesis and release of FSH and LH from the pituitary gland through negative feedback regulation (46). During the follicular phase, there is a shift in the regulation of GnRH/LH by E2 from negative to positive feedback as the follicle matures, leading to the LH peak prior to ovulation (47). Our analysis revealed that E2 concentrations were significantly elevated in the BB genotype compared to the WW genotype during the follicular phase (Supplementary Figure S4), potentially accounting for the increased expression levels of Luteinizing Hormone Subunit Beta (LHB) in the BB genotype. In follicular, via TGFβ/BMP signaling pathway, mutant FecB gene carrier produced small antral follicles with a reduced number of granulosa cells exhibiting higher sensitive to FSH, leading to follicular maturation in advance as attested by precocious LH receptor expression (48).

Pituitary, as the intermediary organ of the reproductive axis, also regulate follicle–luteal transition by the synergy between endocrine and local auto/paracrine factors (49). The transition from follicle to corpus luteum was essential to the reproductive cycle and establish pregnancy by producing progesterone (50). During the luteal–follicular phase transition, our Network clustering analysis of the DAP revealed that Voltage-dependent calcium channel L type alpha-1C (CACNA1C) serves as a point of convergence between the Oxytocin signaling pathway and the GnRH secretion pathway (Figure 4). The CACNA1C gene, recognized for its role in calcium channel activity, has been pinpointed as a significant candidate gene for the initiation of puberty in Jining grey goats (51). Furthermore, evidence suggests a strong correlation between CACNA1C and GnRH, which collaborate to facilitate ERK activation and subsequent elevation in FSH and LH secretion (52). These observations lead us to propose that CACNA1C plays a pivotal role in the transition from the luteal to the follicular phase. The Oxytocin signaling pathway exhibits enrichment of calcium/calmodulin-dependent protein kinase I (CAMKI), as well as the CAMKI and CACNA1C genes, which are enriched in the calcium signaling pathway. Intracellular Ca²⁺ concentration serves as a crucial signaling molecule in the regulation of exocytosis, controlling the release of neurotransmitters and endocrine hormones (53). The release of pituitary gonadotropins is stimulated by GnRH from the hypothalamus through a Ca²⁺ dependent mechanism (54). Research has demonstrated that the activation of CaMKI by GnRH is essential for the derepression of the FSHB gene through the phosphorylation of various class IIa HDACs. Conversely, the derepression of the LHB-subunit gene does not rely on CaMKI activation (55).

After annotation and enrichment of DEGs of interaction group. we identified the Neuroactive ligand–receptor interaction pathway. This pathway was also enriched in our previous multi–omics analysis of the hypothalamic in different genotypes of FecB mutation (56). This pathway was identified in both the hypothalamus and pituitary in comparative transcriptomic studies of the mammalian and poultry HPG axis (53). Proteomic analysis of the hypothalamus and pituitary revealed that Gamma–aminobutyric acid receptor (GABR) related genes were enriched in the Neuroactive ligand–receptor interaction pathway, including GABRA2 in the pituitary, GABRA1 and GABRB2 in the hypothalamus. Gamma–aminobutyric acid (GABA) is a major neurotransmitter in the central nervous system, and it acts mainly through three receptors, GABA_A, GABA_B and GABA_C. GABA_A is a member of the ligand–gated Cl– channel family, while the GABA_B receptor is coupled to G proteins (57). The GABRB1 gene identified by the DIABLO method is an important featured variable in the pituitary transcriptome, and it is highly expressed in follicular BB ewes than in wild–type, this expressed trend that is consistent with the LHB abundance trend in the proteome. GABA can regulate BMP2 gene expression through the BMP pathway by activating the GABAB receptor (58). BMP2 is a potential ligand for BMPR1B, and results from in vitro experiments in sheep granulosa cells suggest that BMP2 affects estrogen secretion (48). It has been shown that GABA stimulates LH secretion from pituitary cells (59). The knockout of GABA receptors in the hypothalamus may influence the negative feedback regulation of estrogen on the gonadotropin-releasing hormone (GnRH) in mice, consequently impacting the secretion of luteinizing hormone (LH) (60). In summary, the pituitary gland functions as a central component in the HPG axis, integrating signals from the GABA system to regulate ovulation through the modulation of LH or GnRH/LH secretion and interactions with the BMP pathway.

5 Conclusions

The investigation of protein abundance in the pituitary of sheep with varying FecB genotypes at distinct estrus stages revealed the existence of anterior pituitary hormones (FSH, LH, PRL, and TSH). Mutations in FecB have the potential to impact follicular development through the regulation of FSH and LH, consequently leading to alterations in the number of ovulations. CACNA1C and CaMKI are significant marker genes in the pituitary gland that play a crucial role in the regulation of LH and FSH secretion through GnRH. The pituitary integrates signals from the GABA system to modulate ovulation by influencing LH or GnRH/LH secretion and interacting BMP pathway. The integrated analysis of the pituitary gland has yielded novel insights into the significant endocrine regulatory function performed by the pituitary gland in the physiological process of FecB mutation impacting ovulation.

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: National Center for Biotechnology Information Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra/ Accession no. PRJNA782215) and the OMIX, China National Center for Bioinformation / Beijing Institute of Genomics, Chinese Academy of Sciences (61, 62) (https://ngdc.cncb.ac.cn/omix/ Accession no. OMIX006109).

Ethics statement

The animal study was approved by the Animal Ethics Committee of the Institute of Animal Sciences, Chinese Academy of Agriculture Science (No. IAS 2019–49). The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

XW: Validation, Software, Methodology, Investigation, Formal analysis, Conceptualization, Writing – review & editing, Writing – original draft. XG: Validation, Methodology, Investigation, Funding acquisition, Writing – review & editing. XH: Methodology, Investigation, Writing – review & editing. RD: Project administration, Investigation, Writing – review & editing. XZ: Resources, Investigation, Writing – review & editing. JZ: Resources, Investigation, Writing – review & editing. MC: Supervision, Resources, Project administration, Funding acquisition, Conceptualization, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This research was funded by the National Natural Science Foundation of China (31861143012, 31902150), the Earmarked Fund for China Agriculture Research System of MOF and MARA (CARS–38), the Central Public-Interest Scientific Institution Basal Research Fund (Y2017JC24, 2018ywf–yb–2), the Agricultural Science and Technology Innovation Program of China (CAAS–ZDRW202106 and ASTIP–IAS13), Jilin Province Youth Growth Technology Project (20230508122RC), the Natural Science Foundation of Tianjin (23JCYBJC00260).

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/fendo.2024.1417530/full#supplementary-material

Supplementary Figure 1 | Distributions of molecular weight, Protein sequence coverage, and number of unique peptides of proteins identified in the pituitary. (A) Distribution of protein molecular weights (kDa). (B) Distribution of protein sequence coverage (%). (C) Venn diagram displaying the overlap of protein identification by proteomic TMT 6-plex experiments.

Supplementary Figure 2 | Illustration of N-integration supervised analysis with DIABLO. (A) Sample plot per data set highlighting the integrated protein and mRNA targets. (B) Diagnostic scatterplot displaying components 1–2 of each dataset (protein and mRNA), (C) Classification performance per component (overall and Balanced Error Rate) for prediction distances (“centroids.dist”, “mahalanobis.dist” and “max.dist”) using repeated stratified cross-validation (3 x 5 fold).

Supplementary Figure 3 | The Pearson’s correlation (A) and Principal component analysis (B) of all samples.

Supplementary Figure 4 | Mean ± SEM of serum P4 and E2, and concentration at vaginal sponge removal (0 h) and at the time of euthanasia (45 h and 216 h) in different FecB genotype ewes. Values with different lowercase letters indicate significant differences between genotypes at the same time. (P ≤ 0.05). (A) Mean ± SEM of serum P4 (pg/mL), and (B) Mean ± SEM of serum E2 (pg/mL).

References

1. Mulsant P, Lecerf F, Fabre S, Schibler L, Monget P, Lanneluc I, et al. Mutation in bone morphogenetic protein receptor-IB is associated with increased ovulation rate in Booroola Mérino ewes. Proc Natl Acad Sci United States America. (2001) 98:5104–9. doi: 10.1073/pnas.091577598

CrossRef Full Text | Google Scholar

2. Souza CJ, Macdougall C, Macdougall C, Campbell BK, Mcneilly AS, Baird DT. The Booroola (FecB) phenotype is associated with a mutation in the bone morphogenetic receptor type 1 B (BMPR1B) gene. J Endocrinol. (2001) 169:1–6. doi: 10.1677/joe.0.169r001

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Wilson T, Wu XY, Juengel JL, Ross IK, Lumsden JM, Lord EA, et al. Highly prolific Booroola sheep have a mutation in the intracellular kinase domain of bone morphogenetic protein IB receptor (ALK-6) that is expressed in both oocytes and granulosa cells. Biol Reprod. (2001) 64:1225–35. doi: 10.1095/biolreprod64.4.1225

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Zhou S, Ding Y, Liu J, Liu Y, Zhao X, Li G, et al. Highly efficient generation of sheep with a defined FecB(B) mutation via adenine base editing. Genet Sel Evol. (2020) 52:35. doi: 10.1186/s12711-020-00554-6

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Zhou S, Yu H, Zhao X, Cai B, Ding Q, Huang Y, et al. Generation of gene-edited sheep with a defined Booroola fecundity gene (FecB(B)) mutation in bone morphogenetic protein receptor type 1B (BMPR1B) via clustered regularly interspaced short palindromic repeat (CRISPR)/CRISPR-associated (Cas) 9. Reprod Fertil Dev. (2018) 30:1616–21. doi: 10.1071/RD18086

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Wang X, Guo X, He X, Liu Q, Di R, Hu W, et al. Effects of fecB mutation on estrus, ovulation, and endocrine characteristics in small tail han sheep. Front Vet Sci. (2021) 8:709737. doi: 10.3389/fvets.2021.709737

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Guo X, Wang X, Di R, Liu Q, Hu W, He X, et al. Metabolic effects of FecB gene on follicular fluid and ovarian vein serum in sheep (Ovis aries). Int J Mol Sci. (2018) 19:539. doi: 10.3390/ijms19020539

PubMed Abstract | CrossRef Full Text | Google Scholar

8. McNatty KP, Hudson NL, Shaw L, Condell LA, Ball K, Seah SL, et al. GnRH-induced gonadotrophin secretion in ovariectomized Booroola ewes with hypothalamic-pituitary disconnection. J Reprod Fertil. (1991) 91:583–92. doi: 10.1530/jrf.0.0910583

PubMed Abstract | CrossRef Full Text | Google Scholar

9. McNatty KP, Hudson NL, Lun S, Heath DA, Shaw L, Condell L, et al. Gonadotrophin-releasing hormone and the control of ovulation rate by the FecB gene in Booroola ewes. J Reprod Fertil. (1993) 98:97–105. doi: 10.1530/jrf.0.0980097

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Juengel JL, Davis GH, McNatty KP. Using sheep lines with mutations in single genes to better understand ovarian function. Reproduciton. (2013) 146:R111–23. doi: 10.1530/REP-12-0509

CrossRef Full Text | Google Scholar

11. Plant TM. 60 YEARS OF NEUROENDOCRINOLOGY: The hypothalamo-pituitary–gonadal axis. J Endocrinol. (2015) 226:T41–54. doi: 10.1530/JOE-15-0113

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Li X, Li C, Wei J, Ni W, Xu Y, Yao R, et al. Comprehensive expression profiling analysis of pituitary indicates that circRNA participates in the regulation of sheep estrus. Genes (Basel). (2019) 10:90. doi: 10.3390/genes10020090

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Zheng J, Wang Z, Yang H, Yao X, Yang P, Ren C, et al. Pituitary Transcriptomic Study Reveals the Differential Regulation of lncRNAs and mRNAs Related to Prolificacy in Different FecB Genotyping Sheep. Genes (Basel). (2019) 10:157. doi: 10.3390/genes10020157

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Li X, Li C, Wureli H, Ni W, Zhang M, Li H, et al. Screening and evaluating of long non-coding RNAs in prenatal and postnatal pituitary gland of sheep. Genomics. (2020) 112:934–42. doi: 10.1016/j.ygeno.2019.06.009

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Li X, Li C, Xu Y, Yao R, Li H, Ni W, et al. Analysis of pituitary transcriptomics indicates that lncRNAs are involved in the regulation of sheep estrus. Funct Integr Genomics. (2020) 20:563–73. doi: 10.1007/s10142-020-00735-y

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Yang H, Ma J, Wang Z, Yao X, Zhao J, Zhao X, et al. Genome-wide analysis and function prediction of long noncoding RNAs in sheep pituitary gland associated with sexual maturation. Genes (Basel). (2020) 11:320. doi: 10.3390/genes11030320

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Liu Q, Di R, Ren C, He X, Wang X, Xia Q, et al. Screening of Differentially Expressed Genes and miRNAs in Hypothalamus and Pituitary Gland of Sheep under Different Photoperiods. Genes (Basel). (2022) 13:1091. doi: 10.3390/genes13061091

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Wan Z, Yang H, Cai Y, Ma J, Cheng P, Wang Z, et al. Comparative transcriptomic analysis of Hu sheep pituitary gland prolificacy at the follicular and luteal phases. Genes (Basel). (2022) 13:440. doi: 10.3390/genes13030440

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Lee AH, Shannon CP, Amenyogbe N, Bennike TB, Diray-Arce J, Idoko OT, et al. The, Dynamic molecular changes during the first week of human life follow a robust developmental trajectory. Nat Commun. (2019) 10:1092. doi: 10.1038/s41467-019-08794-x

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Zhu Y, Xu H, Chen H, Xie J, Shi M, Shen B, et al. Proteomic analysis of solid pseudopapillary tumor of the pancreas reveals dysfunction of the endoplasmic reticulum protein processing pathway. Mol Cell Proteomics. (2014) 13:2593–603. doi: 10.1074/mcp.M114.038786

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Smith PK, Krohn RI, Hermanson GT, Mallia AK, Gartner FH, Provenzano MD, et al. Measurement of protein using bicinchoninic acid. Anal Biochem. (1985) 150:76–85. doi: 10.1016/0003-2697(85)90442-7

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Wiśniewski JR, Zougman A, Nagaraj N, Mann M. Universal sample preparation method for proteome analysis. Nat Methods. (2009) 6:359–62. doi: 10.1038/nmeth.1322

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Chen S, Guo X, He X, Di R, Zhang X, Zhang J, et al. Transcriptome analysis reveals differentially expressed genes and long non-coding RNAs associated with fecundity in sheep hypothalamus with different FecB genotypes. Front Cell Dev Biol. (2021) 9:633747. doi: 10.3389/fcell.2021.633747

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:e47–7. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Zhang Z, Tang J, Di R, Liu Q, Wang X, Gan S, et al. Identification of prolificacy-related differentially expressed proteins from sheep (Ovis aries) hypothalamus by comparative proteomics. Proteomics. (2019) 19:1900118. doi: 10.1002/pmic.201900118

CrossRef Full Text | Google Scholar

26. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. (2014) 15:550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Singh A, Shannon CP, Gautier B, Rohart F, Vacher M, Tebbutt SJ, et al. DIABLO: an integrative approach for identifying key molecular drivers from multi-omics assays. Bioinformatics. (2019) 35:3055–62. doi: 10.1093/bioinformatics/bty1054

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. (2009) 26:139–40. doi: 10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene Ontology: tool for the unification of biology. Nat Genet. (2000) 25:25–9. doi: 10.1038/75556

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. (2000) 28:27–30. doi: 10.1093/nar/28.1.27

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. . Innovation (Camb). (2021) 2:100141. doi: 10.1016/j.xinn.2021.100141

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Peterson AC, Russell JD, Bailey DJ, Westphall MS, Coon JJ. Parallel reaction monitoring for high resolution and high mass accuracy quantitative, targeted proteomics. Mol Cell Proteomics. (2012) 11:1475–88. doi: 10.1074/mcp.O112.020131

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Tang J, Hu W, Chen S, Di R, Liu Q, Wang X, et al. The genetic mechanism of high prolificacy in small tail han sheep by comparative proteomics of ovaries in the follicular and luteal stages. J Proteomics. (2019) 204:103394. doi: 10.1016/j.jprot.2019.103394

PubMed Abstract | CrossRef Full Text | Google Scholar

34. MacLean B, Tomazela DM, Shulman N, Chambers M, Finney GL, Frewen B, et al. Skyline: an open source document editor for creating and analyzing targeted proteomics experiments. Bioinformatics. (2010) 26:966–8. doi: 10.1093/bioinformatics/btq054

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Bardou P, Mariette J, Escudié F, Djemiel C, Klopp C. jvenn: an interactive Venn diagram viewer. BMC Bioinf. (2014) 15:293. doi: 10.1186/1471-2105-15-293

CrossRef Full Text | Google Scholar

36. Stamou MI, Cox KH, Crowley WF Jr. Discovering genes essential to the hypothalamic regulation of human reproduction using a human disease model: adjusting to life in the “-omics” Era. Endocr Rev. (2015) 36:603–21. doi: 10.1210/er.2015-1045

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Juengel JL. How the quest to improve sheep reproduction provided insight into oocyte control of follicular development. J R Soc New Z. (2018) 48:143–63. doi: 10.1080/03036758.2017.1421238

CrossRef Full Text | Google Scholar

38. Alvira-Iraizoz F, Gillard BT, Lin P, Paterson A, Pauža AG, Ali MA, et al. Multiomic analysis of the Arabian camel (Camelus dromedarius) kidney reveals a role for cholesterol in water conservation. Commun Biol. (2021) 4:779. doi: 10.1038/s42003-021-02327-3

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Sun Z, Liu Y, He X, Di R, Wang X, Ren C, et al. Integrative proteomics and transcriptomics profiles of the oviduct reveal the prolificacy-related candidate biomarkers of goats (Capra hircus) in estrous periods. Int J Mol Sci. (2022) 23:14888. doi: 10.3390/ijms232314888

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Pan Z, Zeng Y, An J, Ye J, Xu Q, Deng X. An integrative analysis of transcriptome and proteome provides new insights into carotenoid biosynthesis and regulation in sweet orange fruits. J Proteomics. (2012) 75:2670–84. doi: 10.1016/j.jprot.2012.03.016

PubMed Abstract | CrossRef Full Text | Google Scholar

41. de Sousa Abreu R, Penalva LO, Marcotte EM, Vogel C. Global signatures of protein and mRNA expression levels. Mol Biosyst. (2009) 5:1512–26. doi: 10.1039/b908315d

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Fan J, Lv J. A selective overview of variable selection in high dimensional feature space. Stat Sin. (2010) 20:101–48.

PubMed Abstract | Google Scholar

43. Vázquez-Borrego MC, Gahete MD, Martínez-Fuentes AJ, Fuentes-Fayos AC, Castaño JP, Kineman RD, et al. Multiple signaling pathways convey central and peripheral signals to regulate pituitary function: Lessons from human and non-human primate models. Mol Cell Endocrinol. (2018) 463:4–22. doi: 10.1016/j.mce.2017.12.007

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Dorton AM. The pituitary gland: embryology, physiology, and pathophysiology. Neonatal Netw. (2000) 19:9–17. doi: 10.1891/0730-0832.19.2.9

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Clarke IJ, Arbabi L. New concepts of the central control of reproduction, integrating influence of stress, metabolic state, and season. Domest Anim Endocrinol. (2016) 56:S165–79. doi: 10.1016/j.domaniend.2016.03.001

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Webb R, Campbell BK. Development of the dominant follicle: mechanisms of selection and maintenance of oocyte quality. Soc Reprod Fertil Suppl. (2007) 64:141–63. doi: 10.5661/RDR-VI-141

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Caraty A, Evans NP, Fabre-Nys CJ, Karsch EJ. The preovulatory gonadotrophin-releasing hormone surge: a neuroendocrine signal for ovulation. J Reprod Fertil Suppl. (1995) 49:245–55.

PubMed Abstract | Google Scholar

48. Fabre S, Pierre A, Mulsant P, Bodin L, Di Pasquale E, Persani L, et al. Regulation of ovulation rate in mammals: contribution of sheep genetic models. Reprod Biol Endocrinol. (2006) 4:20. doi: 10.1186/1477-7827-4-20

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Joseph C, Hunter MG, Sinclair KD, Robinson RS. The expression, regulation and function of secreted protein, acidic, cysteine-rich in the follicle–luteal transition. Reproduction. (2012) 144:361–72. doi: 10.1530/REP-12-0099

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Robinson R, Hammond A, Wathes D, Hunter M, Mann G. Corpus luteum–endometrium–embryo interactions in the dairy cow: underlying mechanisms and clinical relevance. Reprod Domest Anim. (2008) 43:104–12. doi: 10.1111/j.1439-0531.2008.01149.x

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Su F, Guo X, Wang Y, Wang Y, Cao G, Jiang Y. Genome-wide analysis on the landscape of transcriptomes and their relationship with DNA methylomes in the hypothalamus reveals genes related to sexual precocity in jining gray goats. Front Endocrinol (Lausanne). (2018) 9:501. doi: 10.3389/fendo.2018.00501

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Dang AK, Murtazina DA, Magee C, Navratil AM, Clay CM, Amberg GC. GnRH evokes localized subplasmalemmal calcium signaling in gonadotropes. Mol Endocrinol. (2014) 28:2049–59. doi: 10.1210/me.2014-1208

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Li X, Lin B, Zhang X, Shen X, Ouyang H, Wu Z, et al. Comparative transcriptomics in the hypothalamic-pituitary-gonad axis of mammals and poultry. Genomics. (2022) 114:110396. doi: 10.1016/j.ygeno.2022.110396

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Bates MD, Conn PM. Calcium mobilization in the pituitary gonadotrope: Relative roles of intra-and extracellular sources. Endocrinology. (1984) 115:1380–5. doi: 10.1210/endo-115-4-1380

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Lim S, Luo M, Koh M, Yang M, M.N. bin Abdul Kadir JH, Ye Z, et al. Distinct mechanisms involving diverse histone deacetylases repress expression of the two gonadotropin β-subunit genes in immature gonadotropes, and their actions are overcome by gonadotropin-releasing hormone. Mol Cell Biol. (2007) 27:4105–20. doi: 10.1128/MCB.00248-07

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Wang X, Guo X, He X, Di R, Zhang X, Zhang J, et al. Integrated proteotranscriptomics of the hypothalamus reveals altered regulation associated with the FecB mutation in the BMPR1B gene that affects prolificacy in small tail han sheep. Biology. (2023) 12:72. doi: 10.3390/biology12010072

CrossRef Full Text | Google Scholar

57. Owens DF, Kriegstein AR. Is there more to gaba than synaptic inhibition? Nat Rev Neurosci. (2002) 3:715–27. doi: 10.1038/nrn919

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Muhammad SI, Maznah I, Mahmud R, Zuki AB, Imam MU. Upregulation of genes related to bone formation by γ-amino butyric acid and γ-oryzanol in germinated brown rice is via the activation of GABAB-receptors and reduction of serum IL-6 in rats. Clin Interv Aging. (2013) 8:1259–71. doi: 10.2147/CIA

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Mijiddorj T, Kanasaki H, Sukhbaatar U, Oride A, Kyo S. DS1, a delta subunit-containing GABAA receptor agonist, increases gonadotropin subunit gene expression in mouse pituitary gonadotrophs1. Biol Reprod. (2015) 92:45. doi: 10.1095/biolreprod.114.123893

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Lee K, Porteous R, Campbell RE, Lüscher B, Herbison AE. Knockdown of GABAA receptor signaling in gnRH neurons has minimal effects upon fertility. Endocrinology. (2010) 151:4428–36. doi: 10.1210/en.2010-0314

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Chen T, Chen X, Zhang S, Zhu J, Tang B, Wang A, et al. The genome sequence archive family: toward explosive data growth and diverse data types. Genomics Proteomics Bioinf. (2021) 19:578–83. doi: 10.1016/j.gpb.2021.08.001

CrossRef Full Text | Google Scholar

62. Members C-N. and partners, database resources of the national genomics data center, China national center for bioinformation in 2022. Nucleic Acids Res. (2021) 50:D27–38. doi: 10.1093/nar/gkab951

CrossRef Full Text | Google Scholar