Weighted Gene Co-expression Network Analysis Identifies Specific Modules and Hub Genes Related to Subacute Ruminal Acidosis

Weighted gene co-expression network analysis (WGCNA) was used to understand the pathogenesis of subacute ruminal acidosis (SARA) and identify potential genes related to the disease. Microarray data from dataset GSE143765, which included 22 cows with and nine cows without SARA, were downloaded from the NCBI Gene Expression Omnibus (GEO). Results of WGCNA identified highly correlated modules of sample genes, and Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses allowed further biological insights into SARA-related modules. The protein-protein interaction (PPI) network, modules from the PPI network, and cistron annotation enrichment of modules were also analyzed. A total of 14,590 DEGs were used for the WGCNA. Construction of a protein-protein network identified DCXR, MMP15, and MMP17 as hub genes. Functional annotation showed that DCXR mainly exhibited L-xylulose reductase (NADP+) activity, glucose metabolic process, xylulose metabolic process, and carbonyl reductase (NADPH) activity, which are involved in the pentose and glucuronate interconversion pathways. MMP15 and MMP17 mainly have a collagen catabolic process. Overall, the results of this study aid the clarification of the biological and metabolic processes associated with SARA at the molecular level. The data highlight potential mechanisms for the future development of intervention strategies to reduce or alleviate the risk of SARA.

Introduction

Subacute ruminal acidosis (SARA) is the most common nutritional disorder afflicting dairy cows. SARA occurs when the rumen pH is between 5.0 and 5.5 for a period of 3 h or more each day (1). The onset of SARA is associated with increased concentrations of short-chain fatty acids (SCFA), such as propionic acid, acetic acid, and butyric acid, created during feed fermentation. Particularly, the quantities of these SCFA exceed the absorption capacity of the rumen wall. Under these conditions, the large accumulation of SCFA in the rumen causes a gradual pH decrease. Hence, when SARA occurs, the chronic low pH environment causes a large number of microorganisms in the rumen to lyse and die, releasing a large amount of lipopolysaccharide (LPS), and damaging the rumen epithelial tissue (2, 3).

In intensive dairy production, the incidence of SARA ranges from 11 to 26%, and can reach >40% in some herds (4). According to estimates, the economic loss due to reduced milk production alone amounts to $400 per cow per one lactation (5, 6). Despite the grave economic losses caused by SARA in the international dairy industry, there is little scientific understanding of the pathophysiology that triggers SARA. Additionally, current research theories about using bioinformatics technology to explore the mechanism of SARA are lacking (7). The use of microarray technology in animal disease research has increased steadily, and bioinformatics data processing of gene expression profiles is widely used to discover potential new diagnostic and therapeutic targets that could help prevent disease (8).

In the past decades, SARA has been extensively studied, and some potential new therapeutic targets have been identified (9), including metabolic pathways of SCFA (10), immune cell suppression and inflammatory response (11), and lipocalin through the blood as preliminary diagnostic indicators of SARA (12). Particularly, the inflammatory response and the immune system have been recognized to have important effects on rumen epithelial health in lactating high-yielding dairy cows, in particular, those related to the systemic immune response associated with SARA (13, 14). Hence, it is important to better understand the genes and pathways that may be associated with systemic immunity and inflammation to further understand the pathophysiology of SARA (15).

Recent genome-wide association studies (GWAS), along with classical studies on metabolism in dairy cows, have increased our understanding of SARA. However, the molecular mechanisms underlying the regulation of the metabolic responses remain unclear (16). The weighted gene co-expression network analysis (WGCNA) algorithmic program clusters sequences into modules as a function of co-expression similarities across samples, leading to a cluster of genes with similar functions. The modules can then be related to one another, to whole-animal phenotypes, and help determine tissue-specific biomarkers and pathophysiological pathways (17, 18).

In the current study, we created correlation networks using publicly accessible resources based on the microarray dataset GSE143765 from the Gene Expression Omnibus (GEO). This study aimed to construct a sequence co-expression network to predict clusters of candidate genes concerned within pathological process associated with SARA. We screened for differentially expressed genes and constructed co-expression networks for all genes in the sample using WGCNA. Subsequently, cistron modules associated with SARA were identified. We used gene ontology (GO) and the genes and genomes (KEGG) pathway enrichment analyses to further uncover biological insights from modules that were highly correlated. Modules of protein-protein interaction (PPI) networks associated with SARA were also screened. In the present study, all potential genes were analyzed to ensure that our results were complete and reliable. The results of this study may facilitate elucidation of the pathophysiological features of SARA development at the molecular level. Potential molecular targets for the development of interventions to prevent or reduce the risk of SARA could be identified.

Methods

Collection of Microarray Data

The mRNA microarray expression profile dataset GSE143765 was downloaded from the gene expression omnibus (GEO) database. Data from the liver samples from 22 cows with SARA 9 cows without SARA were obtained from a study by Kizaki et al. (19) (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE143765). The GEO database is a publicly accessible functional genomic database that contains high-throughput data, including microarrays. The GSE143765 dataset (GPL22092 [Agilent-072598]) was generated using the Bovine_cusotm8x15K platform (Agilent Probe version). We downloaded the raw TXT and the probe annotation files. The probes were assigned the corresponding gene symbols in steps using the annotation data on the platform. All information is freely available on the manufacturer's website. The present study did not include any experiments involving humans or animals.

Data Preprocessing and DEG Screening

After downloading the raw data in TXT format, the Affymetrix package (20) (http://www.bioconductor.org/packages/release/bioc/html/affy.html) in R software (version 4.1.1, https://www.r-project.org/) was used for data preprocessing, and a robust multiarray average (RMA) was obtained after removing batch differences and performing data background correction, normalization, and summarization. To characterize DEGs, SARA and non-SARA groups were analyzed using the LIMMA (linear models for microarray data) package (21) within the R/Bioconductor platform. Benjamini–Hochberg's false discovery rate was used to reduce the likelihood of Type 2 errors and to adjust P-values < 0.05.

Construction of Weighted Gene Co-expression Networks

As systems biology methodology, the development of cistron co-expression networks and the identification of cistron clusters or modules are particularly helpful in characterizing transcriptional alterations in multigene diseases. This is particularly true for phenotypes created by the convergence of small changes in gene expression rather than isolated single-gene effects (22, 23). The weighted factor co-expression network analysis (WGCNA) package (Version one.70-3, https://cran.rstudio.com/web/packages/WGCNA/index.html) in the R software package was used to construct groups of powerfully co-expressed genes into co-expression networks. WGCNA of all expressed genes (14,590 genes). The selection of soft-threshold power is a vital step during the construction of a WGCNA (24). Data preprocessing uses the WGCNA built-in goodSamplesGenes function to check missing items in gene data, items with weights below the threshold, and zero-variance genes. During module choice by cluster analysis, the contiguity matrix was used to calculate the topological overlap measure (TOM) and strength between all sequence pairs (25, 26). Modules of co-expressed genes were generated using stratified cluster dendrograms with completely different colors. The module structure was generated using topological overlap matrix plots. Finally, to ensure the reliability of the results of each module, the minimum number of genes in the module was set to five and the similarity module merging distance was 0.8. This allowed the identification of candidate genes with the highest correlation with SARA for further analysis and validation.

Functional Enrichment Analysis of Genes in Key Modules of Cows With SARA

Gene modules with a P-value <0.05 and the highest correlation with SARA were selected. Analysis of cellular component (CC), molecular function (MF), and potential biological processes (BP) was performed to identify and interpret advanced biological functions supported by the GO terms and KEGG pathway annotation within the co-expression modules. Genes of every handpicked module were submitted to the web-based software DAVID (https://david.ncifcrf.gov/conversion.jsp) for practical and pathway enrichment analyses. DAVID is a helpful online tool for gene expansion (27–29), which provides the practicality of performing concurrent GO and KEGG pathway analyses. P-values <0.05 were used to uncover vital variations.

Integration of Genetics and Highly Connected Hubs in Modules

The top-ranked genes within the modules were considered as hub genes. To consistently analyze hub genes in modules and module eigengenes, genes obtained from every module were mapped using the net search tool STRING information (30) (STRING, V11.5; https://string-db.org/), which can play a key role in identifying protein-protein networks (PPI). A combined score ≥0.4 of PPI pairs was defined as vital. The CytoHubba plugin-supported Cytoscape package (31) (http://www.cytoscape.org/version3.7.1; Institute for Systems Biology, Seattle, WA, USA) was used to construct and visualize the transcriptional regulatory network modules.

Results

Identification of DEGs Associated With Non-SARA and SARA Samples

Thirty-one raw tissue sample files (TXT format) were downloaded from the GEO information database. After batch normalization using the SVA and LIMMA packages in the R language, totally 15,068 probes on the comprehensive dataset GSE143765 were extracted with 14,646 gene expressions.

Weighted Co-expression Network Construction

In this study, soft threshold power was screened once the degree of scale independence was set to β = 7 (scale-free R2 = 0.95). When the degree of scale independence was set as β = 7, the measures decreased steeply with increasing soft threshold power, and subsequently, cistron co-expression module similarity and closeness matrices were obtained using the WGCNA algorithmic rule. After pretreatment, the number of genes was filtered from 14,590 to 2,464, we used WGCNA to spot the modules containing extremely correlated genes (Figures 1A,B). We set MEDissThres to 0.8 to merge similar modules, leading to 23 practical modules (Figure 2). There were 15 genes within the black module, 338 genes within the blue module, 228 genes within the brown module, 9 genes within the cyan module, 5 genes in dark green, dark red, dark-turquoise, and royal blue modules, 73 genes within the green module, 12 genes within the green-yellow module, 7 genes within the gray 60 module, 8 genes within the light-cyan module, 6 genes within the light-green module, 6 genes within the light-yellow module, 14 genes within the light-yellow module, 9 genes within the midnight-blue module, 15 genes within the pink module, 12 genes within the purple module, 38 genes within the red module,11 genes within the salmon module, 11 genes within the tan module, 558 genes within the turquoise module, and 122 genes within the yellow module. The 967 genes that may not be enclosed in any module were placed in the gray module, which was reserved for genes considered as not co-expressed.

FIGURE 1

Figure 1. Clustering of samples and determination of soft-thresholding power. (A) Analysis of the scale-free fit index for various soft-thresholding powers (β). The red line represents the merging threshold. (B) Analysis of the mean connectivity for various soft-thresholding powers. In all, 7 was the most fit power value.

FIGURE 2

Figure 2. Construction of coexpression modules by the WGCNA package in R. The cluster dendrogram of genes in GSE143765. Branches of the cluster dendrogram of the most connected genes gave rise to 23 gene coexpression modules.Genes that could not be clustered into one of these modules were assigned to the gray module. Every gene represents a line in the hierarchical cluster. The distance between two genes is shown as the height on the y-axis.

Correlation Between Modules and Identification of Key Modules

Interactions between the 23 co-expression modules were analyzed, and every gene is depicted within the network heatmap (Figure 3). Interestingly, the results revealed that a number of these sequence modules were valid, particularly, the light-green and inexperienced modules, with a high level of independence among the modules and relative independence of genes expressed in every module. To grasp the link between these 23 co-expression modules, we analyzed the genetic linkage of those 23 co-expression modules. Combined (Figure 4), we ascertained that these 23 modules were classified into two main clusters: one enclosed five modules (green, light yellow, light-cyan, cyan, and dark-red modules), whereas the other enclosed 18 modules (light-green, green-yellow, tan, dark-turquoise, royal-blue, pink, salmon, turquoise, yellow, dark-green, red, black, brown, grey60, purple, blue, magenta, and midnight-blue modules). Furthermore, the results were incontestable by the heatmap of eigengenes (Figure 5). We found that the proximity between modules in the two clusters was very high. The proximity between the cyan module and dark-red module was the highest in the five modules cluster, and the proximity between the tan module and green-yellow module, pink module, and royal-blue module was the highest in the cluster of 23 modules. The final correlation analysis and heat map plotting of the characteristic genes in all modules with the clinical information of SARA and non-SARA samples showed that (Figure 6) only the green and lightweight modules were significantly correlated with SARA (P < 0.05) and negatively correlated with non-SARA (P < 0.05).

FIGURE 3

Figure 3. Interaction relationship analysis of coexpressed genes. Different colors of the horizontal axis and vertical axis represent different modules. The brightness of yellow in the middle represents the degree of connectivity of different modules. There was no significant difference in interactions among different modules, indicating a high degree of independence among these modules.

FIGURE 4

Figure 4. Hierarchical clustering of module hub genes that summarize the modules yielded in the clustering analysis.

FIGURE 5

Figure 5. Heat map based on the results of C plots, with darker red being higher correlations.

FIGURE 6

Figure 6. Correlation analysis of all modules with SARA samples and Non-SARA samples, positive correlations in red and negative correlations in green, values in parentheses are p-values, correlations were considered significant when p-values were < 0.05.

Functional Enrichment Analysis of Genes in the Module

As only the green and light-yellow modules were highly correlated with SARA, GO, and KEGG enrichment analyses were performed to explore the biological functions of the genes within the green and light-yellow modules. As shown in Figure 7A, genes in the green and light-yellow modules were mainly enriched in the cytoplasm to function with the molecular functions of L-xylulose reductase (NADP+) activity and carbonyl reductase (NADPH) activity, which are involved in the negative regulation of the collagen catabolic process, Rho protein signal transduction, glucose metabolic process, xylulose metabolic process, and positive regulation of pseudopodium assembly. KEGG pathway analysis (Figure 7B) showed proximal tubule bicarbonate reclamation, pentose and glucuronate interconversions, and aldosterone-regulated sodium reabsorption. Pyruvate metabolism and endocrine and other factor-regulated calcium reabsorption. These results showed that these genes were closely related to the nutritional metabolism and absorption functions such as carbohydrate, pyruvate, and calcium ion reabsorption.

FIGURE 7

Figure 7. GO analysis and KEGG pathway enrichment results of genes in green and light yellow module. (A) GO enrichment results for the green and light yellow module, the bubble size shows the number of enriched GO terms containing genes from the green and light yellow module, the color represents the significance score of the enrichment, and the bubble nodes represent the ratio of the enriched genes to the GO terms. (B) KEGG enrichment of green and light yellow modules, bubble size shows the number of genes enriched to the biological pathway lexicon containing green and light yellow modules, color represents the color represents the significance score of enrichment, and bubble nodes represent the ratio of enriched genes to the KEGG lexicon.

PPI Network Construction and Module Analysis

The PPI network for all genes in the green and light yellow modules (Figure 8A) was constructed based on the STRING database and Cytoscape software. A high degree was calculated for the hub genes of the green and light-yellow modules using the cytoHubba plugin. As shown in Figure 8B, we crossed the genes in the green and light yellow modules with the differential genes and found that TPRG1L, BDKRB2, CLP1, KLK4, RBP4, and DCXR were the common genes. As shown in Table 1, we crossed genes in the green and light-yellow module with genes enriched in GO, and finally MMP15, MMP17, and DCXR genes were found common. Therefore, dicarbonyl and L-xylulose reductase (DCXR) is the crucial important HUB gene, because it has differential, modular, and HUB genes, whereas MMP15 (matrix metallopeptidase 15) and MMP17 (matrix metallopeptidase 17) have both modular genes and the HUB gene. DCXR, MMP15, and MMP17 are enriched in biological processes and molecular functions, such as collagen catabolic process, L-xylulose reductase (NADP+) activity, glucose metabolism process, xylulose metabolism process, and carbonyl reductase (NADPH) activity.

FIGURE 8

Figure 8. PPI network of hub genes in green and light yellow modules. (A) Green and light yellow module, with 79 nodes and 119 edges. Triangular nodes represent hub genes. (B) The blue circle is the differential gene, the red circle is the blue module and the light yellow module gene, and the common gene is the junction of the two circles.

TABLE 1

Table 1. Most significant GO terms of the turquoise module.

Discussion

Subacute ruminal acidosis (SARA) is a prevalent nutritional metabolic disorder in large-scale dairy farms (32, 33). Despite the progress made regarding the nutritional metabolism of SARA, further research is needed to investigate the relationship between SARA and the pathways related to immunosuppression and inflammatory responses. In recent years, microarray technology has become a popular technique that is often used to identify genes with altered expression during disease pathology (34). The GSE143765 dataset may be important for identifying SARA case physiology and biomarkers. Several researchers have previously analyzed this dataset, and the raw data of GSE143765 provided by Kizaki and Kim demonstrated that SARA differential genes in the liver were mainly enriched for genes involved in oxidation, suggesting that SARA is associated with high oxidative stress. Dong et al. divided sheep into a 40% low concentrate feed group and a 60% high concentrate feed group fed continuously for 9 weeks (35). The same was true for Abaker et al., who divided mid-lactation cows into 40% low and 60% high concentrate groups and found that the high concentration group developed SARA and had greater NF-kB gene expression and an inflammatory response, leading to increased oxidative stress (36). In the present study, WGCNA was used to analyze the molecular mechanisms of SARA samples compared to non-SARA samples. A total of 14,590 original genes were used to establish co-expression networks and to identify high co-expression gene sets. Moreover, 23 different co-expression modules were identified based on the genes processed by the original gene, and the functions of the modules were analyzed based on the co-expression network to obtain the modules most related to SARA (green and light yellow modules).

Over the past 20 years, research has largely established that when SARA occurs in cows, metabolic disorders, inflammatory responses, and immunosuppression primarily occur (37). Many studies have shown that compared with the non-SARA group, a large amount of VFA is produced in the rumen of dairy cows in the SARA group. When VFA is higher than the limit of VFA absorption in the rumen, the pH value in the rumen is reduced to damage the rumen barrier, affect the absorption of nutrients in the rumen, and finally accumulate in the liver, causing nutritional metabolic regulation disorders (38). The green and light-yellow module genes were mainly involved in nutritional metabolism. The enrichment of glucose and xylose in the sugar metabolism is due to the fact that when the lactating cows have SARA symptoms, the energy provided by the feed cannot be absorbed due to the damage of the rumen. Therefore, the liver mobilizes the decomposition of adipose tissue and reduces insulin sensitivity to provide energy for the body (39). The reason for the enrichment of L-xylulose reductase (NADP+) and carbonyl reductase (NADPH) activities is that NADPH is the main component in the process of obtaining energy in dairy cows (40). The collagen catabolic process, Rho protein signal transduction, and positive regulation of pseudopodium assembly are mainly enriched because LPS is produced in the rumen tissue when SARA occurs in dairy cows, and liver damage and systemic inflammation are produced by blood circulation (41). The genes in the green and light-yellow modules were mainly enriched in functions related to carbohydrate metabolism, energy metabolism, and immune regulation, which indicates that the genes in the green and light-yellow modules are involved in the nutritional metabolism and inflammatory response process in dairy cows suffering from SARA diseases. Through the analysis of the PPI network, differential gene screening, and functional enrichment of genes in the green and light-yellow modules, the genes with intersection of the two and above were identified as DCXR, and MMP15 and MMP17 were the key genes. Many studies have shown that DCXR plays a crucial role in fat and glucose metabolism. Palombo et al. also used DCXR as a candidate marker gene for fat metabolism in genome-wide association studies of milk fatty acid composition in Italian Holstein cows (42). MMP15 and MMP17, as HUB genes, are involved in embryonic development, reproduction, and tissue damage recovery during normal physiological processes, however, they promote inflammation, such as arthritis and fibrosis, when the disease occurs. (43), Therefore, although the results of MMP17 and MMP15 are low, it is not difficult to find that both of them increase gene expression after SARA disease, hence, they should be involved in the inflammatory response and damage repair of the rumen and liver tissues (44).

Pathway enrichment analysis showed that genes in green and light-yellow modules were enriched in proximal tubule bicarbonate reclamation, pentose and glucuronate interconversions, aldosterone-regulated sodium reabsorption, pyruvate metabolism, and endocrine and other factor-regulated calcium reabsorption pathways. The inflammatory response is typically mediated by the release of large amounts of inflammatory mediators, activation of inflammatory cells such as monocytes and macrophages, and release of pro-inflammatory mediators, such as TNF-α and interleukin-1B (IL-1B) (45). In recent years, it has been found that the SARA-induced inflammatory response mainly occurs in rumen tissue, mainly because when suffering from SARA, the pH value in the rumen decreases, leading to bacterial fragmentation, releasing a large amount of LPS to damage the rumen tissue barrier, which eventually leads to LPS translocation to the peripheral blood circulation (46). Because the liver is an important organ of blood circulation, translocated LPS and LPS binding protein (LBP) eventually enter the liver, leading to liver disease (47). However, in recent years, in-depth studies on the inflammatory response of SARA have found that when dairy cows suffer from SARA, stress is induced and calcium, non-esterified fatty acids, and glucose concentrations are reduced, and glucocorticoid receptor protein in the liver is significantly downregulated, which is consistent with the enrichment results of this experiment (48). Although LPS translocation eventually leads to liver inflammation, liver inflammation after a series of reactions has rarely been explored.

This study had certain limitations. First, this study focused on two aspects of data mining and data analysis, both of which are genetic methodologies, and the results were not validated experimentally. Second, the samples were obtained from peripartal Holstein cow liver samples; hence, only the liver gene expression analysis during SARA was performed and association analysis with blood samples was not performed. Therefore, the results of this study should be understood within the context of its limitations.

Conclusions

Our study adopted a systems biology-based WGCNA methodology and uncovered various helpful molecular targets for future investigations of the mechanisms and choice of SARA biomarkers. Some essential biological processes and pathways, collagen catabolic process, L-xylulose reductase (NADP+) activity, glucose metabolism, xylulose metabolism, NADPH activity, pentose and glucuronate interconversions, and other pathways, as well as hub genes that play a role in these processes, may help to elucidate the process of SARA production and development. Moreover, potential molecular mechanisms include Rho protein signal transduction, positive regulation of pseudopodium assembly and proximal tubule bicarbonate reclamation, aldosterone-regulated sodium reabsorption, pyruvate metabolism, and endocrine and other factor-regulated calcium reabsorption signaling pathways. Our findings may facilitate the development of new experiments to better understand the pathophysiology and progression of SARA. Further studies should be conducted to validate the value of the resulting genes in the context of SARA.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

Author Contributions

QW, BG, and CX conceived the study. QW, BG, and XY participated in method development and validation. BG and YC carried out the data analysis. QW and BG prepared the original draft. JL, XD, and XW reviewed and modified the manuscript. All the authors have read and approved the final version of the manuscript.

Funding

This work was supported by Major Science and Technology Project of Heilongjiang Province (Grant No. 2021ZX12B03), the Provincial Institute Cooperation Project of the Heilongjiang Science and Technology Plan (Grant Nos. YS20B04 and YS19B01), the Development Project of Local Universities in Heilongjiang Bayi Agricultural University (ZRCQC201803 and ZRCLG201904), the Postdoctoral Scientific Research Developmental Fund of Heilongjiang Province (LBH-Q20161), and the National Natural Science Foundation of China, Beijing, China (Grant Nos. 32125038 and 32072931).

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.

Abbreviations

SARA, subacute ruminal acidosis; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; BP, biological process; MF, molecular function; CC, cellular component. PPI, protein-protein interaction.

References

1. Kleen JL, Hooijer GA, Rehage J, Noordhuizen JP. Subacute ruminal acidosis (SARA): a review. J Vet Med A Physiol Pathol Clin Med. (2003) 50:406–14. doi: 10.1046/j.1439-0442.2003.00569.x

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Gozho GN, Krause DO, Plaizier JC. Ruminal lipopolysaccharide concentration and inflammatory response during grain-induced subacute ruminal acidosis in dairy cows. J Dairy Sci. (2007) 90:856–66. doi: 10.3168/jds.S0022-0302(07)71569-2

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Enemark JM. The monitoring prevention and treatment of sub-acute ruminal acidosis (SARA): a review. Vet J. (2008) 176:32–43. doi: 10.1016/j.tvjl.2007.12.021

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Minami NS, Sousa RS, Oliveira FLC, Dias MRB, Cassiano DA, Mori CSET AL. Subacute ruminal acidosis in zebu cattle: clinical and behavioral aspects. Animals (Basel). (2020) 11:21. doi: 10.3390/ani11010021

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Kebreab E, Dijkstra J, Bannink A, France J. Recent advances in modeling nutrient utilization in ruminants. J Anim Sci. (2009) 87:E111–22. doi: 10.2527/jas.2008-1313

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Plaizier JC, Krause DO, Gozho GN, McBride BW. Subacute ruminal acidosis in dairy cows: the physiological causes, incidence and consequences. Vet J. (2008) 176:21–31. doi: 10.1016/j.tvjl.2007.12.016

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Morgante M, Stelletta C, Berzaghi P, Gianesella M, Andrighetto I. Subacute rumen acidosis in lactating cows: an investigation in intensive Italian dairy herds. J Anim Physiol Anim Nutr (Berl). (2007) 91:226–34. doi: 10.1111/j.1439-0396.2007.00696.x

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Yan Q, Wang K, Han X, Tan Z. The regulatory mechanism of feeding a diet high in rice grain on the growth and microrna expression profiles of the spleen, taking goats as an artiodactyl model. Biology (Basel). (2021) 10:832. doi: 10.3390/biology10090832

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Gholizadeh M, Fayazi J, Asgari Y, Zali H, Kaderali L. Reconstruction and analysis of cattle metabolic networks in normal and acidosis rumen tissue. Animals (Basel). (2020) 10:469. doi: 10.3390/ani10030469

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Zhao C, Bobe G, Wang Y, Zhang X, Zhao Z, Zhang S, et al. Potential role of SLC5A8 expression in the etiology of subacute ruminal acidosis. Front Vet Sci. (2020) 7:394. Published 2020 Jul 30. doi: 10.3389/fvets.2020.00394

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Sato S. Subacute ruminal acidosis (SARA) challenge ruminal condition and cellular immunity in cattle. Jpn J Vet Res. (2015) 63:S25–36.

PubMed Abstract | Google Scholar

12. Maeda M, Kawasumi K, Sato S, Arai T. Evaluation of blood adiponectin levels as an index for subacute ruminal acidosis in cows: a preliminary study. Vet Res Commun. (2019) 43:215–24. doi: 10.1007/s11259-019-09760-0

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Gozho GN, Plaizier JC, Krause DO, Kennedy AD, Wittenberg KM. Subacute ruminal acidosis induces ruminal lipopolysaccharide endotoxin release and triggers an inflammatory response. J Dairy Sci. (2005) 88:1399–403. doi: 10.3168/jds.S0022-0302(05)72807-1

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Stefanska B, Człapa W, Pruszynska-Oszmałek E, Szczepankiewicz D, Fievez V, Komisarek J, et al. Nowak W. Subacute ruminal acidosis affects fermentation and endotoxin concentration in the rumen and relative expression of the CD14/TLR4/MD2 genes involved in lipopolysaccharide systemic immune response in dairy cows. J Dairy Sci. (2018) 101:1297–310. doi: 10.3168/jds.2017-12896

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Humer E, Aditya S, Zebeli Q. Innate immunity and metabolomic responses in dairy cows challenged intramammarily with lipopolysaccharide after subacute ruminal acidosis. Animal. (2018) 12:2551–60. doi: 10.1017/S1751731118000411

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Mahmoudi B, Fayazi J, Roshanfekr H, Sari M, Bakhtiarizadeh MR. Genome-wide identification and characterization of novel long non-coding RNA in Ruminal tissue affected with sub-acute Ruminal acidosis from Holstein cattle. Vet Res Commun. (2020) 44:19–27. doi: 10.1007/s11259-020-09769-w

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Do DN, Dudemaine PL, Fomenky BE, Ibeagha-Awemu EM. Integration of miRNA weighted gene co-expression network and miRNA-mRNA co-expression analyses reveals potential regulatory functions of miRNAs in calf rumen development. Genomics. (2019) 111:849–59. doi: 10.1016/j.ygeno.2018.05.009

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Sheybani N, Bakhtiarizadeh MR, Salehi A. An integrated analysis of mRNAs, lncRNAs, and miRNAs based on weighted gene co-expression network analysis involved in bovine endometritis. Sci Rep. (2021) 11:18050. doi: 10.1038/s41598-021-97319-y

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Tsuchiya Y, Ozai R, Sugino T, Kawashima K, Kushibiki S, Kim YH, et al. Changes in peripheral blood oxidative stress markers and hepatic gene expression related to oxidative stress in Holstein cows with and without subacute ruminal acidosis during the periparturient period. J Vet Med Sci. (2020) 82:1529–36. doi: 10.1292/jvms.20-0426

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Tusher VG, Tibshirani R, Chu G. Significance analysis of microarrays applied to the ionizing radiation response [published correction appears in Proc Natl Acad Sci USA 2001 Aug 28; 98:10515]. Proc Natl Acad Sci USA. (2001) 98:5116–21. doi: 10.1073/pnas.091062498

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Bakhtiarizadeh MR, Mirzaei S, Norouzi M, Sheybani N, Vafaei Sadi MS. Identification of gene modules and hub genes involved in mastitis development using a systems biology approach. Front Genet. (2020) 11:722. doi: 10.3389/fgene.2020.00722

PubMed Abstract | CrossRef Full Text | Google Scholar

23. de Oliveira PSN, Coutinho LL, Cesar ASM, Diniz WJDS, de Souza MM, Andrade BG, et al. Co-expression networks reveal potential regulatory roles of miRNAs in fatty acid composition of nelore cattle. Front Genet. (2019) 10:651. doi: 10.3389/fgene.2019.00651

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Zhao W, Langfelder P, Fuller T, Dong J, Li A, Hovarth S. Weighted gene coexpression network analysis: state of the art. J Biopharm Stat. (2010) 20:281–300. doi: 10.1080/10543400903572753

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. (2008) 9:559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Wang J, Chai Z, Deng L, Wang J, Wang H, Tang Y, et al. Detection and integrated analysis of lncRNA and mRNA relevant to plateau adaptation of Yak. Reprod Domest Anim. (2020) 55:1461–9. doi: 10.1111/rda.13767

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Dennis G Jr, Sherman BT, Hosack DA, Yang J, Gao W, Lane HC, et al. DAVID: database for annotation, visualization, and integrated discovery. Genome Biol. (2003) 4:P3. doi: 10.1186/gb-2003-4-5-p3

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Huang DW, Sherman BT, Tan Q, Kir J, Liu D, Bryant D, et al. DAVID bioinformatics resources: expanded annotation database and novel algorithms to better extract biology from large gene lists. Nucleic Acids Res. (2007) 35:W169–75. doi: 10.1093/nar/gkm415

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Chen R, Ge T, Jiang W, Huo J, Chang Q, Geng J, et al. Identification of biomarkers correlated with hypertrophic cardiomyopathy with co-expression analysis. J Cell Physiol. (2019) 234:21999–2008. doi: 10.1002/jcp.28762

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Szklarczyk D, Morris JH, Cook H, Kuhn M, Wyder S, Simonovic M, et al. The STRING database in 2017: quality-controlled protein-protein association networks made broadly accessible. Nucleic Acids Res. (2017) 45:D362–8. doi: 10.1093/nar/gkw937

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. (2003) 13:2498–504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Tsuchiya Y, Chiba E, Sugino T, Kawashima K, Kushibiki S, Kizaki K, et al. Liver transcriptome response to periparturient hormonal and metabolic changes depends on the postpartum occurrence of subacute ruminal acidosis in Holstein cows. Physiol Genomics. (2021) 53:285–94. doi: 10.1152/physiolgenomics.00048.2021

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Zhao C, Wang Y, Peng Z, Sun X, Sun G, Yuan X, et al. Subacute ruminal acidosis suppressed the expression of MCT1 in rumen of cows. J Cell Physiol. (2019) 234:11734–45. doi: 10.1002/jcp.27829

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Kim YH, Toji N, Kizaki K, Kushibiki S, Ichijo T, Sato S. Effects of dietary forage and calf starter on ruminal pH and transcriptomic adaptation of the rumen epithelium in Holstein calves during the weaning transition. Physiol Genomics. (2016) 48:803–9. doi: 10.1152/physiolgenomics.00086.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Dong H, Wang S, Jia Y, Ni Y, Zhang Y, Zhuang S, et al. Long-term effects of subacute ruminal acidosis (SARA) on milk quality and hepatic gene expression in lactating goats fed a high-concentrate diet. PLoS ONE. (2013) 8:e82850. doi: 10.1371/journal.pone.0082850

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Abaker JA, Xu TL, Jin D, Chang GJ, Zhang K, Shen XZ. Lipopolysaccharide derived from the digestive tract provokes oxidative stress in the liver of dairy cows fed a high-grain diet. J Dairy Sci. (2017) 100:666–78. doi: 10.3168/jds.2016-10871

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Esposito G, Irons PC, Webb EC, Chapwanya A. Interactions between negative energy balance, metabolic diseases, uterine health and immune response in transition dairy cows. Anim Reprod Sci. (2014) 144:60–71. doi: 10.1016/j.anireprosci.2013.11.007

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Aditya S, Humer E, Pourazad P, Khiaosa-Ard R, Zebeli Q. Metabolic and stress responses in dairy cows fed a concentrate-rich diet and submitted to intramammary lipopolysaccharide challenge. Animal. (2018) 12:741–9. doi: 10.1017/S1751731117002191

PubMed Abstract | CrossRef Full Text | Google Scholar

39. McFadden JW. Review: Lipid biology in the periparturient dairy cow: contemporary perspectives. Animal. (2020) 14:s165–75. doi: 10.1017/S1751731119003185

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Chen J, Yang Z, Dong G. Niacin nutrition and rumen-protected niacin supplementation in dairy cows: an updated review. Br J Nutr. (2019) 122:1103–12. doi: 10.1017/S0007114519002216

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Hafez SMNA, Elbassuoni E. Dysfunction of aged liver of male albino rats and the effect of intermitted fasting) Biochemical, histological, and immunohistochemical study. Int Immunopharmacol. (2022) 103:108465. doi: 10.1016/j.intimp.2021.108465

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Palombo V, Milanesi M, Sgorlon S, Capomaccio S, Mele M, Nicolazzi E, et al. Genome-wide association study of milk fatty acid composition in Italian Simmental and Italian Holstein cows using single nucleotide polymorphism arrays. J Dairy Sci. (2018) 101:11004–19. doi: 10.3168/jds.2018-14413

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Zheng J, Zeng P, Zhang H, Zhou Y, Liao J, Zhu W, et al. Long noncoding RNA ZFAS1 silencing alleviates rheumatoid arthritis via blocking miR-296-5p-mediated down-regulation of MMP-15. Int Immunopharmacol. (2021) 90:107061. doi: 10.1016/j.intimp.2020.107061

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Yip C, Foidart P, Noël A, Sounni NE. MT4-MMP: The GPI-anchored membrane-type matrix metalloprotease with multiple functions in diseases. Int J Mol Sci. (2019) 20:354. doi: 10.3390/ijms20020354

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Khafipour E, Krause DO, Plaizier JC. A grain-based subacute ruminal acidosis challenge causes translocation of lipopolysaccharide and triggers inflammation. J Dairy Sci. (2009) 92:1060–70. doi: 10.3168/jds.2008-1389

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Chang G, Zhuang S, Seyfert HM, Zhang K, Xu T, Jin D, et al. Hepatic TLR4 signaling is activated by LPS from digestive tract during SARA, and epigenetic mechanisms contribute to enforced TLR4 expression. Oncotarget. (2015) 6:38578–90. doi: 10.18632/oncotarget.6161

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Pan X, Nan X, Yang L, Jiang L, Xiong B. Thiamine status, metabolism and application in dairy cows: a review. Br J Nutr. (2018) 120:491–9. doi: 10.1017/S0007114518001666

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Chen Q, Wu C, Yao Z, Cai L, Ni Y, Mao S, Zhao R. Whole transcriptome analysis of RNA expression profiles reveals the potential regulating action of long noncoding RNA in lactating cows fed a high concentrate diet. Anim Nutr. (2021) 7:1315–28. doi: 10.1016/j.aninu.2021.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar