Skip to main content
Download PDF

Transcriptomics and metabolomics of blood, urine and ovarian follicular fluid of yak at induced estrus stage

BMC Genomics volume 25, Article number: 201 (2024) Cite this article

A Correction to this article was published on 08 April 2024

This article has been updated

Abstract

To gain a deeper understanding of the metabolic differences within and outside the body, as well as changes in transcription levels following estrus in yaks, we conducted transcriptome and metabolome analyses on female yaks in both estrus and non-estrus states. The metabolome analysis identified 114, 13, and 91 distinct metabolites in urine, blood, and follicular fluid, respectively. The Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis highlighted an enrichment of pathways related to amino acid and lipid metabolism across all three body fluids. Our transcriptome analysis revealed 122 differentially expressed genes within microRNA (miRNA) and 640 within long non-coding RNA (lncRNA). Functional enrichment analysis of lncRNA and miRNA indicated their involvement in cell signaling, disease resistance, and immunity pathways. We constructed a regulatory network composed of 10 lncRNAs, 4 miRNAs, and 30 mRNAs, based on the targeted regulation relationships of the differentially expressed genes. In conclusion, the accumulation of metabolites such as amino acids, steroids, and organic acids, along with the expression changes of key genes like miR-129 during yak estrus, provide initial insights into the estrus mechanism in yaks.

Peer Review reports

Introduction

Yak is a significant source of income for herdsmen in the Qinghai Plateau region, serving as the main economic breed [1]. They provide daily necessities such as meat, milk, and fur to local herdsmen [2]. Additionally, yaks were historically used as a means of transportation in ancient times when modern transportation was not developed. However, the development of the yak industry is hindered by breeding issues such as late sexual maturity, long calving interval, and unclear estrus period [3]. After calving, yaks typically exhibit estrus in the second year or later, rather than the year of calving, resulting in a prolonged postpartum period of low fertility [4]. The seasonal breeding characteristics of yaks also contribute to their low breeding efficiency, limiting the growth of the yak breeding industry [5]. Therefore, studying the yak estrus mechanism to improve breeding efficiency and promote yak production is of significant importance.

Estrus and ovulation in female animals are regulated by the neuroendocrine and reproductive hormones in the hypothalamic-pituitary–gonadal (HPG) axis system [6]. Within the HPG axis, gonadotropin-releasing hormone (GnRH) neurons are primarily located in the hypothalamus and act on the anterior pituitary to stimulate the release of luteinizing hormone (LH) and follicle-stimulating hormone (FSH) [7]. FSH and LH hormones then reach the ovaries via blood circulation, promoting ovarian development and inducing ovarian estrogen production [8]. Studies have shown that the regulation of hypothalamic-pituitary processes can be influenced by the use of exogenous GnRH, which in turn can affect physiological processes related to reproduction [9].

Omics, including genomics, transcriptomics, and other omics, are used to study the genome, transcriptome, and their objects, respectively [10]. In essence, the addition of "omics" to a molecular term implies a comprehensive or global assessment of a set of molecules [11]. Metabolomics, a branch of omics, involves the quantitative analysis of all metabolites in the body and aims to investigate the relationship between metabolites and physiological and pathological changes. As a relatively recent scientific research method, metabolomics has gained significant attention in yak research, particularly in the study of yak milk [12], yak blood [13], and yak meat [14]. Some studies have also explored the potential mechanisms of yak estrus through metabolomics. For example, Jia Zhou et al. used metabolomics and transcriptomics to study the follicles of yaks supplemented with NCG in their diets, revealing that an increased number of large ovarian follicles (diameter > 10 mm) were associated with 889 genes and 94 metabolites, which may contribute to the development of these follicles [15].

Understanding the molecular mechanisms underlying yak estrus development has been a prominent focus in yak genetic breeding research. Zhou et al. conducted metabolomic and transcriptomic analyses to identify the gene regulatory network and metabolic processes related to steroid hormone biosynthesis in ovarian tissue during follicular development, providing valuable insights into the regulation of metabolite changes and follicular development in yak estrus [15]. Hormone levels can objectively reflect the estrus state of yaks, making it essential to investigate hormone changes at the molecular level to study the estrus mechanism in yaks. Granulosa cells (GCs), as a crucial organ involved in animal estrus, secrete steroid hormones such as estrogen and progesterone, creating a crucial microenvironment for follicular growth [16]. In recent years, lncRNAs and miRNAs have emerged as hot topics in research. Constructing an interaction network among lncRNAs, miRNAs, and target mRNAs can shed light on their potential roles in GCs proliferation and differentiation.

In this study, we conducted a comprehensive analysis of the effects of sympathetic estrus treatment on yak estrus using transcriptomic and metabolomic data, aiming to uncover the underlying molecular mechanisms of the treatment on yak estrus.

Materials and methods

Ethics statement

The terms of all individual laboratory animals are approved by the Animal Policy and Welfare Committee of Northwest A&F University College (FAPWC-NWAFU, Protocol number WAFAC1008).

Materials

For this study, we selected 10 non-pregnant female yaks of similar age and randomly divided them into two groups with 5 yaks per group. We employed a synchronous estrus treatment called Ovsynch, On day 0 of the experiment, the animals were injected with 100 µg of gonadotropin-releasing hormone (GnRH) and administered a progesterone vaginal insert. On day 7, they received an injection of 0.5 mg PG, and the vaginal insert was removed. On day 9, another 100 µg of GnRH was injected, and the yaks' estrus condition was observed [17].

Ten female yaks of the same age were raised in the yak breeding base in Gangcha County, Qinghai Province. After calving in May 2021, the postpartum yaks underwent a three-month isolation and weaning process starting in August 2021. Following a 10–12 day period of weaning and isolation for the calves, an estrus induction procedure was implemented. Then, we randomly selected five postpartum female yaks to undergo the Ovsynch protocol. We utilized an ovulation meter to determine the estrus status of yaks, where estrus values ranged from 200 to 300. Simultaneously, in conjunction with ultrasound diagnosis, we observed whether dominant follicles were present on the ovaries. If such follicles were identified, it indicated the occurrence of estrus. On the 10th day of the Ovsynch protocol, the estrus in yaks was confirmed, and the slaughtering process was carried out at Qinghai Xiahua Halal Meat Products Co., Ltd. We obtained ovaries from five yaks in estrus through the slaughtering process. Within 2 h of collection, the ovaries were placed in 35℃ physiological saline (0.9% w/v) and transported to the laboratory. After sterilization, follicular fluid was aspirated from the ovaries using a needle and syringe. Subsequently, granulosa cells were obtained by centrifugation at 1500 rpm [18]. Urine and blood samples were collected from live animals. We collected blood, urine, and follicular fluid from the estrus group (n = 5) and the control group (n = 5) under the same feeding conditions. Samples were rapidly cooled with liquid nitrogen and stored in a refrigerator at -80℃ for metabolomics sequencing and analysis. Additionally, we collected follicular granulosa cells from the control group (n = 3) and the estrus group (n = 3) and used the same storage method for RNA-seq.

Untargeted metabolomics study

The metabolic profiling of urine, follicular fluidand bloodsamples was performed using the 2777C UPLC system (Waters, UK) coupled to the Xevo G2-XS QTOF (Waters, UK). ACQUITY UPLC HSS T3 column (100 mm*2.1 mm, 1.8 μm, Waters, UK) was used for separation.

ACQUITY UPLC HSS T3 column (100 mm*2.1 mm, 1.8 μm, Waters, UK) was used for separation. The column temperature was 40 °C and the flow rate was 0.5 ml/min. The mobile phase A consisted of water and 0.1% formic acid. B mobile phase consists of acetonitrile and 0.1% formic acid. The metabolites were eluted with the following gradient: 0–1 min, 99% mobile phase A; 1–3 min, 1–15% mobile phase B; 3–6 min, 15–50% mobile phase B; 6–9 min, 50–95% mobile phase B; 9–10 min, 95% mobile phase B; At 10.1–12 min, 99% mobile phase A was obtained. The loading volume of each sample was 5 μl.The small molecules eluted from the chromatographic column were respectively collected in positive and negative ion mode by high resolution tandem mass spectrometry Xevo G2-XS QTOF (Waters, UK).

Metabolomics data analysis

We used Progenesis QI (version 2.2) to perform peak alignment, peak extraction, normalization, deconvolution, and compound identification. To calibrate the real sample signal, we utilized information from the QC sample and performed local polynomial regression fitting [19]. OmicStudio (https://www.omicstudio.cn/tool) was used for Principal component analysis (PCA), Partial Least Squares Discriminant (PLSDA), and volcano analysis. When PCA fell short in effectively distinguishing between inter-group samples, PLS-DA emerged as a compelling alternative. Going beyond its capacity for dimensionality reduction, PLS-DA also facilitated the efficient separation of sample categories. In addition to its role in reducing data dimensions, PLS-DA empowered the prediction of sample categories, making it well-suited for classification purposes. By constructing a classification prediction model, PLS-DA delved deeper into identifying the affiliations of additional samples, a task that exploratory PCA methods were unable to accomplish. We determined significantly regulated metabolites between groups by applying the following criteria: |log2(fold change)|≥ 1 and P valve < 0.05. The identified metabolites were then enriched using KEGG and Gene Ontology (GO) enrichment analysis on the MetaboAnalyst website (https://www.metaboanalyst.ca/).

RNA extraction, sequencing of small RNA and LncRNA

We extracted total RNA using Trizol reagent (Ambion, United States) as per the manufacturer's instructions. Ligation PCR products were sequenced on the BGISEQ-500 platform (BGI Group, Shenzhen, China).

MiRNA analysis and identification

To remove low-quality and adaptor sequences, we processed the data using SOAPnuke [20]. High-quality data obtained from the previous step were mapped to the Bos Taurus reference genome (ARS-UCD1.2_Btau5.0.1Y) using HISAT2 [21]. Due to the absence of miRNA sequences for yaks in the miRNA database, and considering the need for a consistent species in data analysis, the reference genome of yellow cattle was chosen. The miRNA sequences for yellow cattle were used for analysis. Data were compared with Rfam database (https://rfam.org/) to filter rRNA, scRNA, snoRNA, snRNA, tRNA and other types of non-coding small RNAs. We calculated the expression of miRNA using Bos Taurus in miRbase release 22.1 as a reference [22, 23]. Read counts were normalized by the total count of mappable reads of each sample, referred to as RPM, for unbiased comparisons among samples. Normalized expressions were further standardized by z-scores for cluster and PCA [24]. To study the expression of miRNA in the control group and the experimental group, we identified DE miRNAs using the edgeR package [25]. After differential expression analysis, we considered DE miRNAs with a range of |log2 (fold change)|≥ 1 and P valve < 0.05. We predicted the target genes of miRNA using miRanda [26], qTar, and RNAhybrid [27]. We performed the enrichment of GO and KEGG using DAVID Online with the target genes of miRNA [28].

LncRNA analysis and identification

We processed and mapped lncRNA data using SOAPnuke and HISAT2, following the same approach we used for miRNA data. We used samtools mapping to convert the sam file to a bam file [29]. The transcript assembly and quantification were accomplished by merging all bam files into a complete GTF file using the default parameters and merge capabilities of StringTie [30]. StringTie then calculated the lncRNA expression FPKM (Fragments per kilo-base of exon per million fragments mapped) of the assembled GTF file using the "-e" and "-B" parameters. The identification of lncRNA was based on existing studies [31]. Gffcompare in StringTie was used to get the transcript assembly [32]. We used gffcompare in StringTie to preserve the intergenic transcripts with class codes "u," "x," "i," "j," and "o," according to StringTie specifications. We deleted transcripts with a single exon and less than 200 bp in length. Based on the location information of the exons in the remaining transcripts, we extracted and fused the genome sequences of corresponding exons into a complete transcript sequence. We evaluated the protein-coding potential of the transcripts using CNCI [33], CPC2 [34], and PLEK [35]. Transcripts that encode proteins were removed to preserve those that do not encode proteins. We translated the previously screened transcripts into possible protein domainjjs using Transeq [36] and compared the corresponding proteins in the Pfam database [37]. We removed those with significant hits (E-value < 1e-5) using HMMER [38].

Prediction and annotation of lncRNA targets

LncRNAs function by acting on protein-coding genes through cis-acting elements and trans-acting factors. LncRNA trans target gene prediction is usually carried out when the sample size is large. We used Bedtools to predict the nearest target gene 100 kb upstream and downstream of the lncRNA [39]. We performed the enrichment of GO and KEGG using DAVID Online by the target genes of lncRNA.

Constructing the targeting relationship between miRNA and lncRNA

By using Cytoscape_v3.7.2 software, differentially expressed mirnas were combined with lncRN, and the targeting relationship between relevant lncrnas and mirnas was constructed.

ceRNA network construction

The Competing Endogenous RNA (ceRNA) regulatory network is a regulatory mechanism based on the mutual competition of RNA molecules. In this network, different types of RNA, such as mRNA, miRNA, and lncRNA, compete with each other for miRNA binding sites, thereby influencing each other's expression levels. This competitive relationship forms a complex regulatory network, believed to play a crucial role in the regulation of gene expression. We used Cytoscape_v3.7.2 software to construct a related cceRNA regulatory network by combining differentially expressed miRNAs and lncRNAs with mRNA predicted by the previous software.

Joint analysis of metabolites and genes.

Building and Analyzing Fully Connected Metabolite and mRNA Networks Using MetScape v3.1.3, a Plugin for Cytoscape Software [40, 41]. MetScape aids in constructing networks of metabolites and genes, tracing their interconnections, and visualizing the composite network.

Results

Analysis of differences in metabolic profiles of blood, urine, follicular fluid

The data of serum, follicular fluid, and urine samples obtained through LC–MS were integrated and analyzed using principal component analysis A (Fig. 1). In principal component analysis, the scoring maps of serum, urine, and follicular fluid samples showed significant changes in metabolites. On the basis of PCA analysis, the control group and blank control group of the three samples were subjected to multivariate data statistical analysis, and the PLS-DA model was used to further analyze and predict the inter group differences between different products (Fig. 2). The blank control group and control group were significantly separated, and the quality of PLS-DA was evaluated by two parameters R2 and Q2 (Fig. 3). Their explanatory power was good (R2), but their predictive power was poor (Q2), so the Variable Importance in Projection (VIP) cannot be used for screening differential metabolites.

Fig. 1
figure 1

PCA of metabolites in urine, blood and follicular fluid in estrus group and control group. A Urine B Blood C Follicular fluid

Full size image
Fig. 2
figure 2

PLSDA of metabolites in urine, blood and follicular fluid in estrus group and control group. A Urine B Blood C Follicular fluid

Full size image
Fig. 3
figure 3

PLSDA replacement test of metabolites in urine, blood and follicular fluid in estrus group and control group. A Urine B Blood C Follicular fluid

Full size image

Screening of differential metabolites

Differential metabolites were screened based on the following conditions: |log2 (fold change)|≥ 2 and p-value < 0.05. A volcano plot (Fig. 4) revealed 1239, 2, and 378 up-regulated metabolites in urine, blood, and follicular fluid, respectively, along with 2899, 30, and 971 down-regulated metabolites. However, these differential metabolites include isomers, which require enrichment into pathways for screening using the metaboanalyst website. Further screening is done based on the P-value (P < 0.05) obtained from one-way analysis of variance during pathway enrichment. The metabolites annotated into significantly enriched pathways are considered as the differential metabolites of interest. A total of 114, 13, and 91 different metabolites were screened in urine, blood, and follicular fluid, respectively, with most of them being amino acids, steroids, and organic acids (Table S1, Table S2, Table S3).

Fig. 4
figure 4

Volcanic plot in urine, blood and follicular fluid of estrus group and control group. A Urine B Blood C Follicular fluid

Full size image

Pathway enrichment analysis of differential metabolites

KEGG is a powerful tool for analyzing metabolism and conducting metabolic network research in vivo. Pathway enrichment analysis can help identify key biochemical metabolic pathways and signal transduction pathways associated with differential metabolites. In this study, KEGG pathway enrichment analysis was performed on differentially expressed metabolites in urine, blood, and follicular fluid from both the control group and the experimental group [42]. The results were visualized as bubble plots in Fig. 5. In the plots, the horizontal axis represents the number of differential metabolites in the corresponding metabolic pathway relative to the total number of identified metabolites in the pathway. A higher value indicates a greater enrichment of differential metabolites in the pathway. The color of the bubbles represents the P-value of the hypergeometric test, with smaller values (i.e., larger -log10(P-value)) indicating higher reliability and statistical significance. The size of the bubbles represents the number of differentiated metabolites in the corresponding pathway, with larger bubbles indicating a higher count of differentiated metabolites. In urine, the differential metabolites were mainly enriched in Arachidonic acid metabolism, Metabolism of xenobiotics by cytochrome P450, Tyrosine metabolism and other pathways. In follicular fluid, the differential metabolites were mainly enriched in Pentose and glucuronate interconversions, Arginine and proline metabolism, Purine metabolism and other pathways. In blood, the differential metabolites were mainly enriched in Pentose and glucuronate interconversions, D − Glutamine and D − glutamate metabolism, Alanine, aspartate and glutamate metabolism, and other pathways (Table S4).

Fig. 5
figure 5

KEGG enrichment in urine, blood and follicular fluid of estrus group and control group

Full size image

Overview of RNA sequencing data

We isolated follicular granulosa cells from yaks during estrus (n = 3, designated as group A) and from control yaks (n = 3, designated as group B) for RNA-seq analysis. The raw read counts for miRNA sequencing ranged from 28,915,662 to 43,636,363 reads, while the clean read counts ranged from 23,870,839 to 26,827,341 reads. The mapping read ratio varied from 70.09% to 85.08% across the six samples, and the Uniquely mapping ratio varied from 44.22% to 77.22%. For all samples, the Q20 score for clean tags was > 97.9% (RNA length is a crucial parameter for sRNA-seq, and we found that the majority of sRNAs were distributed between 20-23nt, with an average length range of 19.85–22.09 nt (Table S5). Regarding lncRNA sequencing, the clean read counts for the six samples ranged from 30,644,868 to 48,337,836 reads, and the GC content ranged from 51.52% to 55.69%. The Q20 score for clean tags was > 96.35% for all samples, the mapping read ratio varied from 91.42% to 92.72% across the six samples (Table S6).

Identification of differentially expressed miRNAs

After conducting a thorough comparison of the Mirbase database, we obtained counts and performed differential expression analysis using the R-DESeq2 v1.14.1 version, employing a threshold of |log2(Fold Change)|> 1 and a significance level of P valve < 0.05. By visually representing the results with a volcanic plot (Fig. 6A), we identified 44 miRNAs that were up-regulated and 78 miRNAs that were down-regulated.

Fig. 6
figure 6

A Volcanic Plot of miRNA B Venn diagrams predicting miRNA target genes using three softwares. C KEGG analysis of target genes for differentially expressed miRNA

Full size image

Analysis of target genes for differentially expressed miRNAs

The target genes of miRNAs were predicted using three different prediction software (miRanda, qTar, and RNAhybrid; Fig. 6B). Subsequently, GO and KEGG enrichment analysis were performed using the DAVID online website. The enrichment results of the KEGG pathway indicate that differentially expressed miRNA target genes are mainly enriched in Endocytosis, Type I diabetes mellitus, Cell adhesion molecules, Viral myocarditis and other pathways (Fig. 6C, Table S7). The GO functional annotation results show that in the biological process category, the main ones are cytoplasm, nucleus, cytosol and so on. In the cellular component, the main ones are immune response, integral component of lumenal side of endoplasmic reticulum membrane, antigen processing and presentation of endogenous peptide antigen via MHC class Ib and so on. In the molecular function, there is only a kinase binding (Fig. 7A, Table S8).

Fig. 7
figure 7

A GO analysis of target genes for differentially expressed miRNA. B GO analysis of target genes for differentially expressed lncRNA

Full size image

Identification of differentially expressed novel lncRNAs

Following the quantification of counts obtained from the featurecounts software [43], we conducted differential expression analysis using DESeq R package version 1, with a threshold of |log2 (Fold Change)|> 1 and a significance level of P valve < 0.05. The outcomes were visually represented using a volcano plot (Fig. 8A), which enabled us to identify 307 up-regulated and 333 down-regulated lncRNAs.

Fig. 8
figure 8

A Volcanic Plot of lncRNA B KEGG analysis of target genes for differentially expressed lncRNA

Full size image

Analysis of target genes for differentially expressed lncRNAs

The target gene of lncRNA was predicted using the prediction software Bedtools. Subsequently, GO and KEGG enrichment analysis was performed using the same method as miRNA target gene analysis. The enrichment results of the KEGG pathway indicate that differentially expressed lncRNA target genes are mainly enriched in Allograft rejection, Type I diabetes mellitus, Graft − versus − host disease, Cellular senescence, Cell adhesion molecules and so on (Fig. 8B, Table S9). The GO functional annotation results show that in the biological process category, the main ones are cytoplasm, nucleus, cytosol. In the cellular component, the main ones are ositive regulation of transcription from RNA polymerase II promoter, nucleolus, ubiquitin-dependent protein catabolic process. In the molecular function, the main ones are BMP signaling pathway, calmodulin binding, positive regulation of protein phosphorylation (Fig. 7B, Table S10).

Relationship between miRNA and lncRNA

To better elucidate the targeting relationship between lncRNA and miRNA, we chose the top 5 miRNAs exhibiting the most significant differential expression. These were then paired with differentially expressed lncRNAs to construct a focused targeting relationship map, encompassing 5 miRNAs and 14 lncRNAs. This strategic approach aims to provide a clear and concise visualization of the intricate interactions within this subset of differentially expressed molecules (Fig. 9).

Fig. 9
figure 9

Relationship between miRNA and lncRNA

Full size image

Construction of ceRNA network

lncRNAs can act as miRNA sponges, reducing their inhibitory effect on mRNA expression by sequestering miRNAs. A competing endogenous RNA (ceRNA) network involving lncRNAs and miRNAs was constructed based on differentially expressed lncRNA, miRNA, and predicted mRNA sequences. The ceRNA network comprised 44 nodes, including 10 lncRNAs, 4 miRNAs, and 30 mRNAs, and was visualized using Cytoscape software (Fig. 10). The ceRNA network revealed several mRNAs associated with estrus, such as GDF6 (Growth Differentiation Factor 6), TNFAIP8L2 (TNF Alpha Induced Protein 8 Like 2), COQ7 (Coenzyme Q7), SLC35G1 (Solute Carrier Family 35 Member G1), and AQP8 (Aquaporin 8). Additionally, several lncRNAs that interacted with miRNAs to exert ceRNA effects were identified, including MSTRG.960.1, MSTRG.960.7, MSTRG.4065.8, MSTRG.6961.7, MSTRG.8020.3, and others. These lncRNAs competitively inhibited miRNAs from binding to their corresponding mRNAs, thereby regulating mRNA expression.

Fig. 10
figure 10

CeRNA Regulatory Network in Yak granulosa cells

Full size image

Network analysis of differential metabolites and genes.

We used Metscape to construct a fully connected network of metabolites and mRNAs in urine, blood, and follicular fluid. Metabolites and mRNAs in the blood are mainly concentrated in the following pathways (Fig. 11): Glycine, serine, alanine and threonine metabolism, Glycolysis and Gluconeogenesis, Histidine metabolism, Purine metabolism, Urea cycle and metabolism of arginine, proline, glutamate, aspartate and asparagine, Vitamin B9 metabolism. Urine metabolites and mRNAs are mainly enriched to the following pathways (Fig. 12): Androgen and estrogen biosynthesis and metabolism, Arachidonic acid metabolism, Bile acid biosynthesis, C21-steroid hormone biosynthesis and metabolism, Glycine, serine alanine and threonine metabolism, Glycolysis and Gluconeogenesis and etc. on. Metabolites and mRNAs in follicular fluid are mainly enriched to the following pathway (Fig. 13): Aminosugars metabolism. Androgen and estrogen biosynthesis and metabolism, Bile acid biosynthesis, Biopterin metabolism, C21-steroid hormone biosynthesis and metabolism, Galactose metabolism, Purine metabolism, Pyrimidine metabolism, Tryptophan metabolism.

Fig. 11
figure 11

A network of metabolites and mRNAs in the blood

Full size image
Fig. 12
figure 12

A network of metabolites and mRNAs in the urine

Full size image
Fig. 13
figure 13

A network of metabolites and mRNAs in the follicular fluid

Full size image

Discussions

The yak is a unique breed of cattle that has been domesticated by highlanders living in the inhospitable high-altitude terrain of the central and inner Himalayas, serving as an indispensable breed for local herders in China. Despite its importance in the yak industry, the regulation of estrus in yaks is not well understood, with limited knowledge about the major endocrine and physiological mechanisms that may promote estrus in anestrus animals [44]. However, recent research on yak estrus is gaining momentum, with scholars conducting comprehensive analyses of the ceRNA network expression profiles of ovarian granulosa cells and cumulus cells in yaks, as well as transcriptome-wide analysis of N6-methyladenosine (m6A) to reveal the regulatory role of m6A in yak ovaries [45, 46]. In this study, we investigated the differentially expressed genes and differentially accumulated metabolites in the transcriptome and metabolome of estrous and unestrus yaks to explore the underlying mechanisms of estrus in yaks.

Changes in primary metabolites related to estrus

Studies have demonstrated that estrus in mammals occurs due to the specific effects of ovarian steroid hormones on behavioral centers in the brain [47]. Metabolic levels of amino acids [48], fatty acids [49], proteins [50], and other substances undergo changes in animals during estrus. Urine, blood, and follicular fluid are important body fluid components and serve as crucial environments for detecting metabolites. Taking urine as an example, mammalian urine contains nonvolatile substances such as uric acid, urea, and protein, as well as volatile and semi-volatile substances such as short-chain fatty acids and species-specific signal communication molecules called pheromones [51]. Urine may play a significant role as a pathway for communication during mammalian estrus. Estrus-specific sex pheromones are believed to carry out chemical signaling in the presence of urinary protein ligands [52]. In our research, most metabolites are amino acids, steroids, and organic acids. In urine, the differential metabolites were mainly enriched in Arachidonic acid metabolism, Metabolism of xenobiotics by cytochrome P450, Tyrosine metabolism and other pathways. In follicular fluid, the differential metabolites were mainly enriched in Pentose and glucuronate interconversions, Arginine and proline metabolism, Purine metabolism and other pathways. In blood, the differential metabolites were mainly enriched in Pentose and glucuronate interconversions, D − Glutamine and D − glutamate metabolism, Alanine, aspartate and glutamate metabolism, and other pathways. Many studies have verified the effects of different metabolites in follicular fluid on follicular growth through oocyte culture in vitro.SINCLAIR et al. reported that glycine and alanine were the most abundant amino acids found in follicular fluid [53]. MATOBA et al. conducted a metabolomics analysis of in vitro cultured oocytes under normal and stop-growth conditions and found that glycine, L-alanine, and L-glutamic acid were positively correlated with follicular growth [54]. Fatty acids and steroids also play essential roles in animal estrus. Studies have shown that estradiol, progesterone, and arachidonic acid regulate estrus in pigs through PKA and MAPK signaling pathways [55]. Improving the nutritional status of female mammals through dietary fatty acid content has been shown to enhance ovarian function and follicle development in goats [56]. These findings suggest that primary metabolic pathways, such as lipid metabolism and amino acid metabolism, play a crucial role in yak estrus.

A series of changes occur in follicle development during yak estrus. The follicular fluid in the follicular cavity is the microenvironment of the oocyte. Follicular fluid provides nutrients for follicle development and contains hormones, growth factors, proteins, and phospholipids produced in part by granulosa cells (GC) [57]. Similarly, metabolites in follicular fluid, such as advanced glycosylation end products, have been found to improve the developmental potential of oocytes. Steroid metabolism in follicles appears to be closely related to ovarian function and oocyte quality [58, 59]. In addition, studies have shown that vitamin b6 supplementation can significantly increase luteinizing hormone (LH) production in sows during estrus [60, 61]. This is similar to the types of differential metabolites and enrichment pathways in our findings.

Effect of miRNA on yak estrus

MiRNAs, which are a class of small non-coding RNAs, have been found to be closely associated with the growth and development of plants and animals in various studies [62]. Numerous studies have investigated the relationship between miRNAs and mammalian estrus. For instance, the interaction between miR-370-3p and target genes has been shown to play a crucial role in the reproductive regulation of small-tailed Han sheep, participating in the post-transcriptional processes of follicular and luteal stages in different genotypes of small-tailed Han sheep [63]. Yongfu La and their colleagues have also identified miRNAs that are associated with testicular development in male yaks [64], shedding light on the role of miRNAs in regulating testis development and improving reproductive performance in male yaks. Additionally, A E Zielak-Steciwko et al. have identified miRNAs that are related to bovine follicular development and found that miR-383 may serve as an important biological regulatory factor controlling the growth and proliferation of dominant follicles [65].

In our research, we investigated the miRNAs of follicular cells involved in yak estrus using the RNA-seq technique. Interestingly, we observed that the small RNA lengths were predominantly concentrated in the range of 21–23 nucleotides, which is similar to the typical size of miRNAs. Through a volcanic plot, we identified 44 up-regulated miRNAs and 78 down-regulated miRNAs, some of which are highly expressed during mammalian estrus, such as miR-129, miR-129-5P, and miR-183, suggesting their potential widespread involvement in yak estrus. The roles of these highly expressed miRNAs in estrus have been studied, and it has been found that a group of differentially expressed miRNAs identified by miR-182-5P may regulate follicle formation through apoptosis, proliferation, and differentiation of follicular cells, ovarian steroid production, and oxidative stress [66]. Furthermore, the miR-183 cluster has been found to be significantly differentially expressed in preovulation dominant follicles, potentially playing a role in the post-transcriptional regulation of bovine follicular development genes during the late follicular phase of the estrous cycle [67]. These highly expressed miRNAs are believed to play a crucial role in yak estrus and should be given more attention in future research. Most miRNAs bind to mRNAs and regulate gene expression, and by predicting the target genes of these miRNAs, we have identified some target genes that play important roles in estrus.

The enrichment of KEGG and GO is a common analysis method used to demonstrate the effects of miRNA target genes. GO enrichment analysis revealed that the target genes of most differentially expressed miRNAs were involved in seven molecular functions, eight cellular components, and 10 biological processes, including microtubule, ATPase activity, and protein phosphorylation, among others, suggesting strong regulatory effects of these miRNAs on cell metabolism, motility, intracellular and intercellular signaling, molecular function, and protein catalytic activity during ovarian development and hormone secretion during estrus. Additionally, the enrichment results of the KEGG pathway analysis indicated that differentially expressed miRNA target genes were mainly enriched in pathways such as Endocytosis, Type I diabetes mellitus, Allograft rejection, Viral myocarditis, and Cell adhesion molecules. After ovulation, progesterone triggers decidualization of endometrial stromal cells, resulting in cytoskeletal changes. This critical process regulates embryonic invasion and establishes the necessary cytokine and immunomodulator environment in the stromal tissue during invasion. Consequently, cell adhesion molecules play a pivotal role in ovarian maintenance [68].

Effect of lncRNA on yak estrus

lncRNA, a type of RNA molecule that exceeds 200 nucleotides and does not encode proteins, has been recognized for its diverse roles in various biological processes, making it a prominent area of research in the field of RNA in recent years. Recent studies have shed light on how lncRNAs can regulate gene expression at different levels, including epigenetic, transcriptional, and post-transcriptional regulation, and participate in important regulatory processes such as genome imprinting [69]. Growing evidence supports the presence and role of lncRNAs in reproductive tissues. For instance, Su and colleagues utilized high-throughput sequencing to identify differentially expressed lncRNAs during early pregnancy in Meishan and Yorkshire pigs, revealing that XLOC-2222497 and its target gene AKR1C1 can interact with progesterone in pig endometrium to regulate pregnancy maintenance [70]. Chen et al. identified 40 differentially expressed lncRNAs in the hypothalamus of sheep with different FecB genotypes during follicular and luteal phases [71]. In our own research, we identified 44 up-regulated and 78 down-regulated miRNAs, which were visualized using volcano plots. Studies have shown that lncRNA expression can regulate the expression of neighboring mRNAs [72]. Thus, we employed bedtools software to search for target genes located within 100 kb upstream and downstream of lncRNAs to predict the function of these lncRNAs. It is possible that lncRNAs act on adjacent or complementary target genes, playing a critical role in yak estrus. For example, MSTRG.73, a highly expressed lncRNA, targets the TIAM1 gene, which has been implicated in the invasive behavior of ovarian cancer, and down-regulation of TIAM1 is associated with increased aggressiveness of ovarian cancer cells [73]. Another target gene, MAP4K4, of MSTRG.2494, has been linked to the development of premature ovarian failure due to gene rearrangement [74]. Based on the function of these target genes, it is believed that these highly expressed lncRNAs may play a crucial role in yak estrus, and further investigation by researchers is anticipated in the future. Through the function of target genes, these highly expressed lncRNAs are believed to play a key role in yak estrus, which is expected to be further studied by scholars in the future.

The GO functional annotation results show that in the biological process category, the main ones are cytoplasm, nucleus, cytosol. In the cellular component, the main ones are ositive regulation of transcription from RNA polymerase II promoter, nucleolus, ubiquitin-dependent protein catabolic process. In the molecular function, the main ones are BMP signaling pathway, calmodulin binding, positive regulation of protein phosphorylation. Animals in estrus secrete hormones, and enzymes are needed for various physiological processes in estrus organisms, resulting in enhanced transcription, which is consistent with the biological processes of go enrichment analysis. The enrichment results of the KEGG pathway indicate that differentially expressed lncRNA target genes are mainly enriched in Allograft rejection, Type I diabetes mellitus, Graft − versus − host disease, Cellular senescence, Cell adhesion molecules and so on. The role of estrogen in the body relies on the transmission of signals between neurons and cells. Research indicates that estrogen can stimulate endocytosis of the postsynaptic membrane of neurons, resulting in increased density of small particles in the membrane and subsequently reducing the number of small particles. This process seems to be the physiological mechanism by which neurons in females change during the estrus cycle [75].

Recent literature reports suggest that the mutual regulatory network involving miRNA, lncRNA, and mRNA may also participate in the estrus process of yaks [45]. In the ceRNA network of yak follicular cells during estrus, a total of 10 lncRNAs, 4 miRNAs, and 30 mRNAs were identified. miR-129-5p is a representative miRNA [66], and lncRNA of MTSRG.26652 serves as a representative of the lncRNA group [76]. MRNA of GDF6 is also involved in animal estrus [77]. miR-129 is a key miRNA, and the nodes in the ceRNA network directly or indirectly participate in yak estrus. The role of the miR-129 family in animal ovarian estrus has been a research hotspot in recent years, and mir-129-5p, as a member of this family, has been shown to affect follicular development in sheep under the stimulation of FSH [66].

Association analysis of transcriptomics and metabolome

Metabolomics and transcriptomics are potent omics technologies that yield comprehensive metabolite and transcript profiles [78]. The study of transcriptomics and metabolomics is in full swing in animal husbandryrespectively. In a recent investigation of Qinchuan cattle, a holistic examination involving blood transcriptomics and metabolomics was undertaken to unveil the molecular mechanisms regulating backfat thickness. The study revealed 1106 differentially expressed genes and 86 differentially expressed metabolites in cattle with varying backfat thickness. Noteworthy among these were the functional genes SMPD3 and CERS1, along with the metabolite sphingosine 1-phosphate, identified as pivotal factors influencing phenotype variations. These results underscore the efficacy of integrating metabolomics and transcriptomics in livestock research, offering valuable insights for enhancing beef performance through informed breeding strategies [79]. In a recent study on pigs, researchers investigated the impact of miR-149-5p on intramuscular fat deposition using a comprehensive approach involving metabolomics and transcriptomics. The study revealed that miR-149-5p promoted the proliferation of porcine intramuscular preadipocytes while concurrently reducing their differentiation ability. Metabolomics analysis further demonstrated significant alterations in lipid, organic acid, and organic oxygen metabolites upon the overexpression of miR-149-5p in porcine intramuscular preadipocytes [80]. In a separate inquiry, an integrated analysis of liver metabolomics and transcriptomics was conducted to unveil oxidative stress in piglets affected by intrauterine growth restriction. The study highlighted a range of metabolic abnormalities in intrauterine growth-restricted piglets, including mitochondrial dysfunction, imbalances in fatty acid composition, disruptions in one-carbon unit supply sources, and abnormal galactose conversion. These irregularities are suggested to contribute to oxidative stress in the liver [81].

In our study, we employed MetScape software to construct and analyze metabolic and mRNA networks in urine, blood, and follicular fluid. In the blood (Fig. 11), the two most enriched pathways for metabolites and mRNA were the urea cycle and the metabolism of arginine, proline, glutamic acid, aspartic acid, and asparagine, along with purine metabolism. The role of purine metabolic pathways in reproduction remains unclear [82]. In the present investigation, purines and pyrimidines served as precursors for DNA and RNA biosynthesis [83], playing a crucial role as genetic material in normal cellular functions [84, 85]. During the estrous cycle of female animals, we hypothesized that cells might have needed to synthesize new DNA and RNA to support the development of ova and the division of reproductive cells, potentially achieved through the activation of purine metabolism. Purines could undergo a series of reactions to be broken down into uric acid, which was then excreted through urine. In our differential metabolite analysis, the levels of purines in the estrus group decreased compared to the control group, aligning with our speculative hypothesis. In the urine (Fig. 12), the gene most prominently associated in our study is the 17β-hydroxysteroid dehydrogenase family. The 17β-hydroxysteroid dehydrogenase family (17β-HSD) constitutes a group of enzymes that catalyze the conversion of 17-ketosteroid steroids into biologically more potent estrogens and androgens. For instance, they facilitate the conversion of estrone (E1) to the more active estradiol (E2) and the transformation of androstenedione into testosterone [86]. The 17β-hydroxysteroid dehydrogenase 12 (HSD17B12) has been demonstrated to play a role in the regulation of the prostaglandin synthesis pathway, ovarian function, and fertility. The HSD17B12 enzyme exhibits widespread expression in both human and murine ovaries. In mice, heterozygous deletion of the HSD17B12 gene results in subfertility, indicating its significant role in ovarian function. Heterozygous female mice with a deficiency in HSD17B12 experience a significantly prolonged estrous cycle, reduced fertility compared to wild-type females, and a noticeable decrease in litter size. Additionally, abnormal formation of the mitotic spindle in immature follicles suggests a defect in cell meiotic division arrest in the ovaries of HSD17B12 heterozygotes [87]. Therefore, the transcription levels of the family genes would affect the mating behavior of animals. These family genes are linked to estradiol 17β-dehydrogenase, and these specific substances play a crucial role in the estrous behavior of female animals.Estradiol 17β-dehydrogenase is an enzyme that catalyzes the chemical reaction that converts estradiol-17β and NAD(P)+ to ketone and NAD(P)H+H+ [88]. During the estrous cycle of animals, estrogens (such as estradiol) play an important role, affecting sexual behavior and the development of reproductive organs [89, 90]. Therefore, the correlation between metabolites and mRNA in our urine may play a role in the mating process of animals, offering theoretical support that can be validated by future scholars.

In follicular fluid (Fig. 13), the mRNA associated with metabolites includes NOS3 and so on. The NOS3 gene encodes an enzyme called endothelial nitric oxide synthase (eNOS), which synthesizes nitric oxide (NO) in the organism. Nitric oxide is a biological signaling molecule involved in various physiological processes, including neural transmission, antimicrobial activity, anti-tumor activity, and reproductive processes [91]. Research suggests that different subtypes of NOS are found in the oviducts of mice and humans at various stages of the estrous and menstrual cycles by assessing NOS levels. Increased expression of iNOS has been observed in the oviducts of women with ectopic pregnancies [92]. The metabolite associated with NOS3 is Nitric-oxide synthase, a molecule that plays a particularly important role in the estrous cycle of female animals.

In summary, we have investigated the correlation between differential metabolites in urine, blood, and follicular fluid and the potentially targeted differential mRNA in the transcriptome of oocyte granulosa cells. This study provides a reference for future research on estrus-related studies in yaks.

Conclusion

In conclusion, this study has identified differences in metabolites among three body fluids of female yaks after estrus treatment, shedding light on changes in body fluid metabolites during yak estrus and providing fundamental data for understanding the mechanisms of yak estrus in vivo. Furthermore, differentially expressed miRNAs and lncRNAs in follicular cells of female yaks have been identified, and a ceRNA regulatory network involving lncRNA-miRNA-mRNA has been constructed. Key genes associated with yak estrus have been discovered in this network, which may serve as potential biomarkers for future studies on yak estrus. These key genes are expected to provide valuable insights for further research in this field.

Availability of data and materials

Sequences are available from GenBank with the Bioproject accession number PRJNA958448 and PRJNA958176.

Change history

References

  1. Li Z, Jiang M. Metabolomic profiles in yak mammary gland tissue during the lactation cycle. PLoS One. 2019;14(7):e0219220.

    CAS  PubMed  PubMed Central  Google Scholar 

  2. Jiang M, Lee JN, Bionaz M, Deng XY, Wang Y. Evaluation of suitable internal control genes for rt-qpcr in yak mammary tissue during the lactation cycle. PLoS One. 2016;11(1):e0147705.

    PubMed  PubMed Central  Google Scholar 

  3. Prakash BS, Sarkar M, Mondal M. An update on reproduction in yak and mithun. Reprod Domest Anim. 2008;43(Suppl 2):217–23.

    PubMed  Google Scholar 

  4. Long RJ, Zhang DG, Wang X, Hu ZZ, Dong SK. Effect of strategic feed supplementation on productive and reproductive performance in yak cows. Prev Vet Med. 1999;38(2–3):195–206.

    CAS  PubMed  Google Scholar 

  5. Zi XD. Reproduction in female yaks (Bos grunniens) and opportunities for improvement. Theriogenology. 2003;59(5–6):1303–12.

    PubMed  Google Scholar 

  6. Robin N, Laforest J, Lussier J, Guilbault L. Induction of estrus with intramuscular injections of GnRH or PMSG in lactating goats (Capra hircus ) primed with a progestagen during seasonal anestrus. Theriogenology. 1994;42(1):107–16.

    CAS  PubMed  Google Scholar 

  7. Wickramasuriya N, Hawkins R, Atwood C, Butler T. The roles of GnRH in the human central nervous system. Horm Behav. 2022;145:105230.

    CAS  PubMed  PubMed Central  Google Scholar 

  8. Zhang J, Wang C, Li X, Zhang Y, Xing F. Expression and functional analysis of GnRH at the onset of puberty in sheep. Arch Anim Breed. 2022;65(3):249–57.

    PubMed  PubMed Central  Google Scholar 

  9. Zanetti BF, Braga DPAF, Setti AS, Iaconelli A Jr, Borges E Jr. Effect of GnRH analogues for pituitary suppression on oocyte morphology in repeated ovarian stimulation cycles. JBRA Assist Reprod. 2020;24(1):24–9.

    PubMed  PubMed Central  Google Scholar 

  10. Hao D, Bai J, Du J, Wu X, Thomsen B, Gao H, Su G, Wang X. Overview of metabolomic analysis and the integration with multi-omics for economic traits in cattle. Metabolites. 2021;11(11):753.

    CAS  PubMed  PubMed Central  Google Scholar 

  11. Hasin Y, Seldin M, Lusis A. Multi-omics approaches to disease. Genome Biol. 2017;18(1):83.

    PubMed  PubMed Central  Google Scholar 

  12. Liu L, Wu P, Chen F, Zhou J, Guo A, Shi K, Zhang Q. Multi-omics analyses reveal that the gut microbiome and its metabolites promote milk fat synthesis in Zhongdian yak cows. PeerJ. 2022;2(10):e14444.

    Google Scholar 

  13. Zhang L, Wang Z, Zhou P, Fu L, Zhang L, Xu C, Loor JJ, Zhang T, Chen Y, Zhou Z, Dong X. Vitamin E supplementation improves post-transportation systemic antioxidant capacity in yak. PLoS One. 2022;17(12):e0278660.

    CAS  PubMed  PubMed Central  Google Scholar 

  14. Wang Q, Dong K, Wu Y, An F, Luo Z, Huang Q, Wei S. Exploring the formation mechanism of off-flavor of irradiated yak meat based on metabolomics. Food Chem X. 2022;4(16):100494.

    Google Scholar 

  15. Zhou J, Yue S, Du J, Xue B, Wang L, Peng Q, Zou H, Hu R, Jiang Y, Wang Z, Xue B. Integration of transcriptomic and metabolomic analysis of the mechanism of dietary N-carbamoylglutamate in promoting follicle development in yaks. Front Vet Sci. 2022;29(9):946893.

    Google Scholar 

  16. Vanderhyden T. Eppig: mouse oocytes promote proliferation of granulosa cells from preantral and antral follicles in vitro. Biol Reprod. 1992;46(6):1196–204.

    CAS  PubMed  Google Scholar 

  17. Williams GL, Stanko RL. Pregnancy rates to fixed-time AI in Bos indicus-influenced beef cows using PGF2α with (Bee Synch I) or without (Bee Synch II) GnRH at the onset of the 5-day CO-Synch + CIDR protocol. Theriogenology. 2020;15(142):229–35.

    Google Scholar 

  18. Chen D, Qi X, Liu J, et al. Improvement of the isolation and culture method for human ovarian granulosa cells(in Chinese). Chin J Tissue Eng Res. 2015;46:7456–60.

    Google Scholar 

  19. Dunn WB, Broadhurst D, Begley P, Zelena E, Francis-McIntyre S, Anderson N, Brown M, Knowles JD, Halsall A, Haselden JN, Nicholls AW, Wilson ID, Kell DB, Goodacre R, C, H. S. M. H. Procedures for large-scale metabolic profiling of serum and plasma using gas chromatography and liquid chromatography coupled to mass spectrometry. Nat Protoc. 2011;6(7):1060–83.

  20. Chen Y, Chen Y, Shi C, Huang Z, Zhang Y, Li S, Li Y, Ye J, Yu C, Li Z, Zhang X, Wang J, Yang H, Fang L, Chen Q. SOAPnuke: a MapReduce acceleration-supported software for integrated quality control and preprocessing of high-throughput sequencing data. Gigascience. 2018;7(1):1–6.

    PubMed  PubMed Central  Google Scholar 

  21. Kim D, Paggi JM, Park C, Bennett C, Salzberg SL. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat Biotechnol. 2019;37(8):907–15.

    CAS  PubMed  PubMed Central  Google Scholar 

  22. Kozomara A, Birgaoanu M, Griffiths-Jones S. miRBase: from microRNA sequences to function. Nucleic Acids Res. 2019;47(D1):D155–62.

    CAS  PubMed  Google Scholar 

  23. Griffiths-Jones S, Grocock RJ, van Dongen S, Bateman A, Enright AJ. miRBase: microRNA sequences, targets and gene nomenclature. Nucleic Acids Res. 2006;34(Database issue):D140-4.

    CAS  PubMed  Google Scholar 

  24. Feng S, Ma J, Long K, Zhang J, Qiu W, Li Y, Jin L, Wang X, Jiang A, Liu L, Xiao W, Li X, Tang Q, Li M. Comparative microRNA transcriptomes in domestic goats reveal acclimatization to high altitude. Front Genet. 2020;31(11):809.

    Google Scholar 

  25. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.

    CAS  PubMed  Google Scholar 

  26. John B, Enright AJ, Aravin A, Tuschl T, Sander C, Marks DS. Human MicroRNA targets. PLoS Biol. 2004;2(11):e363.

    PubMed  PubMed Central  Google Scholar 

  27. Krüger J, Rehmsmeier M. RNAhybrid: microRNA target prediction easy, fast and flexible. Nucleic Acids Res. 2006;34(Web Server issue):W451-4.

    PubMed  PubMed Central  Google Scholar 

  28. Sherman BT, Hao M, Qiu J, Jiao X, Baseler MW, Lane HC, Imamichi T, Chang W. DAVID: a web server for functional enrichment analysis and functional annotation of gene lists. Nucleic Acids Res. 2022;50(W1):W216–21.

    CAS  PubMed  PubMed Central  Google Scholar 

  29. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, 1000 Genome Project Data Processing Subgroup. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.

    PubMed  PubMed Central  Google Scholar 

  30. Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. 2016;11(9):1650–67.

    CAS  PubMed  PubMed Central  Google Scholar 

  31. Chen L, Shi G, Chen G, Li J, Li M, Zou C, Fang C, Li C. Transcriptome analysis suggests the roles of long intergenic non-coding RNAs in the growth performance of weaned piglets. Front Genet. 2019;18(10):196.

    Google Scholar 

  32. Pertea G, Pertea M. GFF Utilities: GffRead and GffCompare. F1000Res. 2020;9:ISCB Comm J-304.

    PubMed  PubMed Central  Google Scholar 

  33. Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, Liu Y, Chen R, Zhao Y. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013;41(17):e166.

    CAS  PubMed  PubMed Central  Google Scholar 

  34. Kang YJ, Yang DC, Kong L, Hou M, Meng YQ, Wei L, Gao G. CPC2: a fast and accurate coding potential calculator based on sequence intrinsic features. Nucleic Acids Res. 2017;45(W1):W12–6.

    CAS  PubMed  PubMed Central  Google Scholar 

  35. Li A, Zhang J, Zhou Z. PLEK: a tool for predicting long non-coding RNAs and messenger RNAs based on an improved k-mer scheme. BMC Bioinformatics. 2014;15(1):311.

    PubMed  PubMed Central  Google Scholar 

  36. Rice P, Longden I, Bleasby A. EMBOSS: the European Molecular Biology Open Software Suite. Trends Genet. 2000;16(6):276–7.

    CAS  PubMed  Google Scholar 

  37. Mistry J, Chuguransky S, Williams L, Qureshi M, Salazar GA, Sonnhammer ELL, Tosatto SCE, Paladin L, Raj S, Richardson LJ, Finn RD, Bateman A. Pfam: the protein families database in 2021. Nucleic Acids Res. 2021;49(D1):D412–9.

    CAS  PubMed  Google Scholar 

  38. Potter SC, Luciani A, Eddy SR, Park Y, Lopez R, Finn RD. HMMER web server: 2018 update. Nucleic Acids Res. 2018;46(W1):W200–4.

    CAS  PubMed  PubMed Central  Google Scholar 

  39. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26(6):841–2.

    CAS  PubMed  PubMed Central  Google Scholar 

  40. Karnovsky A, Weymouth T, Hull T, Tarcea VG, Scardoni G, Laudanna C, Sartor MA, Stringer KA, Jagadish HV, Burant C, Athey B, Omenn GS. Metscape 2 bioinformatics tool for the analysis and visualization of metabolomics and gene expression data. Bioinformatics. 2012;28(3):373–80.

    CAS  PubMed  Google Scholar 

  41. Basu S, Duren W, Evans CR, Burant CF, Michailidis G, Karnovsky A. Sparse network modeling and metscape-based visualization methods for the analysis of large-scale metabolomics data. Bioinformatics. 2017;33(10):1545–53.

    CAS  PubMed  PubMed Central  Google Scholar 

  42. Ogata H, Goto S, Sato K, Fujibuchi W, Bono H, Kanehisa M. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 1999;27(1):29–34. https://doi.org/10.1093/nar/27.1.29.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30(7):923–30.

    CAS  PubMed  Google Scholar 

  44. Xie J, Kalwar Q, Yan P, Guo X. Expression and characterization of the serum proteome from yak induced into estrus by improved nutrition. Anim Biotechnol. 2022;33(5):930–40.

    CAS  PubMed  Google Scholar 

  45. Zhao L, Pan Y, Wang M, Wang J, Wang Y, Han X, Wang J, Zhang T, Zhao T, He H, Cui Y, Yu S. Integrated analysis of the expression profiles of the lncRNA-miRNA-mRNA ceRNA network in granulosa and cumulus cells from yak ovaries. BMC Genomics. 2022;23(1):633.

    CAS  PubMed  PubMed Central  Google Scholar 

  46. Guo S, Wang X, Cao M, Wu X, Xiong L, Bao P, Chu M, Liang C, Yan P, Pei J, Guo X. The transcriptome-wide N6-methyladenosine (m6A) map profiling reveals the regulatory role of m6A in the yak ovary. BMC Genomics. 2022;23(1):358.

    CAS  PubMed  PubMed Central  Google Scholar 

  47. Roelofs J, López-Gatius F, Hunter RH, van Eerdenburg FJ, Hanzen Ch. When is a cow in estrus? Clinical and practical aspects. Theriogenology. 2010;74(3):327–44.

    CAS  PubMed  Google Scholar 

  48. Jones DB, Stahly TS. Impact of amino acid nutrition during lactation on luteinizing hormone secretion and return to estrus in primiparous sows. J Anim Sci. 1999;77(6):1523–31.

    CAS  PubMed  Google Scholar 

  49. Oliveira JD, Fonseca JF, Souza-Fabjan JM, Esteves LV, Feres LF, Rodrigues CA, Torres Filho RA, de Oliveira J, Brandão FZ. Protected fatty acid supplementation during estrus synchronization treatment on reproductive parameters of dairy goats. Anim Sci J. 2017;88(2):254–8.

    CAS  PubMed  Google Scholar 

  50. Zhao C, Shu S, Bai Y, Wang D, Xia C, Xu C. Plasma protein comparison between dairy cows with inactive ovaries and estrus. Sci Rep. 2019;9(1):13709.

    PubMed  PubMed Central  Google Scholar 

  51. Soso SB, Koziel JA. Characterizing the scent and chemical composition of Panthera leo marking fluid using solid-phase microextraction and multidimensional gas chromatography-mass spectrometry-olfactometry. Sci Rep. 2017;7(1):5137.

    PubMed  PubMed Central  Google Scholar 

  52. Beynon RJ, Hurst JL. Urinary proteins and the modulation of chemical scents in mice and rats. Peptides. 2004;25(9):1553–63.

    CAS  PubMed  Google Scholar 

  53. Sinclair KD, Lunn LA, Kwong WY, Wonnacott K, Linforth RS, Craigon J. Amino acid and fatty acid composition of follicular fluid as predictors of in-vitro embryo development. Reprod Biomed Online. 2008;16(6):859–68.

    CAS  PubMed  Google Scholar 

  54. Matoba S, Bender K, Fahey AG, Mamo S, Brennan L, Lonergan P, Fair T. Predictive value of bovine follicular components as markers of oocyte developmental potential. Reprod Fertil Dev. 2014;26(2):337–45.

    CAS  PubMed  Google Scholar 

  55. Skowronska A, MÅ‚otkowska P, Wojciechowicz B, Okrasa S, Nielsen S, Skowronski MT. Progesterone, estradiol, arachidonic acid, oxytocin, forskolin and cAMP influence on aquaporin 1 and 5 expression in porcine uterine explants during the mid-luteal phase of the estrous cycle and luteolysis: an in vitro study. Reprod Biol Endocrinol. 2015;18(13):7.

    Google Scholar 

  56. Nugroho P, Wiryawan KG, Astuti DA, Manalu W. Stimulation of follicle growth and development during estrus in Ettawa Grade does fed a flushing supplement of different polyunsaturated fatty acids. Vet World. 2021;14(1):11–22.

    CAS  PubMed  PubMed Central  Google Scholar 

  57. Zhang X, Wang T, Song J, Deng J, Sun Z. Study on follicular fluid metabolomics components at different ages based on lipid metabolism. Reprod Biol Endocrinol. 2020;18(1):42.

    CAS  PubMed  PubMed Central  Google Scholar 

  58. Huang Y, Tu M, Qian Y, Ma J, Chen L, Liu Y, Wu Y, Chen K, Liu J, Ying Y, Chen Y, Ye Y, Xing L, Zhang F, Hu Y, Zhang R, Ruan YC, Zhang D. Age-Dependent Metabolomic Profile of the Follicular Fluids From Women Undergoing Assisted Reproductive Technology Treatment. Front Endocrinol (Lausanne). 2022;16(13):818888.

    Google Scholar 

  59. Jinno M, Nagai R, Takeuchi M, Watanabe A, Teruya K, Sugawa H, et al. Trapa bispinosa Roxb. Extract lowers advanced glycation end-products and increases live births in older patients with assisted reproductive technology: a randomized controlled trial. Reprod Biol Endocrinol. 2021;19(1):149.

    CAS  PubMed  PubMed Central  Google Scholar 

  60. Dalto DB, Roy M, Audet I, Palin M-F, Guay F, Lapointe J, et al. Interaction between vitamin B6 and source of selenium on the response of the selenium-dependent glutathione peroxidase system to oxidative stress induced by oestrus in pubertal pig. J Trace Elem Med Biol. 2015;32:21–9.

    CAS  PubMed  Google Scholar 

  61. Zhang J, Liu M, Ke S, Huang X, Fang S, He M, Fu H, Chen C, Huang L. Gut and vagina microbiota associated with estrus return of weaning sows and its correlation with the changes in serum metabolites. Front Microbiol. 2021;19(12):690091.

    Google Scholar 

  62. Ullah Y, Li C, Li X, Ni W, Yao R, Xu Y, Quan R, Li H, Zhang M, Liu L, Hu R, Guo T, Li Y, Wang X, Hu S. Identification and profiling of pituitary microRNAs of sheep during anestrus and estrus stages. Animals. 2020;10(3):402.

    PubMed  PubMed Central  Google Scholar 

  63. Chen Y, Liu Y, Chu M. miRNA-mRNA analysis of sheep adrenal glands reveals the network regulating reproduction. BMC Genom Data. 2022;23(1):44.

    CAS  PubMed  PubMed Central  Google Scholar 

  64. La Y, Ma X, Bao P, Chu M, Guo X, Liang C, Yan P. Identification and profiling of microRNAs during yak’s testicular development. BMC Vet Res. 2023;19(1):53.

    CAS  PubMed  PubMed Central  Google Scholar 

  65. Zielak-Steciwko AE, Browne JA, McGettigan PA, Gajewska M, Dzięcioł M, Szulc T, Evans AC. Expression of microRNAs and their target genes and pathways associated with ovarian follicle development in cattle. Physiol Genomics. 2014;46(19):735–45.

    CAS  PubMed  Google Scholar 

  66. Yuan H, Lu J, Xiao SY, Han XY, Song XT, Qi MY, Liu GS, Yang CX, Yao YC. miRNA expression analysis of the sheep follicle during the prerecruitment, dominant, and mature stages of development under FSH stimulation. Theriogenology. 2022;15(181):161–9.

    Google Scholar 

  67. Gebremedhn S, Salilew-Wondim D, Ahmad I, Sahadevan S, Hossain MM, Hoelker M, Rings F, Neuhoff C, Tholen E, Looft C, Schellander K, Tesfaye D. MicroRNA expression profile in bovine granulosa cells of preovulatory dominant and subordinate follicles during the late follicular phase of the estrous cycle. PLoS One. 2015;10(5):e0125912.

    PubMed  PubMed Central  Google Scholar 

  68. Baracat MC, Serafini PC, Simões Rdos S, Maciel GA, Soares JM Jr, Baracat EC. Systematic review of cell adhesion molecules and estrogen receptor expression in the endometrium of patients with polycystic ovary syndrome. Int J Gynaecol Obstet. 2015;129(1):1–4.

    CAS  PubMed  Google Scholar 

  69. Herman AB, Tsitsipatis D, Gorospe M. Integrated lncRNA function upon genomic and epigenomic regulation. Mol Cell. 2022;82(12):2252–66.

    CAS  PubMed  PubMed Central  Google Scholar 

  70. Su T, Yu H, Luo G, Wang M, Zhou C, Zhang L, Hou B, Zhang C, Liu M, Xu D. The interaction of lncRNA XLOC-2222497, AKR1C1, and progesterone in porcine endometrium and pregnancy. Int J Mol Sci. 2020;21(9):3232.

    CAS  PubMed  PubMed Central  Google Scholar 

  71. Chen S, Guo X, He X, Di R, Zhang X, Zhang J, Wang X, Chu M. 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;20(9):633747.

    Google Scholar 

  72. Guttman M, Amit I, Garber M, French C, Lin MF, Feldser D, Huarte M, Zuk O, Carey BW, Cassady JP, Cabili MN, Jaenisch R, Mikkelsen TS, Jacks T, Hacohen N, Bernstein BE, Kellis M, Regev A, Rinn JL, Lander ES. Chromatin signature reveals over a thousand highly conserved large non-coding RNAs in mammals. Nature. 2009;458(7235):223–7.

    CAS  PubMed  PubMed Central  Google Scholar 

  73. Li J, Liang S, Jin H, Xu C, Ma D, Lu X. Tiam1, negatively regulated by miR-22, miR-183 and miR-31, is involved in migration, invasion and viability of ovarian cancer cells. Oncol Rep. 2012;27(6):1835–42.

    CAS  PubMed  Google Scholar 

  74. Ledig S, Röpke A, Wieacker P. Copy number variants in premature ovarian failure and ovarian dysgenesis. Sex Dev. 2010;4(4–5):225–32.

    CAS  PubMed  Google Scholar 

  75. Leedom L, Lewis C, Garcia-Segura LM, Naftolin F. Regulation of arcuate nucleus synaptology by estrogen. Ann N Y Acad Sci. 1994;14(743):61–71.

    Google Scholar 

  76. Liu Z, Dai S, Bones J, Ray S, Cha S, Karger BL, Li JJ, Wilson L, Hinckle G, Rossomando A. A quantitative proteomic analysis of cellular responses to high glucose media in Chinese hamster ovary cells. Biotechnol Prog. 2015;31(4):1026–38.

    CAS  PubMed  Google Scholar 

  77. Chen Q, Wang Y, Liu Z, Guo X, Sun Y, Kang L, Jiang Y. Transcriptomic and proteomic analyses of ovarian follicles reveal the role of VLDLR in chicken follicle selection. BMC Genomics. 2020;21(1):486.

    CAS  PubMed  PubMed Central  Google Scholar 

  78. Goldansaz SA, Guo AC, Sajed T, Steele MA, Plastow GS, Wishart DS. Livestock metabolomics and the livestock metabolome: a systematic review. PLoS One. 2017;12(5):e0177675. https://doi.org/10.1371/journal.pone.0177675. PMID:28531195;PMCID:PMC5439675.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  79. Yu H, Raza SHA, Pan Y, Cheng G, Mei C, Zan L. Integrative analysis of blood transcriptomics and metabolomics reveals molecular regulation of backfat thickness in Qinchuan Cattle. Animals (Basel). 2023;13(6):1060.

    PubMed  Google Scholar 

  80. Liu Y, Wei Y, Dou Y, Li C, Song C, Zhang Z, Qi K, Li X, Qiao R, Wang K, Li X, Yang F, Han X. Effect of miR-149-5p on intramuscular fat deposition in pigs based on metabolomics and transcriptomics. BMC Genomics. 2023;24(1):293.

    CAS  PubMed  PubMed Central  Google Scholar 

  81. Gao H, Chen X, Zhao J, Xue Z, Zhang L, Zhao F, Wang B, Wang L. Integrative analysis of liver metabolomics and transcriptomics reveals oxidative stress in piglets with intrauterine growth restriction. Biology (Basel). 2022;11(10):1430.

    CAS  PubMed  Google Scholar 

  82. Chen M, Zhang B, Cai S, Zeng X, Ye Q, Mao X, Zhang S, Zeng X, Ye C, Qiao S. Metabolic disorder of amino acids, fatty acids and purines reflects the decreases in oocyte quality and potential in sows. J Proteomics. 2019;30(200):134–43.

    CAS  Google Scholar 

  83. Lee CC, Yang YC, Goodman SD, Chen S, Huang TY, Cheng WC, Lin LI, Fang WH. Deoxyinosine repair in nuclear extracts of human cells. Cell Biosci. 2015;8(5):52.

    Google Scholar 

  84. Unnikrishnan A, Freeman WM, Jackson J, Wren JD, Porter H, Richardson A. The role of DNA methylation in epigenetics of aging. Pharmacol Ther. 2019;195:172–85.

    CAS  PubMed  Google Scholar 

  85. Mercer TR, Dinger ME, Mattick JS. Long non-coding RNAs: insights into functions. Nat Rev Genet. 2009;10(3):155–9.

    CAS  PubMed  Google Scholar 

  86. Wang P, Zheng D, Peng W, Wang Y, Wang X, Xiong W, Liang R. Characterization of 17β-hydroxysteroid dehydrogenase and regulators involved in estrogen degradation in Pseudomonas putida SJTE-1. Appl Microbiol Biotechnol. 2019;103(5):2413–25.

    CAS  PubMed  Google Scholar 

  87. Kemiläinen H, Adam M, Mäki-Jouppila J, Damdimopoulou P, Damdimopoulos AE, Kere J, Hovatta O, Laajala TD, Aittokallio T, Adamski J, Ryberg H, Ohlsson C, Strauss L, Poutanen M. The hydroxysteroid (17β) dehydrogenase family gene HSD17B12 is involved in the prostaglandin synthesis pathway, the ovarian function, and regulation of fertility. Endocrinology. 2016;157(10):3719–30.

    PubMed  Google Scholar 

  88. Hao P, Pan H, Lv Z, Zhang J, Wang L, Zhu Y, Basang W, Gao Y. Characterization of 17β-estradiol-degrading enzyme from Microbacterium sp. MZT7 and its function on E2 biodegradation in wastewater. Microb Cell Fact. 2023;22(1):116.

    CAS  PubMed  PubMed Central  Google Scholar 

  89. Hilborn E, Stål O, Jansson A. Estrogen and androgen-converting enzymes 17β-hydroxysteroid dehydrogenase and their involvement in cancer: with a special focus on 17β-hydroxysteroid dehydrogenase type 1, 2, and breast cancer. Oncotarget. 2017;8(18):30552–62.

    PubMed  PubMed Central  Google Scholar 

  90. Nazari E, Suja F. Effects of 17β-estradiol (E2) on aqueous organisms and its treatment problem: a review. Rev Environ Health. 2016;31(4):465–91.

    CAS  PubMed  Google Scholar 

  91. Luo Y, Zhu Y, Basang W, Wang X, Li C, Zhou X. Roles of nitric oxide in the regulation of reproduction: a review. Front Endocrinol (Lausanne). 2021;19(12):752410.

    Google Scholar 

  92. Hu J, Ma S, Zou S, Li X, Cui P, Weijdegård B, Wu G, Shao R, Billig H, Feng Y. The regulation of nitric oxide synthase isoform expression in mouse and human fallopian tubes: potential insights for ectopic pregnancy. Int J Mol Sci. 2014;16(1):49–67.

    PubMed  PubMed Central  Google Scholar 

Download references

Funding

This research has been supported by the Key Technology Support Program of Qinghai Province (2021-ZJ-736).

Author information

Author notes

  1. Huangqing Zhao and Yongzhen Huang are co-first authors.

Authors and Affiliations

  1. Qinghai University, Xining, China

    Huangqing Zhao, Shi Shu, Guowen Wang, Changqi Fu, Rong Huang, Jun Zhang & Wei Peng

  2. Northwest A&F University, Yangling, Shaanxi, China

    Huangqing Zhao, Yongzhen Huang, Chuzhao Lei, Lei Du & Jiahao Zhao

  3. China Agricultural University, Beijing, China

    Huawei Su & Yang He

Authors
  1. Huangqing Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Yongzhen Huang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Shi Shu
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Guowen Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Changqi Fu
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Rong Huang
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Jun Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Huawei Su
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Yang He
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Chuzhao Lei
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Lei Du
    View author publications

    You can also search for this author in PubMed Google Scholar

  12. Jiahao Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

  13. Wei Peng
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

HZ and YH were responsible for analyzing data and writing, SS, GW, CF, YH, WP, and JZ were responsible for experimental design and sample collection, HS and YH were responsible for article data management, CL, LD, and JZ were responsible for writing and modifying.

Corresponding author

Correspondence to Wei Peng.

Ethics declarations

Ethics approval and consent to participate

The yaks in this study were all handled in strict accordance with the guidelines established by the Council for Animal Welfare of China. Additionally, the protocols for sample collection and animal handling were thoroughly reviewed and approved by the Faculty of Animal Policy and Welfare Committee of Northwest A&F University (FAPWCNWAFU, NWAFAC 1008).

Furthermore, the study was carried out in compliance with the ARRIVE guidelines.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

The original version of this article was revised: the ‘Availability of data and materials’, ‘Funding’ and ‘Ethics approval and consent to participate’ declarations were missing.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhao, H., Huang, Y., Shu, S. et al. Transcriptomics and metabolomics of blood, urine and ovarian follicular fluid of yak at induced estrus stage. BMC Genomics 25, 201 (2024). https://doi.org/10.1186/s12864-024-10079-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-024-10079-7

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us