All experimental procedures on live animals have been strictly reviewed and approved by the Institutional Animal Care and Use Committee at Southwest Minzu University (Chengdu, Sichuan, China), and all the experiments were carried out in accordance with the requirements of the directory of the Ethical Treatment of Experimental Animals of China. All Maiwa yaks were raised in Hongyuan County, Aba Tibetan and Qiang Autonomous Prefecture, Sichuan Province. All yaks were naturally grazed without feed supplementation, and calves fed on breast milk and natural forage. Fifteen Maiwa yaks were randomly selected from five different growth stages, including 0 d (0 day of age), 20 d (20 days of age), 60 d (60 days of age), 15 m (15 months of age), and adult (3 years of age). Each age group consisted of 3 healthy yaks with similar body weight, and the gender was selected randomly. Because of the large size of yaks, it is not suitable to use anesthetic in the sampling process, so a more extensive and humane method of electric shock was adopted in this experiment. To ameliorate suffering, yaks were electrically stunned (120 V dc, 12 s) before exsanguination, and sacrificed by bloodletting through the carotid artery and jugular vein while in a coma. After slaughter, the contents of each stomach chamber were emptied immediately and washed with saline several times until the attached contents are completely removed. Tissue blocks of approximately 1 cm × 1 cm × 1 cm were immediately taken from 4 stomachs with high-temperature sterilizing scissors, washed with normal saline and DEPC water successively, then put into the RNA-free frozen storage tube and frozen in liquid nitrogen until RNA was extracted.

All Maiwa yaks were accurately weighed and recorded before slaughter. After slaughter, the ruminant stomach was removed, and the connected parts of rumen, omasum, reticulum and abomasum were ligated, and the four gastric compartments were separated. In the rumen, we found no solid contents in the 0 d group, very little pasture and a lot of milk in the 20 d group, some pasture and a lot of milk in the 60 d group, and only pasture in the 15 m and adult groups. The weights of the rumen, reticulum, omasum and abomasum were separately measured and recorded after completely emptying digesta and fluid. The stomach index was calculated using the following formula: Stomach index = single stomach weight (g)/ruminant stomach weight (g) × 100%.

Total RNA was extracted from each sample using the mirVana miRNA Separation kit (Ambion-1561) according to the manufacturer’s instructions. RNA concentration and integrity were measured and assessed using the NanoDrop 2000 (Thermo Fisher Scientific) and the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States). The samples with RNA Integrity Number (RIN) ≥ 7 were further analyzed.

The libraries were constructed following the manufacturer’s protocol at Shang-hai OE Biotech (Shanghai, China), using TruSeq Stranded mRNA LTSample Prep Kit (Illumina, San Diego, CA, United States). In short, these mRNA samples were enriched with the Oligo (dT) magnetic beads and fragmented using the breaking reagent. Then first-strand cDNA and second-strand cDNA were successively synthesized and purified. The purified double-stranded cDNA was end-repaired, A-tailed, ligated to sequencing adapters, amplified and purified. These libraries were purified using AMPure XP beads and analyzed for purity and size by Agilent 2100 Bioanalyzer. These libraries were sequenced on the Illumina HiSeq X Ten and 150 bp paired-end reads were obtained.

Trimmomatic software was used to process the raw data (raw reads) for quality control (17). The clean reads were obtained by removing the raw reads containing ploy-N and the low-quality reads. The high-quality clean reads were then mapped to the yak reference genome (GCF_000298355.1) using hisat2 (18) (Version 2.2.1). ¹ The cufflink software was used to calculate the FPKM value of each gene (19, 20), and the htseq-count software was used to count the number of reads mapped to each gene (21). For quantification of transcription levels, bowtie2 (22) and eXpress (23) were used to calculate FPKM and read counts value of each transcript. The StringTie software was adopted to reassemble the reads (24). On this basis, using cuffCompare software, the reference genome was compared with the existing annotated genes to identify the gene structure extension and novel transcripts. The ASprofile was used to perform the alternatively splicing analysis of differentially regulated transcripts isoforms or exons (25). SNP and INDEL were called using samtools (26) and bcftools (27), and the details were shown on samtools webpage. ² Then the effects of variants on genes were annotated and predicted by the sniff (28).

Based on the DESeq (2012) (29) R package functions estimateSizeFactors and nbinomTest, differential expression analysis was performed and differentially expressed genes (DEGs) were identified. Adjusted p value < 0.05 and fold Change ≥ 2 or fold Change ≤ 0.5 was set as the threshold for significantly differential expression. To explore gene and transcripts expression patterns, hierarchical cluster analysis of DEGs was executed. GO enrichment and KEGG pathway enrichment of DEGs were analyzed according to the methods in the references (30).

The STEM clustering algorithm was used to explore the relationship between temporal gene expression patterns and yak stomach tissues (rumen, reticulum, omasum, and abomasum) during the five developmental stages (31). Time-series analysis mainly studied the dynamic expression of genes, measuring a range of processes strongly related to developmental stages. Significantly enriched model profiles were displayed in different colors (Benjamini-Hochberg-adjusted p-values ≤ 0.05). The corrected p-values were sorted from small to large. Colored model profiles (except white) indicated that the temporal trends of mRNAs were statistically significant. Model profiles with the same color belonged to the same cluster.

The co-expression networks of rumen-reticulum-omasum-abomasum mRNAs were established by the WGCNA package of R software through the developmental stages (32). After eliminating samples with outliers samples, the Pearson’s correlation coefficient between any two genes was determined in the gene set, and the correlation coefficients matrix was builded. The appropriate threshold value (β value) was selected to measure the weighted power exponent of the correlation coefficient matrix, and the adjacency matrix was established. Next, the topological overlap matrix was established and employed to the connections between different genes. Based on related traits, the gene modules were preliminarily divided by hierarchical clustering analysis, and the eigengenes were obtained. According to the similarity of eigengenes, these modules were combined to form a final module for further study (33). The module membership (MM) was calculated by the WGCNA function signed KME, which correlated the module eigengene (ME) with gene expression values, thereby quantifying how close a gene was to a given module. Gene significance (GS) referred to the correlation between a single gene and a biological trait. The summation of adjacency performed for all genes in a particular network was calculated as the intramodular connectivity (K.in). We first selected genes with GS > 0.25 and MM > 0.95. Then, we screened for genes with the top50 of the number of gene interaction nodes as hub genes. Visualization analysis of gene interaction network was performed using cytoscape software v3.8.2.

To verify the reliability of RNA-seq data, 16 genes were randomly selected for RT-qPCR quantification. Total RNA was isolated from the rumen, reticulum, omasum and abomasum of yaks according to the instructions of the mirVana miRNA Isolation Kit (Ambion). Reverse transcription was used to synthesize cDNA according to the instructions of the PrimeScript RT reagent Kit with gDNA Eraser (TaKaRa, Dalian, China). The transcription levels of mRNAs (S100A12, KRT4, KRT6A, ACTG2, A2ML1, S100A9, FTH1, DES, RPLP0, ANXA1, RPS3A, TFF1, GKN1, ATP4B, PGC, and RPS8) were verified by LightCycler 96 System (Roche, United States). The forward and reverse primers for RT-qPCR are shown in Supplementary Table S1. RT-qPCR was performed in a 96-well plate containing 20 μL mixture per well, including 10 μL of TB Green Premix Ex Taq^™ II (Tli RNaseH Plus) (TaKaRa, Dalian, China), 2.0 μL of cDNA template, 0.5 μL of forward and reverse primer (10 μM), and 7.0 μL of RNase Free ddH2O. The running conditions of RT-qPCR were set as follows: 95°C for 30 s (pre-denaturation), 40 cycles of amplification (95°C for 5 s, 60°C for 30 s) and 95°C for 5 s, 60°C for 1 min (Melting Curves). For each time point, the gene validation was conducted in triplicate. The 2^−∆∆Ct method was used to calculate the expression level of each verified gene at each time point.

The increase in volume and weight of the ruminant stomach is one of the important indexes to measure the development stage of yak. This study measured the weight of four stomachs at different developmental stages. At 0 d, the weight of the abomasum was the largest, followed by the rumen, and the weight of the abomasum was markedly heavier than that of the other three stomachs (p < 0.05, Figure 1A). At 20 d, the weight of abomasum was almost equal to that of rumen, and there was no significant difference between them (p > 0.05), but the weight of abomasum was significantly heavier than that of reticulum and omasum (p < 0.05, Figure 1A). At 60 d, the weight of rumen was the heaviest, which was significantly heavier than that of abomasum (p < 0.05), and the abomasum weight was significantly heavier than that of reticulum and omasum (p < 0.05, Figure 1A). At 15 m, the weight of rumen was significantly heavier than that of omasum (p < 0.05), and the weight of reticulum and abomasum was the lightest. At adult group, the weight of rumen was the heaviest, which was significantly heavier than that of omasum (p < 0.05), and the weight of abomasum was heavier than that of reticulum (p > 0.05, Figure 1A). According to the yak stomach index at 0 d and 20 d, the ratio of abomasum was higher than that of the forestomach. At 60 d, the ratio of the rumen was the highest, followed by the abomasum. At 15 m and adult groups, the ratio of rumen was the largest, followed by omasum, abomasum and reticulum (Figure 1B).

To identify the dynamic changes of four stomachs in transcriptome profiles at different growth stages, the transcriptome sequencing of 60 samples was performed, and 2906.51 Mb clean reads was gained (Supplementary Table S2). For each sample, the effective data volume was distributed in 6.08–7.34 G, the Q30 base distribution was 93.36–96.03%, and the average GC content was 51.05%. The genome alignment of each sample was gained by aligning the reads to the reference genome, and the alignment rate was 89.63–96.59%. The above results indicated that the sequencing data were qualified and in-depth enough for further analysis.

To understand the regulatory roles of mRNAs of rumen, reticulum, omasum and abomasum, we analyzed the expression profile of mRNAs in four stomachs at 0 d, 20 d, 60 d, 15 m and adult. The study was able to evaluate the dynamic changes of mRNA expression from non-ruminant to ruminant stage, and screen and identify important mRNAs related to the growth and metabolism in ruminant stomach of yaks. First, we designed two comparative methods, namely, closed groups (20 d vs. 0 d, 60 d vs. 0 d, 15 m vs. 0 d and adult vs. 0 d), and consecutive groups (20 d vs. 0 d, 60 d vs. 20 d, 15 m vs. 60 d, and adult vs. 15 m) to characterize the DEmRNAs during development. Among these comparison groups in the rumen, 3,240 mRNAs were defined as DEmRNAs (Figure 2A).

We attempted to evaluate the potential functions of DEmRNAs between two closed time points in the rumen. Compared with 0 d, 532, 829, 1,441, and 738 mRNAs were up-regulated at 20 d, 60 d, 15 m and adult groups, respectively (Figure 2B). KEGG pathway analysis showed that, compared with 0 d, the significantly up-regulated pathways in the other four stages are all concerned with lipid metabolism, carbohydrate metabolism, and immune response including such processes as steroid biosynthesis, butanoate metabolism, propanoate metabolism, hematopoietic cell lineage (Supplementary Table S3). Among these DEmRNAs, 297 were identified as common DEmRNAs throughout rumen development (Figure 2C, Supplementary Table S4). These common DEmRNAs played essential roles in rumen development, immune response, carbohydrate metabolism and lipid metabolism. KEGG pathway analysis indicated that these common DEmRNAs were primarily enriched in ECM-receptor interaction, viral protein interaction with cytokine and cytokine receptor, cell adhesion molecules, propanoate metabolism, steroid biosynthesis, and hematopoietic cell lineage, etc. (Figure 2D, Supplementary Table S5).

To study the role of DEmRNAs in transforming rumen function, in addition to 20 d vs. 0 d, three comparison groups including 60 d vs. 20 d, 15 m vs. 60 d and adult vs. 15 m were added. Compared to 0 d, the up-regulated pathways at 20 d were primarily related to steroid biosynthesis, steroid hormone biosynthesis, pentose and glucuronate interconversions, hematopoietic cell lineage, butanoate metabolism, mineral absorption, cell cycle, and the down-regulated pathways were observed to be principally related to MAPK signaling pathway, cortisol synthesis and secretion, central carbon metabolism in cancer, and human T-cell leukemia virus 1 infection. Compared to 20 d, up-regulated pathways at 60 d were chiefly related to steroid hormone biosynthesis, drug metabolism—cytochrome P450, retinol metabolism, and pentose and glucuronate interconversions, and the down-regulated pathways were primarily related to aldosterone-regulated sodium reabsorption, pancreatic secretion and sphingolipid metabolism. Compared to 60 d, the up-regulated pathways at 15 m were mainly involved in kaposi sarcoma-associated herpesvirus infection, IL-17 signaling pathway, glycine, serine and threonine metabolism, and NF-kappa B signaling pathway. The down-regulated pathways were primarily involved in such metabolic pathways as ECM-receptor interaction, and protein digestion and absorption. Compared to 15 m, the significantly up-regulated pathways at adult were primarily related to antifolate resistance, ABC transporters, bile secretion, purine metabolism and cAMP signaling pathway, and the down-regulated pathways were engaged mainly in such immune processes as pertussis, systemic lupus erythematosus, and complement and coagulation cascades (Supplementary Table S6). These up- and down-regulated mRNAs at different development stages may be involved in regulating the beginning or ending of growth and development, and/or important physiological processes at specific developmental stages. We also focused on the changing rules of the 20 most abundant mRNAs in the rumen. Notably, several of the most abundant mRNAs in the rumen originate from protein-coding genes with pivotal roles in rumen growth, immune response, and translation (e.g., S100A12, S100A9, RPS2, KRT4, RPS8, KRT6A, DES, A2ML1, EEF1A1, ACTG2, RPS3A, and RPLP0) (Figure 2E, Supplementary Table S7). KEGG analysis was performed to further explore the switching of fatty acid metabolism-related functions at different developmental stages in the rumen. The results showed that 20 days is an important metabolic time point, and many fatty acid pathways are significantly enriched at 20 days, and there are different degrees of enrichment at 60 days, 15 months and adult groups (Supplementary Figure S1, Figure 2F). There were no differences in some metabolic pathways in the 60 d vs. 20 d, 15 m vs. 60 d and adult vs. 15 m groups, including fatty acid biosynthesis, fatty acid elongation, synthesis and degradation of ketone bodies, steroid biosynthesis, pyruvate metabolism, propanoate metabolism, synthesis and degradation of ketone bodies, steroid biosynthesis, pyruvate metabolism, propanoate metabolism, butanoate metabolism (Figure 2F).

The reticulum acts mechanically to reduce the ingesta to fine particles further. In the seven comparison groups, we identified 3,968 DEmRNAs to gain an overview of DEmRNAs in the reticulum (Figure 3A). Among these closed time points, compared to 0 d, we detected 751, 1,029, 1,541, and 1,351 DEmRNAs upregulated at 20 d, 60 d, 15 m and adult (Figure 3B). At these four closed time points, KEGG analysis showed that the up-regulated pathways were substantially enriched in immune-related signaling pathways, including hematopoietic cell lineage, cell adhesion molecules, cytokine-cytokine receptor interaction, viral protein interaction with cytokine and cytokine receptor, etc. Additionally, between 20 d and 60 d, significantly up-regulated pathways were involved in retinol metabolism, bile secretion, pentose and glucuronate interconversions, ascorbate and aldarate metabolism, steroid hormone biosynthesis, and butanoate metabolism, and the down-regulated pathways were primarily involved in IL-17 signaling pathway, viral protein interaction with cytokine and cytokine receptor, neomycin, kanamycin and gentamicin biosynthesis, african trypanosomiasis, and glycerophospholipid metabolism. Compared to 60 d, up-regulated pathways at 15 m were predominantly related to metabolic pathways including mineral absorption, steroid hormone biosynthesis, and arachidonic acid metabolism, and the down-regulated pathways are mainly related to ECM-receptor interaction, protein digestion and absorption, PI3K-Akt signaling pathway, etc. Compared to 15 m, the up-regulated pathways at adult were primarily related to hematopoietic cell lineage, Ras signaling pathway, MAPK signaling pathway, and PI3K-Akt signaling pathway, and the down-regulated pathways were primarily related to peroxisome, primary bile acid biosynthesis, steroid biosynthesis, arginine and proline metabolism, and arachidonic acid metabolism (Supplementary Table S8).

During the entire development process, compared to 0 d, 545 were regarded as the common DEmRNAs of the reticulum (Figure 3C), showing that the KEGG pathways were chiefly focused on lysosome, phagosome, viral protein interaction with cytokine and cytokine receptor, ECM-receptor interaction, cell adhesion molecules, cytokine-cytokine receptor interaction, primary immunodeficiency, TCA cycle, valine, leucine and isoleucine degradation, propanoate metabolism, and hematopoietic cell lineage, etc. (Figure 3D, Supplementary Table S9). Furthermore, among the 20 most abundant mRNAs expressed in the reticulum, several mRNAs mainly derived from protein-coding genes with the critical roles in reticulum development, immune response and mineral absorption are highly expressed (e.g., S100A12, KRT4, RPS2, KRT6A, ACTG2, A2ML1, EEF1A1, RPS8, FTH1, RPLP0) (Figure 3E, Supplementary Table S10).

The omasum is a compact spherical organ which composed of the omasal canal and omasal body, but its function is not fully known. A total of 3,975 DEmRNAs were identified among these comparison groups (Figure 4A). In four closed groups, compared with 0 d, 677, 771, 1,330, and 1,334 DEmRNAs were up-regulated at 20 d, 60 d, 15 m, and adult (Figure 4B), and the number of down-regulated DEmRNAs was 190, 286, 1,222, and 1,203, respectively. The number of upregulated DEmRNAs was 74, 294, and 306 in three consecutive groups, while the downregulated DEmRNAs were 76, 750, and 161, respectively (Figure 4B).

During the entire development stages, compared with 0 d, 428 DEmRNAs were considered as common DEmRNAs, revealing that the KEGG pathways were significantly enriched in the cell adhesion molecules, ECM-receptor interaction, cytokine-cytokine receptor interaction, primary immunodeficiency, staphylococcus aureus infection, steroid biosynthesis, retinol metabolism, arachidonic acid metabolism, hematopoietic cell lineage, and gastric acid secretion, etc. (Figures 4C,D, Supplementary Table S11). Notably, among the most abundant mRNAs expressed in the omasum, several mRNAs mainly were related to the key roles of omasum developmentand immune response (e.g., KRT4, S100A12, S100A9, RPS2, KRT6A, A2ML1, EEF1A1, RPS8, ACTG2, ANXA1) (Figure 4E, Supplementary Table S12).

Abomasum is the real “stomach,” which can secrete digestive enzymes. In addition to the digestive function of protein in feed, the gastric juice of abomasum can also kill the microorganisms in chyme and provide nutrition for ruminants. Among the seven comparison groups, we clustered 2,638 identified DEmRNAs to obtain an overview of DEmRNAs in the abomasum (Figure 5A). We detected 499, 874, 903, and 1,053 upregulated DEmRNAs at 20 d, 60 d, 15 m and adult in four closed time points (Figure 5B). At these four closed time points, compared with 0 d, KEGG analysis showed that the up-regulated pathways were remarkably enriched in the immune-related pathways, including cell adhesion molecules, hematopoietic cell lineage, natural killer cell-mediated cytotoxicity, Th17 cell differentiation, and Th1 and Th2 cell differentiation, etc. However, in two comparison groups (20 d vs. 0 d and 60 d vs. 0 d), the down-regulated KEGG pathways were primarily related to such protein metabolism and disease-related metabolic pathways as protein processing in endoplasmic reticulum, parkinson disease, ribosome biogenesis in eukaryotes, prion disease, chemical carcinogenesis-DNA adducts and glutathione metabolism, etc. The down-regulated KEGG pathways of the other two closed time points (15 m vs. 0 d and adult vs. 0 d) were all related to relaxin signaling pathway, protein digestion and absorption, amoebiasis, and AGE-RAGE signaling pathway in diabetic complications, and ECM-receptor interaction. Moreover, between 20 d and 60 d, remarkably up-regulated pathways were principally involved in antigen processing and presentation, steroid hormone biosynthesis, DNA replication, natural killer cell mediated cytotoxicity and cell cycle, and down-regulated pathways are mainly related to drug metabolism-cytochrome P450, aldosterone-regulated sodium reabsorption, and cGMP-PKG signaling pathway, etc. Compared to 60 d, the up-regulated pathways at 15 m were primarily enriched in drug metabolism-cytochrome P450, bile secretion, and the down-regulated pathways are primarily related to ECM-receptor interaction, cell cycle, protein digestion and absorption, cAMP signaling pathway, and DNA replication. Compared to 15 m, the up-regulated pathways at adult are primarily related to cell cycle, DNA replication, oocyte meiosis and progesterone-mediated oocyte maturation, and down-regulated pathways are primarily related to PPAR signaling pathway, chemical carcinogenesis-DNA adducts, pentose and glucuronate interconversions, retinol metabolism, and steroid hormone biosynthesis (Supplementary Table S13).

During the entire development stages, compared to 0 d, 332 DEmRNAs were identified as common DEmRNAs, indicating that KEGG pathways were primarily involved in the phagosome, viral protein interaction with cytokine and cytokine receptor, ECM-receptor interaction, cell adhesion molecules, primary immunodeficiency, arachidonic acid metabolism, hematopoietic cell lineage, and antigen processing and presentation (Figures 5C,D, Supplementary Table S14). Remarkably, several of the most abundant DEmRNAs in abomasum from protein-coding genes played pivotal roles in abomasum development, gastric acid secretion and digestive enzyme activity (e.g., TFF1, TMSB10, GKN1, ATP4B, RPS2, PGC, TFF2, EEF1A1, TPT1, ATP4A, RPS3A, RPLP0, and RPS8) (Figure 5E, Supplementary Table S15).

There are few reports on association analysis of the mRNAs of the four stomachs in yak. To identify the key mRNAs regulating differences in growth and functional metabolism, we evaluated the expression levels and analyzed differential expression of mRNAs in four stomachs (Figure 6A). Compared to the abomasum, we found that 3,687, 3,861, and 3,943 mRNAs were differentially expressed in rumen, reticulum, and omasum, respectively (Figure 6B). In addition, we also found that there were relatively few differential genes among rumen, reticulum and omasum, including 205 DEmRNAs between rumen and reticulum, 598 DEmRNAs between rumen and omasum, and 585 DEmRNAs between reticulum and omasum (Figure 6B). To further clarify the expression pattern of DEmRNAs among four stomachs, hierarchical cluster analysis was performed, and the expression heatmap of DEmRNAs was generated. The results displayed that the expression patterns of these mRNAs were remarkably different between forestomach and abomasum in adult groups (Figure 6C). KEGG enrichment analysis illustrated that these DEmRNAs were principally focused on axon guidance, focal adhesion, gastric acid secretion, cAMP signaling pathway, cGMP-PKG signaling pathway, arachidonic acid metabolism, PI3K-Akt signaling pathway, and regulation of actin cytoskeleton, etc. (Figure 6D, Supplementary Table S16). Therefore, these DEmRNAs were potential candidates for the regulation of development and regeneration, immune response, lipid metabolism and regulation of actin cytoskeleton.

According to the dynamic expression patterns of the mRNAs at the five developmental periods, all identified mRNAs in rumen, reticulum, omasum and abomasum were classified into 6, 4, 7, and 5 cluster profiles, containing 12, 10, 12, and 9 significantly enriched model profiles, respectively (Figures 7A–D). Among them, modules 9, 19, 38, and 41 were significantly enriched in rumen, reticulum, omasum and abomasum, and modules 46 and 49 were only significantly enriched in the rumen. Modules 9, 38, and 41 displayed a gradual increase or decrease trend. Modules 2 and 15 showed an “wavy” type expression patterns that were firstly reducing, then gradually increasing, and then decreasing, while modules 28 and 42 showed almost the opposite expression patterns. Modules 37, 39, 43, 44, and 46 showed the expression patterns that gradually increased and then progressively reduced, while modules 0, 6, 8, 14, 19, and 24 showed the opposite expression patterns. Finally, we observed colored modules, and considered the most significant three modules (modules 9, 38, and 41). Depending on the degree of significance, GO enrichment analysis of module profile 9 in four stomachs showed that the genes were mainly enriched for heparin binding, extracellular matrix structural constituent, proteinaceous extracellular matrix, extracellular matrix organization, calcium ion binding, cell adhesion (Supplementary Figures S2A–D). Module 38 showed that 20 days was an important node for regulating immune response, metabolism-related hormone levels, and enzymatic activity in four stomachs tissues (Supplementary Figures S2E–H). Module 41 was participated in fatty acid beta-oxidation, tricarboxylic acid cycle and pyruvate metabolic process (rumen), motor activity, mitochondrial matrix, fatty acid beta-oxidation and non-canonical Wnt signaling pathway (reticulum), negative regulation of canonical Wnt signaling pathway, actin filament, AMPA glutamate receptor complex, and channel regulator activity (omasum), immune response, response to sucrose, phosphatidylinositol-4,5-bisphosphate 3-kinase activity, and transmembrane receptor protein tyrosine kinase signaling pathway (abomasum) (Supplementary Figures S2I–L). In addition, 147, 101, 127, and 93 KEGG pathways of DEmRNAs were enriched in profile 41 in rumen, reticulum, omasum and abomasum, respectively (Supplementary Table S17). In rumen and reticulum, DEmRNAs were involved in valine, leucine and isoleucine degradation, citrate cycle (TCA cycle), propanoate metabolism, pyruvate metabolism, butanoate metabolism and PPAR signaling pathway, and these pathways were closely associated with amino acid metabolism, carbohydrate metabolism and lipid metabolism. In omasum, DEmRNAs participated in the focal adhesion, purine metabolism, ECM-receptor interaction, Wnt signaling pathway and hematopoietic cell lineage, and these pathways were related to cell growth, amino acid metabolism and immune regulation. In abomasum, DEmRNAs were mainly involved in the intestinal immune network for IgA production, cytokine-cytokine receptor interaction, Fc epsilon RI signaling pathway, sphingolipid signaling pathway, adipocytokine signaling pathway, fatty acid biosynthesis, and protein digestion and absorption, and these pathways were closely associated with intestinal immune regulation, protein metabolism and lipid metabolism (Supplementary Table S18).

A total of 18,575 genes were collected in the original data. After filtering out genes with low expression fluctuation (standard deviation ≤ 0.5), 5,486 genes were eventually left. By WGCNA analysis, a total of 5,486 genes were categorized into 10 modules (the gray module was excluded) (Figure 8A). In module-sample correlation analysis, we found the turquoise and light green modules had the highest correlations with abomasum and reticulum (Figure 8B). We also found that the turquoise and blue modules were sample-specific. The turquoise module was positively correlated with the abomasum, and negatively correlated with rumen, reticulum and omasum. The blue module was the opposite of the turquoise module. Genes in the turquoise (1568) and blue (1451) modules were expressed more highly than the other modules. Among them, CCKBR, KCNQ1, FER1L6, and A4GNT were the hub genes of the turquoise module (Figure 8C); PAK6, TRIM29, ADGRF4, TGM1, and TMEM79 were the hub genes of the blue module (Figure 8D). Through GO enrichment analysis, turquoise and blue modules were principally participated in the biological process (48.31–63.61%), followed by molecular function (22.15–26.69%), and least by cellular components (14.24–25.00%) (Supplementary Tables S19, S20). In the turquoise module, KEGG enrichment analysis indicated that these co-expressed genes were involved in gastric acid secretion, sphingolipid metabolism, ether lipid metabolism, pancreatic secretion, and glycosphingolipid biosynthesis-globo and isoglobo series. In the blue module, the co-expressed genes were participated in pancreatic secretion, riboflavin metabolism, pantothenate and CoA biosynthesis, proximal tubule bicarbonate reclamation, regulation of actin cytoskeleton, and starch and sucrose metabolism (Supplementary Tables S21, S22).

Sixteen mRNAs—S100A12, KRT4, KRT6A, ACTG2, A2ML1 (rumen)—S100A9, ACTG2, FTH1, DES, RPLP0 (reticulum)—ANXA1, RPS3A, RPLP0, A2ML1, KRT4 (omasum)—TFF1, GKN1, ATP4B, PGC, RPS8 (abomasum)—were randomly selected for reliable verification of RNA-seq data in yak. The selected genes were quantitatively analyzed by RT-qPCR at five developmental periods. The results presented that RT-qPCR data was consistent with RNA-seq data, revealing that RNA-seq data was dependable (Supplementary Figures S3–S6).