Female yaks in good physical condition, with normal reproductive function, free of reproductive tract diseases, and provided with good husbandry management were examined rectally to ensure that there were no female yaks that had been impregnated due to natural mating under grazing conditions, as well as those that had not yet given birth. The simultaneous estrus programme used a combination of progesterone pessary + PGF_2α + GnRH, which was firstly buried with a progesterone pessary (Institute of Animal Husbandry Science and Technology, Aba Tibetan and Qiang Autonomous Prefecture, Sichuan Province, China), and the day of pessary release was set as day 0; the pessary was removed on day 7 and an intramuscular injection of sodium cloprostenol (PGF_2α 0.4 mg/head) (Ningbo Sansheng Biotechnology Co., Ltd., Zhejiang Province, China) was given; an injection of gona/gonadrin (GnRH 100 μg/head) (Ningbo Second Hormone Factory, Zhejiang Province, China) was given on day 9; and estrus was identified on day 10, together with her timed insemination. The artificial insemination was carried out with reference to the Agricultural Industry Standard of the People’s Republic of China “Technical Procedures for Artificial Insemination of Yaks” (NY/T 3979-2020).

Blood samples were tracked during anestrus (Y-A), estrus (Y-E) and pregnancy (Y-P). Whole blood was collected for transcriptome sequencing using anticoagulated (EDTA-containing) blood collection tubes, with 250 μL whole blood + 750 μL TRIzol LS per tube, and whole blood was collected for metabolome sequencing using anticoagulated (EDTA-containing) blood collection tubes, with 1 mL of blood per tube, and the blood samples were processed and stored in liquid nitrogen immediately after collection. To ensure that the blood samples were collected at the correct stage, the blood of six cows in their second year of calving was selected for sequencing. The samples were collected from the artificially inseminated herd of Yak Breeding and Promotion Service Centre of Qinghai Province (E101°22′, N37°15′), and all the animals were adult yaks aged 4–8 years.

Total RNA was extracted from blood samples using Trizol LS reagent (Thermo Fisher Scientific) according to the manufacturer’s instructions. Briefly, samples were homogenized, mixed with chloroform, and centrifuged to separate phases. The aqueous RNA-containing layer was collected, and RNA was precipitated with isopropanol, washed with 75% ethanol, and finally dissolved in nuclease-free water. RNA integrity and concentration were assessed, and qualified samples were stored at −80 °C for subsequent analysis.

RNA sequencing libraries were prepared from 5 μg of total RNA using the NEBNext Ultra II RNA Library Prep Kit (New England Biolabs, United States), in accordance with the manufacturer’s instructions. Final library quality was assessed using an Agilent Bioanalyzer, revealing an average insert size of 300 ± 50 bp. Paired-end sequencing (2 × 150 bp) was conducted on an Illumina Novaseq 6000 platform. The downstream raw data format was Fastq, and the Cutadapt software was used to perform quality control on the downstream raw data, removing the reads with joints, containing N (N indicates that the base information cannot be determined) in a proportion greater than 5%, and low quality (the number of bases with a quality value of Q ≤ 10 accounted for more than 20% of the whole read), to obtain the valid data (clean data). The pre-processed valid data (valid data) were then compared to the reference genome (Bos grunniens, genome version: v101) using the comparison software HISAT2 (13, 14). Mapped reads for each sample were assembled using StringTie (http://ccb.jhu.edu/software/stringtie/, version: stringtie-2.1.6). All transcriptomes were then merged using gffcompare software to reconstruct a comprehensive transcriptome. After the final transcriptome was generated, StringTie and ballgown were used to estimate the expression levels of all transcripts, and the expression levels of all transcripts were estimated by calculating the FPKM (Fragments per kilobase of exon model per million mapped reads, per thousand transcripts sequenced bases per million sequenced bases value, number of sequenced fragments) values were used to calculate mRNA abundance. Genes were screened in terms of both fold difference and significant level. Folds of difference FC ≥2 or FC ≤0.5 (i.e., |log₂FC| ≥1) and q value <0.05 (|log₂FC| ≥1 & q < 0.05) were used as the threshold criterion (no multiplicity of difference for multiple comparisons, and genes screened for q < 0.05 as statistically different between the multiple groups), and genes screened as a result were considered to be differentially expressed genes.

Principal component analysis (PCA) between the three groups of samples was performed using the R package vegan, and gene box plots and correlation clustering heat maps were drawn using the R package stats. The clustering heatmap was used to show the expression of genes in different samples. In order to better visualize the clustering expression pattern, log10 (FPKM+1) was used for non-biological replicates to show the gene expression. For biological replicates, the differential gene FPKM was used for gene expression display by Z-value (Zsample-i = [(FPKMsample-i)-Mean (FPKM of all samples)]/[Standard deviation (FPKM of all samples)]). Differentially expressed genes (DEGs) were enriched by GO (Gene Ontology) using the g:Profiler online website, and KEGG (Kyoto Encyclopedia of Genes and Genomes, Kyoto) was performed on DEGs using Kobas 3.0. Encyclopedia of Genes and Genomes functional analysis with a screening condition of p < 0.05.

The samples were thawed in the ice water bath and mixed by vortices for 30 s. Then a 200 μL aliquot of each individual sample were transferred to a 1.5 mL Eppendorf tube (Supplementary Table S1). After the addition of 20 μL of ISTD solution and 400 μL of methal, the samples were vortexed for 30 s, followed by centrifuged for 10 min, 12,000 rpm [RCF = 13,800 (×g), R = 8.6 cm], and 4 °C. Then 360 μL of water was added in the aliquot of 300 μL supernatant, followed by vortexed for 30 s. After centrifugation (15 min, 12,000 rpm, and 4 °C), a 570 μL aliquot of the supernatant was further purified with SPE. The SPE cartridges were washed with 200 μL of methanol, then equilibrated with 200 μL H₂O. After loading a sample, the SPE cartridges were washed with 200 μL 10% ACN/H₂O (v/v) and 200 μL Hexane, the flow-through fraction was discarded. The cartridge was then rinsed with 24 μL 90% ACN/H₂O (v/v). After this single-step SPE, 36 μL of water were added to the cluent. All the samples were shaked for 3 min, then the clear supernatant was subjected to UHPLC-MS/MS analysis.

Stock solutions were individually prepared by dissolving or diluting each standard substance to give a final concentration of 10 mmol/L. An aliquot of each of the stock solutions was transferred to a 10 mL flask to form a mixed working standard solution. A series of calibration standard solutions were then prepared by stepwise dilution of this mixed standard solution (containing isotopically-labelled internal standard mixture in identical concentrations with the samples).

The UHPLC separation was carried out using an ACQUITY UPLC-I/CLASS (Waters), equipped with a Waters ACQUITY UPLC BCH C8 column (100 × 2.1 mm, 1.7 μm, Waters). The mobile phase A was 0.5 mmol/L NH4F in water, and the mobile phase B was methal. The column temperature was set at 50 °C. The auto-sampler temperature was set at 10 °C and the injection volume was 10 μL. Waters Xevo TQ-S triple quadrupole mass spectrometer (Waters), equipped with an IonDrive Turbo V electrospray ionization (ESI) interface, was applied for assay development. Typical ion source parameters were: Capillary Voltages = 2.50 kV; Cone Voltages = 80 V (pos)/30 V (neg); Desolvation Temp = 550 °C; Desolvation Gas Flow = 1,100 L/Hr; Cone Gas Flow = 150 L/h; Nebuliser Gas Flow = 7.0 Bar. The MRM parameters for each of the targeted analytes were optimized using flow injection analysis, by injecting the standard solutions of the individual analytes, into the API source of the mass spectrometer. Several most sensitive transitions were used in the MRM scan mode to optimize the collision energy for each Q1/Q3 pair (Supplementary Table S2). Among the optimized MRM transitions per analyte, the Q1/Q3 pairs that showed the highest sensitivity and selectivity were selected as “quantifier” for quantitative monitoring. The additional transitions acted as “qualifier” for the purpose of verifying the identity of the target analytes. Water MassLynx V4.1 and Skyline were employed for MRM data acquisition and processing.

Weighted gene co-expression network analysis (WGCNA) target metabolomics part of the data as traits in association with the transcriptome. Based on the gene expression profiles of each breeding period, the selected gene set needs to be screened and filtered before WGCNA analysis to remove the low-quality genes that cause unstable effects on the results and to improve the accuracy of the network construction. The gene-developmental stage co-expression network was constructed using the WGCNA software package (version: 1.69). The correlation coefficients were transformed by a weighting function (adjacency function) to form an Adjacency Matrix, in which the elements are continuous. The fundamental premise is the principle of scale-free networks, which posits that the gene expression network adheres to the power function distribution characteristic of scale-free networks. The gene correlation network was made to conform to the scale-free distribution by successively weighting from 1 to 30 to find a suitable soft threshold. It was determined that the plateau threshold line set for this WGCNA analysis was 0.85 and the soft threshold was 16 (Figure 1A). To be more biologically meaningful, topological overlap measure (TOM) was used to calculate the degree of association between genes, and in addition to analysing the relationship between two genes, the connections between these two genes and other genes were also considered. The division of gene modules is based on the sparsity of connections between genes, and the TOM matrix (Similarity) is transformed into a dissimilarity matrix (Dissimilarity), which in turn constructs a hierarchical clustering tree (Figure 1B), using gene modules generated by the Dynamic Hybrid Shear Algorithm, and optimized merged gene modules [the minimum distance of merged modules (mergeCutHeight) is 0.25]. The minimum number of genes in a module (minModuleSize) was set to 30 for module detection. To identify period-specific modules, the correlation between each module and different breeding periods was explored using the WGCNA package. Modules closely associated with specific time points were screened as period-specific modules.

Bipartite sequencing of 18 samples was conducted on an Illumina Novaseq^™ 6000 platform using PE150 mode, generating 767,684,034 raw reads (average 42,649,113 per sample). Following the filtration process with Fastp, a total of 745,678,206 clean reads were obtained, constituting in excess of 97% of the original data set. Furthermore, it was observed that a minimum of 99.97% of bases exhibited quality scores of at least 20, with over 97% demonstrating scores of at least 30, and a global content of approximately 49%. The results obtained from this study suggest that the data is of a high quality and suitable for subsequent analysis. The alignment of the clean reads against the reference genome (Bos grunniens, v101) revealed that 85.79% were successfully mapped (Supplementary Table S3). Of these, 71% were found to be uniquely mapped to a single genomic location, while 26.05% aligned to multiple locations. The reads that were uniquely matched to a single location on the genome were used for the next data analysis.

In order to assess inter-group variability and reproducibility among the six replicates per group, principal component analysis (PCA) and correlation analysis were performed on all 18 samples. The findings of both analyses demonstrated clear separation between the three groups and high reproducibility among replicates, indicating stable and consistent results (Figures 2A,B). Gene-level expression analysis employing FPKM values (Supplementary Table S4) revealed broadly similar expression profiles across groups. The distribution of gene expression levels was quantified in the following manner: 17.27% of genes demonstrated low expression (FPKM <0.3), 51.99% exhibited moderate expression (0.3 ≤ FPKM < 15), 22.33% displayed intermediate expression (15 ≤ FPKM < 60), and 8.41% were highly expressed (FPKM ≥60). In summary, approximately half of the detected genes demonstrated moderately high expression levels.

The Y-A, Y-E and Y-P groups were divided into three comparison groups: Y-E vs. Y-A, Y-P vs. Y-E and Y-P vs. Y-A. Differentially expressed genes and transcripts were screened using |log₂FC| ≥1 & q < 0.05 as threshold criteria (Figures 2C,D). The volcano diagram illustrates the screened DEGs and their overall distribution (Supplementary Figure S1). Specific analyses showed that a total of 93 DEGs were identified between the Y-E and Y-A group, with 26 up-regulated and 67 down-regulated. Three hundred and nineteen DEGs were identified in the analyses of the Y-P and Y-A group, and 159 DEGs were identified. DEGs, 159 with up-regulated expression and 160 with down-regulated expression. In contrast, a total of 422 DEGs were identified in the comparison of Y-P and Y-E group, with 218 up-regulated and 204 down-regulated. In order to characterise the expression patterns that were observed across the physiological stages, the top 100 most significant differentially expressed genes (lowest q-value) from each comparison were subjected to hierarchical clustering (Figure 2E). The resulting heatmap revealed clear sample clustering by physiological stage and distinct gene expression profiles among groups.

GO enrichment analysis described the distribution of genes on Biological Process (BP), Cellular Component (CC) and Molecular Function (MF). The DEGs of the three control groups were analysed by GO through the g:Profiler online website (p < 0.05), and the detailed results are shown in Figure 2F. In terms of biological processes, DEGs in Y-E vs. Y-A group were significantly enriched in heme biosynthesis process, oxygen transport and erythrocyte differentiation; DEGs in Y-P vs. Y-A group were significantly enriched in oxygen transport, erythrocyte differentiation and cGMP-mediated signalling; DEGs in Y-P vs. Y-E group were significantly enriched in classical pathways such as synaptic pruning, antimicrobial peptide-mediated humoral immune response and complement activation and classical pathways. In terms of molecular function, DEGs from all three groups were significantly associated with oxygen binding, heme binding, iron binding and ammonium transport across membranes.

To further understand the biological functions of different genes, KEGG pathway enrichment analysis was performed to explore the biological processes involved in DEGs at three different physiological periods. The results of the analysis are presented in Figure 2G, and the screening criterion was set at p < 0.05. When comparing the Y-E vs. Y-A group, 93 DEGs were found to exhibit significant enrichment in 17 pathways, which were mainly focused on extracellular matrix (ECM)-receptor interactions. In the Y-P vs. Y-A group, 319 DEGs were significantly enriched in 26 pathways, such as IL-17 signalling pathway, p53 signalling pathway, oxytocin signalling pathway and NF-kB signalling pathway were significantly associated. In the Y-P vs. Y-E group, 422 DEGs were significantly enriched in 18 pathways, which were mainly related to follicular development and ovarian activities, such as PPAR signalling pathway, ECM-receptor interaction, IL-17 signalling pathway, relaxin signalling pathway and arginine biosynthesis were significantly associated. The three groups of DEGs were co-enriched into 71 signalling pathways, which were mainly related to oocyte maturation and ovulation, etc., such as progesterone-mediated oocyte maturation, HIF-1 signalling pathway, VEGF signalling pathway and PPAR signalling pathway.

The raw data contained 18 experimental samples from which 19 metabolites were extracted, and for better analysis of the data, individual metabolites were filtered and only metabolite data with no more than 50% nulls in a single group or no more than 50% nulls in all groups were retained. Missing values in the raw data were simulated, and the numerical simulation method was filled by multiplying the minimum value by a random number between (0.1, 0.5). After preprocessing 15 metabolites were retained.

The PCA scores of all the samples showed that the samples within the groups could be clustered together and there was a clear trend of separation among the three groups of samples, and the variances of PC1 and PC2 were 70.8 and 12.4%, respectively (Figure 3A), which indicated that there was a difference in the metabolites in the blood of the yaks in different breeding periods. And all samples were in the 95% confidence interval, indicating that there were no outliers in the analysed samples and there was a significant trend of separation among the three groups of samples. According to the OPLS-DA (orthogonal partial least squares discriminant analysis) model, the Y-P vs. Y-A and Y-P vs. Y-E models were better than the original model, with good robustness and no overfitting phenomenon; the model of the Y-E vs. Y-A group reached an R²Y of 0.763 in explaining the overall variability, which indicated that the model had successfully captured the variation patterns of steroid hormone metabolites (metabolites) in different physiological stages to a certain extent. Steroid hormone metabolites at different physiological stages (Supplementary Figure S2).

For the 15 steroid hormone target metabolites after sample pretreatment, the results of statistical analyses of metabolites were visualised in the form of volcano plots (Figures 3B–D). In order to improve the reliability and biological significance of the identified differentially expressed metabolites, the metabolites retained after fulfilling the dual criteria of VIP >1 & p < 0.05 were considered to be the more plausible and meaningful differentially expressed metabolites after multiple screening using a combination of VIP (variable importance in projection) values and statistical significance tests. Dehydroepiandrosterone was significantly up-regulated and aldosterone, cortisone, corticosterone, 11-deoxycortisol and 17α-hydroxyprogesterone were significantly down-regulated in the Y-E vs. Y-A group; progesterone, testosterone, dehydroepiandrosterone, 17α-hydroxyprogesterone, androstenedione, 11-deoxycortisol and oestrone were significantly up-regulated and cortisone was significantly down-regulated in the Y-P vs. Y-A group; and in the Y-P vs. Y-E group 11 -deoxycorticosterone, progesterone, testosterone, 17α-hydroxyprogesterone, androstenedione, pregnenolone, dehydroepiandrosterone and aldosterone were significantly up-regulated.

The transcriptome sequencing results of three groups of 17 blood samples were analysed using the WGCNA approach, and 20,537 genes were finally retained for the construction of weighted gene co-expression networks. In order to find the core genes in a more biological sense, the data of steroid hormone target metabolome were jointly analysed with the transcriptome data as traits, and different breeding periods were used as developmental periods to correlate traits with gene modules, and the results of correlation strengths and weaknesses between different modules and different breeding periods were obtained, by analysing the correlation between the traits and the modules and plotting the heatmap to visually show the correlation between the modules and the given trait. We analysed the correlation between traits and modules and drew heat maps to visualize the correlation between modules and given traits. Among them, the dark grey module had the highest correlation with group Y-A (cor = 0.71, p = 0.001); the turquoise module had a high negative correlation with group Y-E (cor = −0.54, p = 0.03); and group Y-P had the highest correlation with the turquoise module (0.85, 1 × 10⁻⁵), dark turquoise module (cor = 0.76, p = 4 × 10⁻⁴) and brown module (cor = 0.75, p = 5 × 10⁻⁴) with a strong positive correlation (Figure 4A). Therefore, each of these four modules was used as a period-specific module with its closely associated time point for the next analysis of the biological functions of the key modules and the screening of the core (Hub) genes. Dark grey module had a high correlation with 11-deoxycortisol, corticosterone, cortisol, and cortisone in the glucocorticoid family; turquoise module and brown module had high correlation with 17α-hydroxyprogesterone, androstenedione, dehydroepiandrosterone, luteinising hormone and testosterone in the steroid hormones (Figure 4B).

In order to further explore the potential regulatory relationships among genes in the period-specific modules and to mine the key regulatory genes in different periods, gene co-expression networks of highly correlated modules in different breeding periods were generated using the WGCNA approach. The period-specific co-expression networks were analysed by Cytoscape software and the genes located in the top 50 connectivity nodes were visualised. In these co-expression networks, the larger the shape where the node is located indicates a higher degree of connectivity, which can be used as a Hub gene in this network. From the visualised networks, it can be seen that the SLC25A39, MARCHF2, ENSBGRG00000026772 and ENSBGRG00000026382 genes are Hub genes in the dark grey module (Figure 5A); the MATR3 gene is a Hub gene in the turquoise module (Figure 5B); ALX1, RGS1, TENT5B, and ENSBGRG00000024929 genes are Hub genes in the turquoise, dark turquoise, and brown modules (Figure 5C).