Yili Liu1†
Qi Min- 1Key Laboratory of Qinghai-Tibetan Plateau Animal Genetic Resource Reservation, College of Animal and Veterinary Sciences, Southwest Minzu University, Chengdu, China
- 2Institute of Qinghai-Tibetan Plateau, Southwest Minzu University, Chengdu, China
The development of the four stomachs of yak is closely related to its health and performance, however the underlying molecular mechanisms are largely unknown. Here, we systematically analyzed mRNAs of four stomachs in five growth time points [0 day, 20 days, 60 days, 15 months and 3 years (adult)] of yaks. Overall, the expression patterns of DEmRNAs were unique at 0 d, similar at 20 d and 60 d, and similar at 15 m and adult in four stomachs. The expression pattern in abomasum was markedly different from that in rumen, reticulum and omasum. Short Time-series Expression Miner (STEM) analysis demonstrated that multi-model spectra are drastically enriched over time in four stomachs. All the identified mRNAs in rumen, reticulum, omasum and abomasum were classified into 6, 4, 7, and 5 cluster profiles, respectively. Modules 9, 38, and 41 were the most significant three colored modules. By weighted gene co-expression network analysis (WGCNA), a total of 5,486 genes were categorized into 10 modules. CCKBR, KCNQ1, FER1L6, and A4GNT were the hub genes of the turquoise module, and PAK6, TRIM29, ADGRF4, TGM1, and TMEM79 were the hub genes of the blue module. Furthermore, functional KEGG enrichment analysis suggested that the turquoise module was involved in gastric acid secretion, sphingolipid metabolism, ether lipid metabolism, etc., and the blue module was enriched in pancreatic secretion, pantothenate and CoA biosynthesis, and starch and sucrose metabolism, etc. Our study aims to lay a molecular basis for the study of the physiological functions of rumen, reticulum, omasum and abomasum in yaks. It can further elucidate the important roles of these mRNAs in regulation of growth, development and metabolism in yaks, and to provide a theoretical basis for age-appropriate weaning and supplementary feeding in yaks.
1. Introduction
The yak (Bos grunniens) is a rare and valuable domesticated ruminant, which is mainly distributed in the Qinghai-Tibet Plateau and the adjacent high mountains above 3,000 meters. According to the National Catalogue of Livestock and Poultry Genetic Resources (2022 Edition), there are currently 21 yak species in China, including 19 local breeds and 2 cultivated breeds. They are mainly distributed in Sichuan, Tibet, Qinghai and other places in China, providing local herders with meat, milk, skin, hair, fuel and income. Compared with ordinary cattle, the meat, milk and processed products of yak have higher nutrient content and higher edible value, so yak is known as the “boat on the plateau” (1, 2). Maiwa yaks, one of the 19 local breeds, mainly rely on natural pasture with limited supplementation of concentrated feed. In nature, yaks graze on natural grassland all year round, but the lack of foraging resources poses a serious challenge to their survival in cold season (3). Therefore, yak ruminant stomach should have evolved a potential regulatory mechanism to adapt to the harsh plateau environment. Thus, transcriptome studies on yak ruminant stomach are of great significance to explore the development and metabolism-related genes and their expression patterns.
The ruminant stomach is composed of four compartments: rumen, reticulum, omasum, and abomasum, and each compartment has different physiological functions (4). It is well known that rumen epithelium contributes greatly to nutrient transport, absorption and metabolism of ruminant livestock (5, 6). Transcriptome sequencing is widely used to study the molecular mechanism of rumen development and metabolism (7–9). The gene expression pattern in the rumen is influenced by many factors, such as physiological state (10), RFI (residual feed intake, a measure of feed efficiency) (11), forage quality and particle size (8, 12), and age (13), etc. Ruminal epithelial cells are prone to acidosis, which further leads to rumen epithelial insufficiency, erosion and ulceration. The signal pathways in the rumen tissues of acidosis are enriched to cell signal transduction and morphogenesis, suggesting that rumen acidosis may affect the development of rumen epithelium in Holstein bull calves (10). Transcriptome analysis of rumen with different RFI has shown that the differentially expressed genes and pathways screened are mainly related to nutrient transport, cell growth and proliferation, immune function, inflammation, and apoptosis in beef cattle (11). In addition, nutrients and particle size are also important factors that alter the expression of genes participated in cell proliferation and apoptosis and complement complex. However, the related research on the reticulum, omasum and abomasum of yaks mainly focuses on the diversity of microorganisms, differences in flora and the crucial roles of these microorganisms in the whole growth and development (14, 15). Few reports indicate that the ruminant stomach, particularly the rumen and abomasum, show markedly different growth and developmental processes after birth and subsequent ruminant stages in Holstein cattle (16). From birth to adulthood, the four stomachs of yaks have significant changes in size and function, which play important physiological functions, respectively (4). As far as we know, there are almost no reports on the combined analysis of gene expression patterns in rumen, reticulum, omasum, and abomasum tissues of yaks at different developmental stages.
In this study, to further elucidate the molecular basis of physiological functions of rumen, reticulum, omasum and abomasum, we first compared mRNA expression levels of ruminant stomach at five developmental stages. At the transcriptome level, we systematically analyzed the differentially expressed genes, spatiotemporal expression patterns, and potential biological functions in four stomach tissues of yaks. STEM analysis was used to investigate the dynamic expression patterns of mRNAs, and a network diagram of rumen-reticulum-omasum-abomasum interactions was constructed with significant changes over time. In addition, screening of mRNAs affecting the development and metabolism of the four stomachs during the whole developmental stage of yaks can significantly help to clarify biological and metabolic functions and relative importance. This finding will help researchers to further clarify the physiological function and regulatory mechanism of ruminant stomach in yaks, and lay a theoretical foundation for subsequent studies on the efficient and precise regulation of ruminant stomach in terms of development and metabolism, and contribute to the further development of animal husbandry on the Tibetan Plateau.
2. Materials and methods
2.1. Animals and sample collection
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.
2.2. Organism index measurement
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%.
2.3. RNA extraction
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.
2.4. Library construction and sequencing
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.
2.5. Read mapping and transcriptome assembly
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).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.2 Then the effects of variants on genes were annotated and predicted by the sniff (28).
2.6. Differential expression analysis and functional enrichment
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).
2.7. Time-series analysis
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.
2.8. Co-expression networks
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.
2.9. Quantitative real-time PCR analysis
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.
3. Results
3.1. Change in yak stomach index
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).
Figure 1. Analysis of weight changes in four stomachs of yaks at five developmental stages. (A) Weight variation of four stomachs in yak at different ages. N = 3/group, the abscissa shows the age, and the ordinate shows the weight (g). (B) Analysis of the ratio of individual stomach weight to total stomach weight of yaks in different stages. N = 3/group, the abscissa represents the age, and the ordinate represents the ratio.
3.2. Sequencing and de novo transcriptome assembly
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.
3.3. Expression patterns of DEmRNAs in rumen tissue
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).
Figure 2. Expression patterns of DEmRNAs in rumen. (A) Heatmap of DEmRNAs in the seven compared groups. (B) Number of DEmRNAs in rumen. (C) The number of common DEmRNAs in rumen. (D) KEGG pathway analysis of common DEmRNAs in rumen. The top 20 enriched KEGG pathways ranked by p-values are shown. (E) Top 20 mRNAs expressed in rumen. (F) KEGG enrichment bubble plot of DEmRNAs involved in fatty acid metabolism in four consecutive groups.
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).
3.4. Expression patterns of DEmRNAs in reticulum tissue
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).
Figure 3. Expression patterns of DEmRNAs in reticulum. (A) Heatmap of DEmRNAs in the seven compared groups (20 d vs. 0 d; 60 d vs. 0 d; 15 m vs. 0 d; adult vs. 0 d; 60 d vs. 20 d; 15 m vs. 60 d; and adult vs. 15 m groups). (B) The number of DEmRNAs. (C) The number of common DEmRNAs. (D) KEGG pathway analysis of common DEmRNAs. The top 20 enriched KEGG pathways ranked by p-values are shown. (E) The top 20 expressed mRNAs.
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).
3.5. Expression patterns of DEmRNAs in omasum tissue
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).
Figure 4. Expression patterns of DEmRNAs in omasum. (A) Heatmap of all DEmRNAs among the seven compared groups. (B) The number of DEmRNAs. (C) The number of common DEmRNAs. (D) KEGG pathway analysis of common DEmRNAs. The top 20 enriched KEGG pathways ranked by p-values are shown. (E) The top 20 expressed mRNAs.
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).
3.6. Expression patterns of DEmRNAs in abomasum tissue
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).
Figure 5. Expression patterns of DEmRNAs in abomasum. (A) Heatmap of all DEmRNAs among the seven groups. (B) The number of DEmRNAs. (C) The number of common DEmRNAs. (D) KEGG pathway analysis of common DEmRNAs. The top 20 enriched KEGG pathways ranked by p-values are shown. (E) The top 20 expressed mRNAs.
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).
3.7. DEmRNAs in ruminant stomach
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.
Figure 6. Expression analysis of mRNAs among four stomachs of yak. (A) Venn chart of DEmRNAs detected in rumen, reticulum, omasum and abomasum. (B) The number of DEmRNAs among four stomachs. (C) Heatmap illustrating the relative expression of DEmRNAs from rumen, reticulum, omasum and abomasum tissues at adult. Rows represent mRNAs and columns represent different tissues. (D) KEGG pathway analysis of DEmRNAs among four stomachs tissues. The top 20 enriched KEGG pathways ranked by p-values are presented.
3.8. Time-series analysis (STEM) of mRNAs in rumen, reticulum, omasum, and abomasum
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).
Figure 7. STEM analysis determined temporal expression profiles of mRNAs in rumen (A), reticulum (B), omasum (C), and abomasum (D). The top panel indicates the module number. The number in the upper right corner represents a statistically significant p-value. The number at the bottom left shows the number of mRNAs in each profile module. Different colors indicate the colored profiles (Bonferroni-adjusted p values < 0.05). The p-value is ranked from smallest to largest. Profiles with the same color belong to the same cluster.
3.9. Co-expression network analysis of mRNAs in rumen, reticulum, omasum, and abomasum
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).
Figure 8. The WGCNA network analysis of 59 samples and the expression profiles of the key modules. (A) Hierarchical clustering dendrogram of mRNAs co-expression modules in rumen, reticulum, omasum and abomasum tissues. Each branch shows a cluster of mRNAs. Dynamic tree cut represents original split module, and merged dynamic represents the final merged modules. (B) Module-sample correlation heatmap. The number on the left shows the number of module genes. The color scale on the right represents correlation from −1 (blue) to 1 (red). Day represents developmental stages, and Sex represents male and female yaks. RumW_SW represents rumen weight/total weight of four stomachs, and RumW_W represents rumen weight/body weight. (C) The co-expression network diagram of mRNAs of four tissues in the turquoise module. (D) The co-expression network diagram of mRNAs of four tissues in the blue module. Color indicates different genes, red: hub gene; blue: related genes of hub gene.
3.10. RT-qPCR quantification of mRNAs
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).
4. Discussion
The development of ruminant stomach of yak is closely related to its health and performance, which has attracted wide attention (11, 12, 15). However, the development mechanism of the whole ruminant stomach remains unclear (16), and further studies are needed to provide theoretical references and new perspectives for production in yak. In this study, 3,240, 3,968, 3,975, and 2,638 DEmRNAs were identified in rumen, reticulum, omasum and abomasum, respectively. In these heatmaps, the DEmRNAs expression patterns at 0 d were significantly different from those at the other four stages. These results coincide with those of others, which found that the gene expression patterns of rumen at birth were completely different from those in youth and adulthood in Simmental beef cattle (34). We also found that the DEmRNAs expression patterns at 20 d and 60 d were similar, as well as at 15 m and in adulthood in four stomachs. This is consistent with the rumen status observed at different development stages during dissection. In rumen, we found no solid contents at 0 d, very little pasture and a lot of milk at 20 d, some pasture and a lot of milk at 60 d, and only pasture in the 15 m and adult groups. Indeed, during the suckling and weaning periods, dietary structure changes may also influence the digestive function. This may be related to DEmRNAs (9, 12). Nevertheless, nearly half of the DEmRNAs were down-regulated at 15 m and adult, indicating that the regulatory roles of these DEmRNAs at these stages were relatively poorer than at other three stages in four stomachs. After 15 months, yaks are almost sexually mature and has a mature body with relatively complete organ, which may be the reason why these DEmRNAs are down-regulated.
In this study, we constructed venn diagrams to classify the common DEmRNAs along five developmental stages in rumen, reticulum, omasum and abomasum. In the common DEmRNAs of rumen, we found that HMGCS2, PPARG, ACADS, ECHS1, CBR4, ABHD6, AKR1A1, and INSIG1 were involved in lipid metabolism, CCNB1, KRT6A, IGF2BP2, E2F8, CORO2B, RPS6KA3, CDKN2C, CRISPLD2, CREB3L1, FSTL3, FZD5, and CKAP2L were related to cell cycle, cell differentiation and proliferation, and CX3CR1, CD1E, ARID5A, C3AR1, CD3D, CHGA, FCER2, IL1R2, IL36G, and IL27RA were participated in immunoregulatory and inflammatory processes. HMGCS2, a rate-limiting enzyme that regulates ketone body synthesis, acts an important role in the VFA metabolism of the rumen epithelium (35–37). ACADS, a member of the acyl-CoA dehydrogenase family, catalyzes the initial step of the mitochondrial fatty acid β-oxidation pathway. Mutation of this gene can cause SCAD deficiency an acute acidosis and muscle weakness in infants and lipid-storage myopathy in adults (38, 39). CCNB1, also known as cyclin B1, is involved in the cell division and cell cycle. CCNB1 also plays critical roles in accelerating cell cycle and rumen development in goats (40). CX3CR1, C-X3-C motif chemokine receptor 1, is the only receptor for the chemokine fractalkine (CX3CL1). CX3CL1 belongs to the CX3C subgroup of chemokines and mediates the chemotaxis and adhesion of inflammatory cells through CX3CR1 (41). In the present study, CX3CR1 were observed in four closed groups, further indicating that changes in dietary structure may trigger the persistence of rumen epithelial local inflammation throughout the development stage. Interestingly, these genes associated with fatty acid metabolism, such as HMGCS2, PPARG, and ACADS could also be observed in the common DEmRNAs of reticulum and omasum, but not in abomasum. In addition, according to the correlation analysis of the four stomachs, the expression pattern of abomasum was remarkably different from that of rumen, reticulum and omasum. This also indicated that the four stomachs, especially the rumen and abomasum, had different developmental pathways after birth and subsequent onset of rumination.
From the KEGG pathway analysis of common DEmRNAs in abomasum, among the most abundant mRNAs expressed in abomasum, several mRNAs (e.g., TFF1, TMSB10, GKN1, ATP4B, RPS2, PGC, TFF2, EEF1A1, TPT1, ATP4A, RPS3A, RPLP0, and RPS8) that influence the important roles of abomasum development, gastric acid secretion and digestive enzyme activity are significantly expressed in abomasum. Gastrokine-1 (GKN1), a member of the gastrokine family, is an anti-amyloid secreted by the stomach, has mitogenic activity, actively acts in protecting the integrity of the gastric mucosal epithelium and regulates diet-induced obesity (42). Pepsinogen C (PGC), a member of the aspartic protease family, is secreted by gastric chief cells. PGC can be activated by pepsin C to digest peptides and amino acids in the abomasum (43). Ribosome biogenesis in eukaryotes requires the participation of 80 ribosomal proteins, and its rate must be synchronized with cellular growth. RPS2, RPS3A, RPLP0, and RPS8 each encode a ribosomal protein that regulates the translation and controls cell growth, and may further promote the development of abomasum (44, 45). These findings suggest that these DEmRNAs may be associated with abomasum development, gastric acid secretion and digestive enzyme activity, which may affect the growth of yaks. To further understand the role and regulation of these DEmRNAs, it is necessary to study the functions of the identified DEmRNAs. Better understanding of these DEmRNAs may help producers to manipulate and further improve yak yield.
We realize that most biological processes are dynamic, so sequential experiments (time series) are critical for researchers to gain insight into these processes (46). STEM analysis revealed that mRNAs expression patterns were dynamic and similar during the development of the four stomachs in yaks. In rumen, the STEM analysis showed that mRNAs profiles 39 and 41 belonged to the same cluster enriched in amino acid metabolism and lipid metabolism. In profile 41, short-chain enoyl-CoA hydratase, short chain 1 (ECHS1) was involved in amino and fatty acid catabolism in mitochondria, mainly catalyzing the hydration step of fatty acids, and its deficiency can lead to Leigh syndrome or exercise-induced dystonia (47, 48). ECHS1 had higher expression levels at 20 d, 60 d, 15 m and adult compared to 0 d. The Acyl Coenzyme A oxidase 2 (ACOX2), a member of the acyl-CoA oxidase family, primarily participated in the degradation of long-branched fatty acids and bile acid intermediates in peroxisomes. Lack of this enzyme leads to the accumulation of branched fatty acids and bile acid intermediates (49). Fatty acid beta oxidation is an essential pathway for energy production in growth-stable newborn livestock, and it has been reported that ACOX2 may also be involved in fatty acid β-oxidation (50). We also found that HSD17B2, ACADSB, ECHS1, HADH, ACSS2, ME1, ME3, and SLC27A2 were involved in the lipid metabolism. Additionally, the KEGG pathway analysis indicated that these genes principally participated to lipid-related pathways, such as steroid hormone biosynthesis, fatty acid degradation pyruvate metabolism and PPAR signaling pathway.
In abomasum, STEM analysis showed that mRNA profiles 38, 41, and 43 belonged to the same cluster, which was mainly involved in immune regulation, signaling molecules and interaction, energy metabolism and lipid metabolism. FASLG, also known as CD95L, had higher expression levels at 20 d, 60 d, 15 m and adult compared to 0 d after birth. This may be due to the hypoplasia of the gastrointestinal tract, lower immunity and lower expression level of related genes in newborn calves (0 d). FASLG has been reported to be critical in triggering apoptosis of some types of cells such as lymphocytes and defects in this gene may be associated with lymphoproliferation and thus a predisposition to autoimmunity (51). We also found that ITK, ITGAL, PIK3CD, SH2D1A, CXCR3 AICDA, IL15, and ZAP70 were involved in the regulation of the immune system. KEGG pathway also indicated that these genes were mainly involved in immune-related pathways, such as natural killer cell mediated cytotoxicity, Th17 cell differentiation and intestinal immune network for IgA production.
By the rumen-reticulum-omasum-abomasum-mRNA co-expression network analysis, CCKBR, KCNQ1, FER1L6, and A4GNT were determined as the hub genes of turquoise module, while PAK6, TRIM29, ADGRF4, TGM1, and TMEM79 were identified as the hub genes of the blue module. The cholecystokinin B receptor (CCKBR), mainly found in the central nervous system and gastrointestinal tract, has a high affinity for both sulfated and nonsulfated CCK analogs (52). Previous studies have indicated that the genetic inactivation of CCKBR causes deficits in the gastrointestinal system, control of food intake, memory and exploration, and anxiety-related behaviors (53–55). However, the functional characterization of CCKBR remains scarce in yaks. A4GNT encodes α-1,4-N-acetylglucosaminyltransferase, which transfers N-acetylglucosamine (GlcNAc) to core 2 branched O-glycans, forming a unique structure GlcNAcα1-4Galβ-R (56). In addition, A4GNT plays a vital role in the biosynthesis of type III mucins. Furthermore, the knockdown of A4GNT significantly influences gastric development, gradually leading to gastric dysplasia, precancerous lesion of gastric cancer, dysplasia (57, 58). Since A4GNT is a hub gene in the co-expression analysis, it is speculated that A4GNT may play an essential regulatory role in maintaining normal metabolic homeostasis in yak stomachs. The p21 protein (Cdc42/Rac)-activated kinase 6 (PAK6) regulates a variety of biological activities, including cell proliferation, apoptosis, invasion, metastasis, cytoskeleton rearrangement and the MAP kinase signaling pathway. It has been found that PAK6 interacts with androgen receptor (AR) and estrogen receptor (ER), two steroid hormone-dependent transcription factors that play important roles in sexual differentiation and development (59). However, the functional characterization of PAK6 remains rare in yak stomachs. The functional and regulatory role of PAK6 in the growth and development of yak stomachs will be further studied.
In addition to hub genes, the gastric acid secretion, hedgehog (Hh) signaling pathway, sphingolipid metabolism, glycine, serine and threonine metabolism, carbohydrate digestion and absorption, mineral absorption, thyroid hormone synthesis and cAMP signaling pathway were also enriched and involved in the development and metabolism of yak stomachs. Hh signaling pathway, a highly conserved signaling pathway, participated in regulating the invertebrate and vertebrate organogenesis and tissue homeostasis such as the gastrointestinal tract (60, 61). Previous studies reveal that Hh signaling may inhibit fat formation and influence the adipogenic differentiation of preadipocytes (62, 63). Hedgehog signaling pathway has also been declared to promote lipolysis in adipose tissue by directly regulating brummer (Bmm)/adipose triglyceride lipase (ATGL) (64). In this study, abundant genes were enriched in the Hh signaling pathway. Consequently, it is speculated that the Hh signaling pathway plays a key role in regulating the lipid metabolism in ruminant stomach, which may provide a reference for improving fat ratio and weight gain of yaks. Thyroid hormone (TH) is an endocrine messenger that regulates metabolic processes necessary for normal growth, development, and function in almost all vertebrates (65–67). It is well determined that TH status is closely related to body weight and energy expenditure (68–70). Previous studies have reported that TH can stimulate both lipogenesis and lipolysis, as well as excess thyroid hormones, promoting a hypermetabolic state characterized by increased resting energy expenditure, elevated cholesterol levels, reduced lipolysis, and gluconeogenesis, etc. (71, 72). Clinically, thyroid dysfunction manifests as dominant or subclinical hypothyroidism that negatively affects lipid metabolism and leads to hypercholesterolemia, which gradually increases the risk of cardiovascular disease and potential mortality (73). The role of TH in energy metabolism and lipid metabolism may be helpful to improve the survival condition of yaks in low oxygen, high cold and high altitude weather. cAMP signaling pathway, a classic pathway, participated in the regulation of key physiological processes, including gene transcription, cell fate, muscle contraction, metabolism, secretion and calcium homeostasis (74–76). cAMP signaling pathway has been declared to control adipocyte differentiation and regulate lipid metabolism in white adipose tissue (77). Previous studies have shown that cAMP signaling pathway can regulate the feed efficiency of yorkshire pigs by participating in the influence of lipid metabolism in adipose tissue. These results show that the hedgehog signaling pathway, thyroid hormone synthesis and cAMP signaling pathway may play a crucial role in maintaining normal homeostasis and growth metabolism in yaks. The ruminant stomachs are the major digestive organs participated in nutrient digestion and absorption in yaks. Living at high altitudes and in cold environments all year round, yaks require normal homeostasis to maintain proper growth. Hence, a deeper understanding of mRNAs co-expression in rumen, reticulum, omasum and abomasum can provide some clues and basic knowledge for improving homeostasis, feed utilization and meat quality, and thus improve the growth and survival of yaks.
5. Conclusion
In summary, we first revealed the expression profile of mRNAs in the rumen, reticulum, omasum and abomasum of yaks in five developmental stages. In this study, 3,240, 3,968, 3,975, and 2,638 DEmRNAs were identified in rumen, reticulum, omasum and abomasum, respectively. Overall, the expression patterns of DEmRNAs were unique at 0 d, similar at 20 d and 60 d, and similar at 15 m and adult in four stomachs. 20 days is an important node for regulating immune response, metabolism-related hormone levels, and enzymatic activity in four stomachs tissues. On the whole, the expression pattern of abomasum was greatly different from that of rumen, reticulum and omasum. STEM and WGCNA analysis demonstrated that many important mRNAs were involved in lipids, protein metabolism, saccharides, and other biological processes. They may perform key functions in health, growth and development, and tissue repair of ruminant stomach. The potential regulatory relationship of mRNAs between the four stomachs has been further explored and confirmed. These findings can provide data reference for further study of the function and mechanisms of important mRNAs in the four stomachs of yaks. It also provided an important theoretical basis for the research of ruminant stomach in other yak breeds. This study provides a molecular basis for further research on age-appropriate weaning and supplementary feeding of yak calves, and also provides an important reference for the research on nutrition and feeding management of yaks. This will help improve our understanding of ruminant physiology and facilitate further development of animal husbandry, while reducing overgrazing pressures on the natural grasslands on the Qinghai-Tibetan Plateau.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE222396.
Ethics statement
The animal studies were approved by The Institutional Animal Care and Use Committee at Southwest Minzu University. The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent was obtained from the owners for the participation of their animals in this study.
Author contributions
YL performed data analysis and drafted the manuscript. QM conducted the experiments and data analysis, with equal contributions. JT and LY participated in the sorting of the test data. TP and XM participated in partial operation of the RT-PCR experiments. MJ designed this study and revised the manuscript. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by the Natural Science Foundation of Sichuan Province (2022NSFSC1665), Southwest Minzu University Research Startup Funds (RQD2022048), Sichuan Science and Technology Program (2023YFQ0076), Sichuan Science and Technology Program (2021YFYZ0001 and 2021YFN0001).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fvets.2023.1204706/full#supplementary-material
Footnotes
References
1. Guo, X, Long, R, Kreuzer, M, Ding, L, Shang, Z, Zhang, Y, et al. Importance of functional ingredients in yak milk-derived food on health of Tibetan nomads living under high-altitude stress: a review. Crit Rev Food Sci Nutr. (2014) 54:292–302. doi: 10.1080/10408398.2011.584134
2. Huang, Q, Dong, K, Wang, Q, Huang, X, Wang, G, An, F, et al. Changes in volatile flavor of yak meat during oxidation based on multi-omics. Food Chem. (2022) 371:131103. doi: 10.1016/j.foodchem.2021.131103
3. Long, R, Dong, S, Wei, X, and Pu, X. The effect of supplementary feeds on the bodyweight of yaks in cold season. Livest Prod Sci. (2005) 93:197–204. doi: 10.1016/j.livprodsci.2004.08.016
4. Grünberg, W, and Constable, PD. Function and dysfunction of the ruminant forestomach. Current veterinary. Therapy. (2009):12–9. doi: 10.1016/B978-141603591-6.10006-5
5. Baaske, L, Gäbel, G, and Dengler, F. Ruminal epithelium: a checkpoint for cattle health. J Dairy Res. (2020) 87:322–9. doi: 10.1017/S0022029920000369
6. Aschenbach, JR, Penner, GB, Stumpff, F, and Gäbel, G. Ruminant nutrition symposium: role of fermentation acid absorption in the regulation of ruminal pH. J Anim Sci. (2011) 89:1092–107. doi: 10.2527/jas.2010-3301
7. Malmuthuge, N, Liang, G, and Guan, LL. Regulation of rumen development in neonatal ruminants through microbial metagenomes and host transcriptomes. Genome Biol. (2019) 20:172–16. doi: 10.1186/s13059-019-1786-0
8. Xue, Y, Yin, Y, Trabi, EB, Xie, F, Lin, L, and Mao, S. Transcriptome analysis reveals the effect of high-grain pelleted and non-pelleted diets on ruminal epithelium of Hu-lamb. Animal. (2021) 15:100278. doi: 10.1016/j.animal.2021.100278
9. Ramos, SC, Jeong, CD, Mamuad, LL, Kim, SH, Kang, SH, Kim, ET, et al. Diet transition from high-forage to high-concentrate alters rumen bacterial community composition, epithelial transcriptomes and ruminal fermentation parameters in dairy cows. Animals. (2021) 11:838. doi: 10.3390/ani11030838
10. Li, W, Gelsinger, S, Edwards, A, Riehle, C, and Koch, D. Transcriptome analysis of rumen epithelium and meta-transcriptome analysis of rumen epimural microbial community in young calves with feed induced acidosis. Sci Rep. (2019) 9:4744. doi: 10.1038/s41598-019-40375-2
11. Kong, RS, Liang, G, Chen, Y, Stothard, P, and Guan le, L. Transcriptome profiling of the rumen epithelium of beef cattle differing in residual feed intake. BMC Genomics. (2016) 17:592. doi: 10.1186/s12864-016-2935-4
12. Wang, B, Wang, D, Wu, X, Cai, J, Liu, M, Huang, X, et al. Effects of dietary physical or nutritional factors on morphology of rumen papillae and transcriptome changes in lactating dairy cows based on three different forage-based diets. BMC Genomics. (2017) 18:353. doi: 10.1186/s12864-017-3726-2
13. Wang, X, Zhang, D, Wang, W, Lv, F, Pang, X, Liu, G, et al. Transcriptome profiling reveals differential gene expression in the rumen of Hu lambs at different developmental stages. Anim Biotechnol. (2021) 34:471–81. doi: 10.1080/10495398.2021.1975728
14. Xin, J, Chai, Z, Zhang, C, Zhang, Q, Zhu, Y, Cao, H, et al. Comparing the microbial Community in Four Stomach of dairy cattle, yellow cattle and three yak herds in Qinghai-Tibetan plateau. Front Microbiol. (2019) 10:1547. doi: 10.3389/fmicb.2019.01547
15. Xue, D, Chen, H, Luo, X, Guan, J, He, Y, and Zhao, X. Microbial diversity in the rumen, reticulum, omasum, and abomasum of yak on a rapid fattening regime in an agro-pastoral transition zone. J Microbiol. (2018) 56:734–43. doi: 10.1007/s12275-018-8133-0
16. Roh, SG, Kuno, M, Hishikawa, D, Hong, YH, Katoh, K, Obara, Y, et al. Identification of differentially expressed transcripts in bovine rumen and abomasum using a differential display method. J Anim Sci. (2007) 85:395–403. doi: 10.2527/jas.2006-234
17. Bolger, AM, Lohse, M, and Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. (2014) 30:2114–20. doi: 10.1093/bioinformatics/btu170
18. Kim, D, Langmead, B, and Salzberg, SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. (2015) 12:357–60. doi: 10.1038/nmeth.3317
19. Roberts, A, Trapnell, C, Donaghey, J, Rinn, JL, and Pachter, L. Improving RNA-Seq expression estimates by correcting for fragment bias. Genome Biol. (2011) 12:R22. doi: 10.1186/gb-2011-12-3-r22
20. Trapnell, C, Williams, BA, Pertea, G, Mortazavi, A, Kwan, G, van Baren, MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. (2010) 28:511–5. doi: 10.1038/nbt.1621
21. Anders, S, Pyl, PT, and Huber, W. HTSeq--a Python framework to work with high-throughput sequencing data. Bioinformatics. (2015) 31:166–9. doi: 10.1093/bioinformatics/btu638
22. Langmead, B, and Salzberg, SL. Fast gapped-read alignment with bowtie 2. Nat Methods. (2012) 9:357–9. doi: 10.1038/nmeth.1923
23. Roberts, A, and Pachter, L. Streaming fragment assignment for real-time analysis of sequencing experiments. Nat Methods. (2013) 10:71–3. doi: 10.1038/nmeth.2251
24. Pertea, M, Pertea, GM, Antonescu, CM, Chang, TC, Mendell, JT, and Salzberg, SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. (2015) 33:290–5. doi: 10.1038/nbt.3122
25. Florea, L, Song, L, and Salzberg, SL. Thousands of exon skipping events differentiate among splicing patterns in sixteen human tissues. F1000Res. (2013) 2:188. doi: 10.12688/f1000research.2-188.v2
26. Li, H, Handsaker, B, Wysoker, A, Fennell, T, Ruan, J, Homer, N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. (2009) 25:2078–9. doi: 10.1093/bioinformatics/btp352
27. Li, H . A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics. (2011) 27:2987–93. doi: 10.1093/bioinformatics/btr509
28. Cingolani, P, Platts, A, Wang, LL, Coon, M, Nguyen, T, Wang, L, et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly (Austin). (2012) 6:80–92. doi: 10.4161/fly.19695
29. Anders, S, and Huber, W. Differential expression of RNA-Seq data at the gene level–the DESeq package, vol. 10. Heidelberg: European Molecular Biology Laboratory (EMBL) (2012). f1000research p.
30. Kanehisa, M, Araki, M, Goto, S, Hattori, M, Hirakawa, M, Itoh, M, et al. KEGG for linking genomes to life and the environment. Nucleic Acids Res. (2008) 36:D480–4. doi: 10.1093/nar/gkm882
31. Ernst, J, and Bar-Joseph, Z. STEM: a tool for the analysis of short time series gene expression data. BMC Bioinformatics. (2006) 7:191. doi: 10.1186/1471-2105-7-191
32. Langfelder, P, and Horvath, S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. (2008) 9:559. doi: 10.1186/1471-2105-9-559
33. Langfelder, P, Zhang, B, and Horvath, S. Defining clusters from a hierarchical cluster tree: the dynamic tree cut package for R. Bioinformatics. (2008) 24:719–20. doi: 10.1093/bioinformatics/btm563
34. Zhang, Y, Cai, W, Li, Q, Wang, Y, Wang, Z, Zhang, Q, et al. Transcriptome analysis of bovine rumen tissue in three developmental stages. Front Genet. (2022) 13:821406. doi: 10.3389/fgene.2022.821406
35. Nakamura, MT, Yudell, BE, and Loor, JJ. Regulation of energy metabolism by long-chain fatty acids. Prog Lipid Res. (2014) 53:124–44. doi: 10.1016/j.plipres.2013.12.001
36. Mirzaei-Alamouti, H, Moradi, S, Patra, AK, and Mansouryar, M. Monensin supplementation downregulated the expression signature of genes involved in cholesterol synthesis in the ruminal epithelium and adipose tissue of lambs. Trop Anim Health Prod. (2022) 54:167. doi: 10.1007/s11250-022-03168-w
37. Wang, W, Li, C, Li, F, Wang, X, Zhang, X, Liu, T, et al. Effects of early feeding on the host rumen transcriptome and bacterial diversity in lambs. Sci Rep. (2016) 6:32479. doi: 10.1038/srep32479
38. He, D, Cai, L, Huang, W, Weng, Q, Lin, X, You, M, et al. Prognostic value of fatty acid metabolism-related genes in patients with hepatocellular carcinoma. Aging (Albany NY). (2021) 13:17847–63. doi: 10.18632/aging.203288
39. He, M, Pei, Z, Mohsen, AW, Watkins, P, Murdoch, G, van Veldhoven, PP, et al. Identification and characterization of new long chain acyl-CoA dehydrogenases. Mol Genet Metab. (2011) 102:418–29. doi: 10.1016/j.ymgme.2010.12.005
40. Gui, H, and Shen, Z. Concentrate diet modulation of ruminal genes involved in cell proliferation and apoptosis is related to combined effects of short-chain fatty acid and pH in rumen of goats. J Dairy Sci. (2016) 99:6627–38. doi: 10.3168/jds.2015-10446
41. Verge, GM, Milligan, ED, Maier, SF, Watkins, LR, Naeve, GS, and Foster, AC. Fractalkine (CX3CL1) and fractalkine receptor (CX3CR1) distribution in spinal cord and dorsal root ganglia under basal and neuropathic pain conditions. Eur J Neurosci. (2004) 20:1150–60. doi: 10.1111/j.1460-9568.2004.03593.x
42. Overstreet, AC, Grayson, BE, Boger, A, Bakke, D, Carmody, EM, Bales, CE, et al. Gastrokine-1, an anti-amyloidogenic protein secreted by the stomach, regulates diet-induced obesity. Sci Rep. (2021) 11:9477. doi: 10.1038/s41598-021-88928-8
43. Shen, S, Jiang, J, and Yuan, Y. Pepsinogen C expression, regulation and its relationship with cancer. Cancer Cell Int. (2017) 17:57. doi: 10.1186/s12935-017-0426-6
44. Jiang, X, Prabhakar, A, van der Voorn, SM, Ghatpande, P, Celona, B, Venkataramanan, S, et al. Control of ribosomal protein synthesis by the microprocessor complex. Sci Signal. (2021) 14. doi: 10.1126/scisignal.abd2639
45. Ishii, K, Washio, T, Uechi, T, Yoshihama, M, Kenmochi, N, and Tomita, M. Characteristics and clustering of human ribosomal protein genes. BMC Genomics. (2006) 7:37. doi: 10.1186/1471-2164-7-37
46. Bar-Joseph, Z, Gitter, A, and Simon, I. Studying and modelling dynamic biological processes using time-series gene expression data. Nat Rev Genet. (2012) 13:552–64. doi: 10.1038/nrg3244
47. Schäff, C, Börner, S, Hacke, S, Kautzsch, U, Albrecht, D, Hammon, HM, et al. Increased anaplerosis, TCA cycling, and oxidative phosphorylation in the liver of dairy cows with intensive body fat mobilization during early lactation. J Proteome Res. (2012) 11:5503–14. doi: 10.1021/pr300732n
48. Kuwajima, M, Kojima, K, Osaka, H, Hamada, Y, Jimbo, E, Watanabe, M, et al. Valine metabolites analysis in ECHS1 deficiency. Mol Genet Metab Rep. (2021) 29:100809. doi: 10.1016/j.ymgmr.2021.100809
49. Zhang, Y, Chen, Y, Zhang, Z, Tao, X, Xu, S, Zhang, X, et al. Acox2 is a regulator of lysine crotonylation that mediates hepatic metabolic homeostasis in mice. Cell Death Dis. (2022) 13:279. doi: 10.1038/s41419-022-04725-9
50. Yeh, CS, Wang, JY, Cheng, TL, Juan, CH, Wu, CH, and Lin, SR. Fatty acid metabolism pathway play an important role in carcinogenesis of human colorectal cancers by microarray-bioinformatics analysis. Cancer Lett. (2006) 233:297–308. doi: 10.1016/j.canlet.2005.03.050
51. Poissonnier, A, Sanséau, D, le Gallo, M, Malleter, M, Levoin, N, Viel, R, et al. CD95-mediated calcium signaling promotes T helper 17 trafficking to inflamed organs in lupus-prone mice. Immunity. (2016) 45:209–23. doi: 10.1016/j.immuni.2016.06.028
52. Nishimura, S, Bilgüvar, K, Ishigame, K, Sestan, N, Günel, M, and Louvi, A. Functional synergy between cholecystokinin receptors CCKAR and CCKBR in mammalian brain development. PLoS One. (2015) 10:e0124295. doi: 10.1371/journal.pone.0124295
53. Nagata, A, Ito, M, Iwata, N, Kuno, J, Takano, H, Minowa, O, et al. G protein-coupled cholecystokinin-B/gastrin receptors are responsible for physiological cell growth of the stomach mucosa in vivo. Proc Natl Acad Sci USA. (1996) 93:11825–30. doi: 10.1073/pnas.93.21.11825
54. Langhans, N, Rindi, G, Chiu, M, Rehfeld, JF, Ardman, B, Beinborn, M, et al. Abnormal gastric histology and decreased acid production in cholecystokinin-B/gastrin receptor-deficient mice. Gastroenterology. (1997) 112:280–6. doi: 10.1016/s0016-5085(97)90000-7
55. Miyasaka, K, Ohta, M, Kanai, S, Yoshida, Y, Sato, N, Nagata, A, et al. Enhanced gastric emptying of a liquid gastric load in mice lacking cholecystokinin-B receptor: a study of CCK-A,B, and AB receptor gene knockout mice. J Gastroenterol. (2004) 39:319–23. doi: 10.1007/s00535-003-1297-2
56. Zhang, MX, Nakayama, J, Hidaka, E, Kubota, S, Yan, J, Ota, H, et al. Immunohistochemical demonstration of alpha1,4-N-acetylglucosaminyltransferase that forms GlcNAcalpha1,4Galbeta residues in human gastrointestinal mucosa. J Histochem Cytochem. (2001) 49:587–96. doi: 10.1177/002215540104900505
57. Iida, M, Desamero, MJ, Yasuda, K, Nakashima, A, Suzuki, K, Chambers, JK, et al. Effects of orally administered Euglena gracilis and its reserve polysaccharide, paramylon, on gastric dysplasia in A4gnt knockout mice. Sci Rep. (2021) 11:13640. doi: 10.1038/s41598-021-92013-5
58. Desamero, MJ, Kakuta, S, Chambers, JK, Uchida, K, Hachimura, S, Takamoto, M, et al. Orally administered brown seaweed-derived β-glucan effectively restrained development of gastric dysplasia in A4gnt KO mice that spontaneously develop gastric adenocarcinoma. Int Immunopharmacol. (2018) 60:211–20. doi: 10.1016/j.intimp.2018.05.002
59. Lee, SR, Ramos, SM, Ko, A, Masiello, D, Swanson, KD, Lu, ML, et al. AR and ER interaction with a p21-activated kinase (PAK6). Mol Endocrinol. (2002) 16:85–99. doi: 10.1210/mend.16.1.0753
60. Kong, JH, Siebold, C, and Rohatgi, R. Biochemical mechanisms of vertebrate hedgehog signaling. Development. (2019) 146. doi: 10.1242/dev.166892
61. van den Brink, GR . Hedgehog signaling in development and homeostasis of the gastrointestinal tract. Physiol Rev. (2007) 87:1343–75. doi: 10.1152/physrev.00054.2006
62. Suh, JM, Gao, X, McKay, J, McKay, R, Salo, Z, and Graff, JM. Hedgehog signaling plays a conserved role in inhibiting fat formation. Cell Metab. (2006) 3:25–34. doi: 10.1016/j.cmet.2005.11.012
63. James, AW, Leucht, P, Levi, B, Carre, AL, Xu, Y, Helms, JA, et al. Sonic hedgehog influences the balance of osteogenesis and adipogenesis in mouse adipose-derived stromal cells. Tissue Eng Part A. (2010) 16:2605–16. doi: 10.1089/ten
64. Zhang, J, Liu, Y, Jiang, K, and Jia, J. Hedgehog signaling promotes lipolysis in adipose tissue through directly regulating Bmm/ATGL lipase. Dev Biol. (2020) 457:128–39. doi: 10.1016/j.ydbio.2019.09.009
65. Brent, GA . Mechanisms of thyroid hormone action. J Clin Invest. (2012) 122:3035–43. doi: 10.1172/JCI60047
66. Cheng, SY, Leonard, JL, and Davis, PJ. Molecular aspects of thyroid hormone actions. Endocr Rev. (2010) 31:139–70. doi: 10.1210/er.2009-0007
67. Mendoza, A, and Hollenberg, AN. New insights into thyroid hormone action. Pharmacol Ther. (2017) 173:135–45. doi: 10.1016/j.pharmthera.2017.02.012
68. Fox, CS, Pencina, MJ, D'Agostino, RB, Murabito, JM, Seely, EW, Pearce, EN, et al. Relations of thyroid function to body weight: cross-sectional and longitudinal observations in a community-based sample. Arch Intern Med. (2008) 168:587–92. doi: 10.1001/archinte.168.6.587
69. Yongqiang, T, Xingxu, Z, Minqiang, W, Zhonglin, L, and Rongchang, Z. Endocrine changes and their relationship with body weight in growing yaks. Anim Sci. (2002) 74:89–94. doi: 10.1017/S1357729800052243
70. Knudsen, N, Laurberg, P, Rasmussen, LB, Bülow, I, Perrild, H, Ovesen, L, et al. Small differences in thyroid function may be important for body mass index and the occurrence of obesity in the population. J Clin Endocrinol Metab. (2005) 90:4019–24. doi: 10.1210/jc.2004-2225
71. Motomura, K, and Brent, GA. Mechanisms of thyroid hormone action. Implications for the clinical manifestation of thyrotoxicosis. Endocrinol Metab Clin N Am. (1998) 27:1–23. doi: 10.1016/s0889-8529(05)70294-2
72. Baranowska-Bik, A, and Bik, W. The Association of Obesity with autoimmune thyroiditis and thyroid function-possible mechanisms of bilateral interaction. Int J Endocrinol. (2020) 2020:8894792–14. doi: 10.1155/2020/8894792
73. Duntas, LH, and Brenta, G. A renewed focus on the association between thyroid hormones and lipid metabolism. Front Endocrinol (Lausanne). (2018) 9:511. doi: 10.3389/fendo.2018.00511
74. Schmidt, M, Dekker, FJ, and Maarsingh, H. Exchange protein directly activated by cAMP (epac): a multidomain cAMP mediator in the regulation of diverse biological functions. Pharmacol Rev. (2013) 65:670–709. doi: 10.1124/pr.110.003707
75. Emamghoreishi, M, Li, PP, Schlichter, L, Parikh, SV, Cooke, R, and Warsh, JJ. Associated disturbances in calcium homeostasis and G protein-mediated cAMP signaling in bipolar I disorder. Biol Psychiatry. (2000) 48:665–73. doi: 10.1016/s0006-3223(00)00884-2
76. Freitas, FZ, de Paula, RM, Barbosa, LC, Terenzi, HF, and Bertolini, MC. cAMP signaling pathway controls glycogen metabolism in Neurospora crassa by regulating the glycogen synthase gene expression and phosphorylation. Fungal Genet Biol. (2010) 47:43–52. doi: 10.1016/j.fgb.2009.10.011
Keywords: yak, stomach, DEmRNAs, STEM, WGCNA
Citation: Liu Y, Min Q, Tang J, Yang L, Meng X, Peng T and Jiang M (2023) Transcriptome profiling in rumen, reticulum, omasum, and abomasum tissues during the developmental transition of pre-ruminant to the ruminant in yaks. Front. Vet. Sci. 10:1204706. doi: 10.3389/fvets.2023.1204706
Edited by:
Mohammad Reza Bakhtiarizadeh, University of Tehran, IranReviewed by:
Qianjun Zhao, Chinese Academy of Agricultural Sciences, ChinaYangchun Cao, Northwest A&F University, China
Copyright © 2023 Liu, Min, Tang, Yang, Meng, Peng and Jiang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Mingfeng Jiang, bWluZ2ZlbmdqaWFuZ0B2aXAuc2luYS5jb20=
†These authors have contributed equally to this work