Abstract
Nutrition impacts the epigenetic signature, including DNA methylation. The aim of this study was to identify genomic regions differentially methylated in the rumen of Italian Mediterranean dairy buffaloes fed green forage [Total Mixed Ration (TMR) + ryegrass green feed (30% of diet)] compared to those receiving a standard TMR diet, through Reduced Representation Bisulfite Sequencing. We found 6571 differentially methylated genomic regions (DMRs), 51.73% hypomethylated and 48.27% hypermethylated. DMRs were uniformly dispersed in genes and intergenic regions and along chromosomes. Genes-associated DMRs were mainly hypomethylated, while intergenic DMRs were mostly hypermethylated. We highlighted 4648 genes associated with DMRs (differentially methylated genes, DMGs), mostly protein-coding genes. Gene Ontology study performed with hypermethylated or hypomethylated DMGs highlighted categories related to response to oxidative stress and inflammation, as well as rumen functionality. The integration of our results with differential expression data identified genes whose expression varies as a function of DNA methylation. This subset of genes included those involved in immune system functioning, inflammation, fatty acid metabolism, and stress response. Our findings highlighted the impact of green forage on rumen DNA methylation, which potentially influences molecular mechanisms relevant to rumen functionality and, then, animal welfare and production, through the modulation of gene expression.
Similar content being viewed by others
Introduction
Buffalo is a globally important livestock specie since it provides economic value from milk, meat, leather, and draft power in both developed and developing countries1. Buffaloes are known worldwide for their higher resistance to hot and humid climate2, parasites3, and to have better adaptability to low-quality and less-digestible roughage4.
Recent research on diverse feeding practices has shifted consumer preferences towards quality more than quantity. There is now a growing consumer interest in products sourced from grass-fed animals, raising questions about the quality distinctions between grass-fed and grain-fed animals. Over the past few decades, it has become evident that different feeding practices can alter the nutritional composition of animal-derived products5,6. In particular, it has been demonstrated that a rich and balanced diet plays a pivotal role in milk and dairy production, although there are some features concerning rumen physiology that should be considered when preparing the feed ration7,8.
The metabolic processes of rumen can influence the composition of milk nutrients. The byproducts of rumen fermentation, including short-chain fatty acids (SCFAs) generated in the forestomachs, are primarily absorbed through the epithelium of both the rumen and omasum9. Alongside other beneficial molecules, they are then transferred into the bloodstream, metabolized by the liver, and ultimately transported to the mammary gland10. Notably, SCFAs fulfill roughly 80% of the energy requirements for these animals, hence playing crucial roles in maintaining energy balance in ruminants11. Therefore, the effect of diet on the rumen physiology should be considered when preparing the feed ration7,8. In these regard, green forage addition in the diet of ruminants increases the content of nutraceutical molecules12 such as vaccenic and rumenic acids13, and betaines and carnitines6.
It is known that nutrition might impact on the epigenetic signature, such as DNA methylation. In vertebrate DNA methylation is an important epigenetic feature involved in the regulation of many biological processes, including embryonic development, chromosome stability, chromatin structure, and transcription14,15,16,17. Nutrients, such as choline, folate, and betaine have a significant effect on DNA methylation, which act as methyl-group donors18,19. Furthermore, it has been demonstrated that nutrient quality and quantity modulate epigenetic landscape in livestock20. However, to date, few studies have been performed on DNA methylation in ruminants21, and only one study has correlated the effect of different diets with changes of DNA methylation and gene expression in cattle22.
Recently, we demonstrated that green forage, in Italian Mediterranean dairy buffaloes, modulates the expression of genes involved in biological processes relevant for rumen functionality and physiology23. Given the crucial role of DNA methylation in gene regulation, and its susceptibility to diet, it is important to understand if green forage might impact on ruminal transcriptional program through the modulation of DNA methylation profiles.
With this study, we aimed to identify genomic regions differentially methylated in the ruminal wall of Italian Mediterranean dairy buffaloes fed green forage [Total Mixed Ration (TMR) + ryegrass green feed (30% of diet)], in comparison with animals that received a standard TMR diet, through Reduced Representation Bisulfite Sequencing (RRBS). Moreover, we investigated the correlations between changes in DNA methylation and gene expression, by integrating DNA methylation data with transcriptomic profiles previously obtained from the same rumen samples23. Our findings shed light on the impact of green forage on one of the most studied epigenetic modifications at the ruminal level, which might influence molecular mechanisms relevant for animal welfare.
Results
Reduced representation bisulfite sequencing and Methylomic analysis in the ruminal wall of dairy buffaloes
We obtained methylomic profiles from the ruminal wall of Italian Mediterranean dairy buffaloes fed green forage (8 animals, Treated group) or standard TMR diet (6 animals, Control group), through RRBS. RRBS yielded an average of 21,524,186 high quality 100 bp reads per sample. Approximately, the 70% of the reads were mapped to the reference genome sequence of Bubalus bubalis (version UOA_WB_1 from NCBI) and, more than 50% of them were uniquely mapped. After deduplication processing for the correction of amplification bias using Unique Molecular Identifiers (UMIs), the uniquely mapped reads percentage was, on average, 40% (Supplementary Table S1). Overall, the average of CpG sites identified per sample was 6,195,365 (Supplementary Table S2), of which, 3,456,649 were shared between all the samples and represent the 2.52% of the CpG content in the Bubalus bubalis genome (data not shown). All samples showed a similar distribution of CpG coverage that spans between 1 and 25 reads, with a large number of CpGs covered by less than 10 reads (Supplementary Figure S1). The average number of methylated CpG sites per sample, obtained assuming that the CpG site with the sequence corresponding to the methylated state was covered by at least one read, was 3,581,857 (Supplementary Table S2). Among the CpG sites shared between all the samples analyzed (3,456,649, see above), 1,815,690 CpGs (52,5%) were methylated in all the samples of the Control group, while 1,614,748 CpGs (46,7%) were methylated in all the samples of the Treated group (data not shown).
Subsequently, the genomic methylation has been compared between the Treated and Control groups. Principal component analysis revealed a significant clustering (p = 0.0156) of the Treated samples, which is less pronounced in the Control group (Fig. 1a). The differential methylation analysis between Treated and Control groups allowed to detect 6571 statistically significant differentially methylated genomic regions (DMRs), containing at least two CpG sites and showing a coverage ≥ 5 reads. Among them, 3399 DMRs showed hypomethylation (51.73%) and 3172 DMRs (48.27%) were hypermethylated in the Treated versus the Control group (Fig. 1b and c, left panel; Supplementary Table S3). The CpG sites included in the DMRs are more than 0.045% of the total number of CpG sites in the Bubalus bubalis genome; 0.03% of these CpGs are included in the hypo-DMRs, while the 0.015% are included in hyper-DMRs (Fig. 1c, right panel).
Effect of green forage on genomic DNA methylation in the rumen of dairy buffaloes. (a) Principal component analysis that illustrates the methylation status of the CpG sites in all the samples. The blue and yellow dots represent samples of the Treated and Control groups, respectively. (b) MA plot showing differential methylation of all DMRs, reported as mean differential methylation versus the mean methylation between the two experimental groups. Red and green dots indicate the DMRs that are significantly hypermethylated and hypomethylated, respectively, in the Treated versus Control groups. Black dots are the DMRs whose differential methylation did not reach the statistical significance between the two experimental groups. (c) Left: Number of total DMRs, hypo-DMRs and hyper-DMRs identified in Treated versus Control animals. Right: Percentage of CpG sites included in total DMRs, hypo-DMRs and hyper-DMRs with respect to the total number of CpG sites in the Bubalus bubalis genome.
Green forage impacts at the genome-wide level on DNA methylation in the ruminal wall of dairy buffaloes
Next, we investigated the chromosomal distribution of the DMRs, to understand whether a preferential distribution of the DMRs on specific chromosomes and/or chromosomal regions occurs. This analysis revealed a uniform distribution of DMRs along all chromosomes, without any preference for chromosomal features (Fig. 2a). Moreover, the number of DMRs per chromosome positively correlates with both the number of genes in each chromosome (Fig. 2b) and the chromosome size (Fig. 2c), as expected.
Genomic distribution of DMRs in the rumen of dairy buffaloes fed green forage or standard diet. (a) Chromosomal distribution of the DMRs identified in the Treated as compared with the Control group. Chromosomal length is displayed on the x axis (Mb). Chromosomes are sorted by size. Red bars indicate hypermethylated DMRs, while green bars show hypomethylated DMRs. The presence of multiple DMRs in the same genomic region is indicated by longer bars. The percentage of DMRs in each chromosome is reported on the right of the figure. (b) and (c) Scatterplots that show the correlation of the number of DMRs with the number of genes on the different chromosomes (b) or with the chromosome size (c). The Pearson’s correlation coefficients and the statistical significance are displayed. (d) Left: Percentage of DMRs localized in intergenic and genic (exons, introns and promoters) regions. Right: Number of hypomethylated and hypermethylated DMRs included in intergenic regions, exons, introns and promoters. For promoters, DMRs located up to 2 kb upstream of the transcriptional start site is considered. (e) Pie chart showing the distribution of the DMGs according to their biotype.
We, further, analyzed the genomic distribution of the DMRs with respect to genes and intergenic regions. We found that DMRs are distributed similarly between genic (55%) and intergenic (45%) regions. Among DMRs located in genic regions, 22% were located at the promoter regions (up to 2 kb upstream of the transcriptional start site), 18% at the introns, and 15% at the exons (Fig. 2d, left). Moreover, intergenic DMRs are predominantly hypermethylated, as opposed to what has been observed for genic DMRs (Fig. 2d, right).
Genes associated with DMRs
Subsequently, we identified genes located in the proximity (distance < 10 kb) of, at least, one DMR, which we named Differentially Methylated Genes (DMGs) (Supplementary Table S4). Overall, we identified 4648 DMGs, 3830 (82.4%) of which are protein-coding genes, 509 (10.95%) are long non-coding genes, 181 (3.89%) are pseudogenes, and 128 (2.75%) belong to other biotypes (Fig. 2e). Several DMGs are associated with more than one DMR, and the same DMR can be associated with more than one gene (Supplementary Table S4).
Using the list of identified DMGs (Supplementary Table S4), we focused the attention on genes associated with significant hyper-DMRs (Table 1) and hypo-DMRs (Table 2) that showed the highest differential methylation (top hyper- and hypo-DMRs). The list of genes associated with the top hyper-DMRs (Table 1), includes Mitochondrial Calcium Uptake 2 (MICU2), Potassium Channel Tetramerization Domain Containing 1 (KCTD1), Nucleolar Protein 4 Like (NOL4L), Protein Kinase D1 (PRKD1), RuvB Like AAA ATPase 1 (RUVBL1), SEC61 Translocon Subunit Alpha 1 (SEC61A1), TBC1 Domain Family Member 14 (TBC1D14), Sprouty RTK Signaling Antagonist 2 (SPRY2), and Translocase Of Outer Mitochondrial Membrane 20 (TOMM20). MICU2 encodes a key regulator of the mitochondrial calcium uniporter; KCTD1 encodes a factor containing a POxvirus and Zinc finger (POZ) protein-protein interaction domain, which is involved in the regulation of the Wnt signaling pathway; PRKD1 encodes an enzyme involved in several biological processes, such as cell migration and differentiation, cell survival and many signaling pathways; RUVBL1 encodes an ATPase that interacts with complexes with a role in the epigenetic regulation of pro-inflammatory factors; SEC61A1 encodes a factor involved in the calcium homeostasis; TBC1D14 encodes a factor that negatively regulates autophagy; SPRY2 encodes a factor belonging to the sprouty family that modulates the tyrosine kinase receptor signaling; TOMM20 encodes a subunit of a translocase important for the recognition and translocation of mitochondrial proteins from the cytosol into the mitochondria.
Among genes associated with the top hypo-DMRs (Table 2), we identified Transmembrane Protein 245 (TMEM245), Catenin Alpha Like 1 (CTNNAL1), Cryptochrome Circadian Regulator 2 (CRY2), Mitogen-Activated Protein Kinase 8 Interacting Protein (MAPK8IP1), Fibronectin Leucine Rich Transmembrane Protein 2 (FLRT2), Cyclin Y (CCNY), Cholinergic Receptor Muscarinic 3 (CHRM3), Proprotein Convertase Subtilisin/Kexin Type 5 (PCSK5), Ring Finger Protein 217 (RNF217), Sestrin 3 (SESN3), Calcium Binding Protein 1 (CABP1), and HECT Domain E3 Ubiquitin Protein Ligase 3 (HECTD3). TMEM245 encodes a factor hypothesized to be an integral component of membrane; CTNNAL1 encodes a factor that shows similarities to vinculin and α-catenin and is predicted to act as a cytoskeletal linker protein; CRY2 encodes a key component of the circadian clock complex; MAPK8IP1 encodes a negative regulator of MAPK8; FLRT2 encodes a cell adhesion molecule implicated in the regulation of embryonic vascular and neural development; CCNY encodes a factor that regulates cyclin-dependent kinases and controls cell cycle through the regulation of wnt signaling pathway; CHRM3 encodes a muscarinic cholinergic receptor controlling the smooth muscle contraction; PCSK5 encodes a convertase involved in the processing of precursors of other proteins into their biologically active products; RNF217 encodes a member of the RING1-IBR-RING24 (RBR) ubiquitin protein ligase family, proposed to have a role in apoptosis signaling; SESN3 encodes a stress-induced factor with an anti-oxidant function, involved also in the glucose homeostasis and lipid storage; CABP1 encodes a calcium binding protein similar to calmodulin, implicated in the regulation of voltage-gated calcium ion channels; HECTD3 encodes an enzyme involved in the ubiquitination of target proteins, thus controlling their degradation.
Afterwards, we investigated DMGs associated with more than three DMRs (Supplementary Table S5). This inspection highlighted, among others, Retinoic Acid Induced 1 (RAI1), Calcineurin Binding Protein 1 (CABIN1), Ubiquitin Conjugating Enzyme E2 Q2 (UBE2Q2), Ras Homolog Family Member F, Filopodia Associated (RHOF), Myotubularin Related Protein 4 (MTMR4), CXXC Finger Protein 5 (CXXC5), Bone Morphogenetic Protein 6 (BMP6), and Sulfatase 2 (SULF2) genes. RAI1, associated with one hypomethylated and two hypermethylated DMRs, encodes a transcription factor implicated in the regulation of the circadian clock. CABIN1, associated with three hypomethylated DMRs, encodes a factor involved in the T-cell receptor-mediated signal transduction pathway and in the negative regulation of P53. The enzyme encoded by UBE2Q2, which is associated with one hypomethylated and two hypermethylated DMRs, plays a role in the protein turnover and stability. The RAS homolog factor encoded by RHOF gene, associated with four hypermethylated DMRs, is a GTP-binding protein, recently associated with a pro-inflammatory activity. MTMR4, linked to three hypermethylated DMRs, encodes a phosphatase that negatively regulate innate immune responses. The zinc finger protein encoded by the CXXC5 gene, associated with two hypomethylated and three hypermethylated DMRs, is an epigenetic factor with a function in several signal transduction pathways. BMP6, associated with one hypermethylated and three hypomethylated DMRs, encodes a secreted ligand of the TGF-beta involved in several processes including iron homeostasis. The arylsulphatase encoded by the SULF2 gene, associated with one hypermethylated and three hypomethylated DMRs, mediates many signal transduction pathways by removing the 6-O-sulfate groups from heparan sulfate, thus modulating its activity. Among DMGs associated with two DMRs, we found Pyruvate Carboxylase (PC) (one hypermethylated and one hypomethylated DMRs), encoding a mitochondrial enzyme that catalyzes the initial step of gluconeogenesis.
Gene ontology analysis of genes associated with DMRs
In order to identify Gene Ontology (GO) terms enriched in genes associated with at least one DMR (Supplementary Table S4), we performed GO Enrichment Analysis (GOEA) using DMGs previously identified (Supplementary Table S4). Genes associated with hypermethylated or hypomethylated DMRs (hyper-DMGs and hypo-DMGs) were functionally associated with 170 or 128 significantly enriched GO terms (FDR-adjusted p-value < 0.05), respectively (Supplementary Tables S6 and S7).
The 17.65% (30) of GO terms enriched in hyper-DMGs were categorized as a molecular function, 10% (17), as a cellular component, and 72.36% (123) as a biological process. The 21.88% (28) of GO terms enriched in hypo-DMGs were classified as a molecular function, 14.84% (19) as a cellular component, and 63.28% (81) as a biological process. Selected GO terms enriched in hyper-DMGs and hypo-DMGs are illustrated in Fig. 3a and b. Among GO terms enriched in hyper-DMGs, we found many terms related to energy metabolism, mitochondrial respiratory function, oxidative phosphorylation, circadian rhythm, muscular function and development, fatty acid metabolism, response to hydrogen peroxide, and inflammation (response to chemokine) (Fig. 3a). Moreover, GO terms enriched in hypo-DMGs are linked to energy metabolism, mitochondrial function, respiration, and ion transport, most of which, are suggestively, related to ATP function. Furthermore, we found enrichment in categories linked to L-methionine biosynthetic process, cellular response to hormone stimuli, and to response to fungi (Fig. 3b).
Gene ontology enrichment analysis (GOEA) of genes associated with DMRs. Selected GO-terms, enriched in genes associated with hyper-DMRs (a) and hypo-DMRs (b) in the ruminal wall of dairy buffaloes fed green forage (Treated) in comparison with animals fed TMR standard diet (Control), are shown. GO categories were classified as molecular function (yellow bars), cellular component (burgundy bars) and biological process (blue bars). Bar graphs indicate the statistical significance of the enrichment, as -log10 (adj p-value). Vertical yellow bars indicate the cut-off level for significance (p < 0.05, adjusted by Benjamini-Hochberg correction). For graphical reasons, for GO terms that display an FDR-adj p-value = 0, an arbitrary value of -log10 (adj p-value) higher than all the other GO terms showed in the same graph has been used (see also Supplementary Table S6 and S7).
Relationship between DNA methylation and gene expression
Overall, our study revealed a remarkable impact of the green forage-based diet on genomic DNA methylation. In order to investigate the potential effect of the observed DNA methylation changes on the transcriptional program in dairy buffaloes fed green forage, we compared data of genomic methylation with the previously obtained data of gene expression23 in the ruminal wall of Treated buffaloes. For studies of correlation between DNA methylation and gene expression, we used only the data obtained with the ruminal samples of the same animals used for the two experiments (five samples from buffaloes fed TMR standard diet and six samples from animals fed green forage). Overall, we did not find a significant correlation between the global gene expression and the DNA methylation in the DMRs neither within the promoter regions (up to 2 kb upstream the transcriptional start site) nor within the gene bodies in rumen of the Treated buffaloes (Fig. 4a e b). However, the comparison of the previously highlighted differential expression23 with changes of DNA methylation in DMGs within their promoter or their gene body allowed to identify a discrete number of genes for which changes in the expression upon green forage administration are associated with changes in DNA methylation (Fig. 4c and d; Tables 3 and 4).
Integration of transcriptomic and methylomic results. (a) and (b) Scatterplots that display the mean gene expression, reported as log2 of normalized counts (FPKM), plotted against the mean DNA methylation, displayed as percentage, in the differentially methylated regions (DMRs) in the promoter (up to 2 kb upstream the transcriptional start site) or in the gene body (b), in rumen of the Treated animals. In each plot, the Pearson’s correlation coefficients and the statistical significance are displayed. (c) and (d) Scatterplots that show the impact of green forage on the transcriptional program and the DNA methylation in the Treated animals in comparison with the Control group. Changes in the gene expression, displayed as log2 Fold Change, are plotted against changes in the DNA methylation, expressed as percentage, in the promoters (c) or the gene bodies (d) of the DMGs. Yellow dots: hypermethylated and up-regulated genes; blue dots: hypermethylated and down-regulated genes; green dots: hypomethylated and up-regulated genes; red dots: hypomethylated and down-regulated genes.
Among DEGs for which a significant differential expression is associated with changed DNA methylation, we found Growth Factor Receptor Bound Protein 10 (GRB10), GREB1 Like Retinoic Acid Receptor Coactivator (GREB1L), SIM BHLH Transcription Factor 2 (SIM2), and Tryptophanyl tRNA Synthetase (WARS), which display differential gene expression and changes in DNA methylation within their promoter (Table 3). Moreover, BMP6, CDC42 Effector Protein 5 (CDC42EP5), Dual Oxidase 1 (DUOX1), Potassium Two Pore Domain Channel Subfamily K Member 10 (KCNK10), Kinesin Family Member 21 A (KIF21A), LOC102409914 (or Zinc Finger Protein 786, ZNF786), Peptidase M20 Domain Containing 2 (PM20D2), Sterile Alpha Motif Domain Containing 11 (SAMD11), Teneurin Transmembrane Protein 3 (TENM3), Transcription Factor AP-2 Epsilon (TFAP2E), Tenascin XB (TNXB), Tripartite Motif Containing 7 (TRIM7), and Zinc Finger FYVE-Type Containing 16 (ZFYVE16) show differential gene expression and changes in DNA methylation within their gene body (Table 4). In addition, LOC102400339 (or Mitotic Arrest Deficient 2 Like 1, MAD2L1) and Uridine Phosphorylase 1 (UPP1) are genes differentially expressed and differentially methylated in DMRs located in both promoter and gene body (Tables 3 and 4). GRB10 is a down-regulated gene associated with one hypo-DMR (19 CpG sites) located in its promoter, which encodes a growth factor receptor-binding protein belonging to the family of adapter proteins that bind receptor tyrosine kinases and signaling molecules, known to interact with insulin-like growth-factor receptors and insulin receptors. GREB1L is an up-regulated gene associated with one hypo-DMR (11 CpG sites) located in its promoter, encoding a poorly characterized protein which displays high similarity with the growth factor Growth Regulating Estrogen Receptor Binding 1, which has, in turn, a role in cell proliferation and in the response to estrogen stimuli. SIM2 is an up-regulated gene associated with two hypo-DMRs in its promoter (4 and 11 CpG sites), encoding a transcription factor that is implicated in the regulation of neurogenesis and the immune response. WARS is a down-regulated gene associated with one hyper-DMR in the promoter region (2 CpG sites), encoding a tryptophanyl-tRNA synthetase that catalyzes the aminoacylation of tRNA with tryptophan, involved in the innate immune response and inflammation. BMP6 gene, encoding a ligand of the TGF-beta (see above), is up-regulated and associated with three hypo-DMRs (8, 4 and 8 CpG sites) and one hyper-DMR (4 CpG sites) in the gene body; two of these hypo-DMRs and the hyper-DMRs are included in an annotated CpG island containing 307 CpG sites that spans the exon 1 of the gene. CDC42EP5 is a down-regulated gene encoding a Rho GTP-ase involved in the polymerization of the actin cytoskeleton, influencing the cellular shape, associated with one hypo-DMR in the gene body (2 CpG sites), which is included in an annotated CpG island containing 681 CpG sites, spanning the 3’-end of the gene. DUOX1 is an up-regulated gene associated with one hypo-DMR in the gene body (14 CpG sites), encoding a member of the NADPH oxidase family known to mediate Reactive Oxygen Species (ROS) production, involved in innate immune defense and anti-microbial function. KCNK10 is an up-regulated gene associated with two hyper-DMRs in the gene body (2 and 4 CpG sites), encoding a factor that regulates calcium homeostasis. KIF21A is a down-regulated gene associated with one hypo-DMR in the gene body (6 CpG sites), encoding a member of the KIF4 sub-family of kinesin-like motor proteins with a possible role in the microtubule-dependent transport. LOC102409914 (or ZNF786) is an up-regulated gene associated with one hypo-DMR (12 CpG sites), encoding a zinc finger protein possibly involved in the transcriptional regulation. PM20D2 is a down-regulated gene associated with one hyper-DMR in the gene body (7 CpG sites), which encodes a peptidase that hydrolyzes dipeptides having basic amino acids lysine, ornithine or arginine, with a possible role in the elimination of dipeptide products of carnosine biosynthesis. SAMD11 is an up-regulated gene associated with two hypo-DMRs in the gene body (9 and 3 CpG sites), encoding a factor predicted to negatively regulate the transcription. TENM3 is a down-regulated gene associated with one hyper-DMR in the gene body (10 CpG sites), which encodes a member of the teneurin transmembrane protein family, with a possible role in the regulation of neuronal development. TFAP2E gene is an up-regulated gene associated with one hyper-DMR in its gene body (6 CpG sites), encoding a transcription factor that binds promoters and enhancers and regulates the expression of genes implicated in the cell cycle, development and differentiation. TNXB is an up-regulated gene associated with one hypo-DMR in its gene body (12 CpG sites), encoding a member of the tenascin family of extracellular matrix glycoproteins, implicated in the expression and deposition of collagen into the extracellular matrix. TRIM7 is an up-regulated gene associated with one hypo-DMR in its gene body (3 CpG sites), which encodes a member of the tripartite motif (TRIM) family with a possible role in the initiation of glycogen biosynthesis. ZFYVE16 is a down-regulated gene associated with a hyper-DMR in its gene body (10 CpG sites), which encodes an endosomal protein belonging to the FYVE zinc finger proteins, which regulate the endosomal membrane trafficking and are involved in the TGF-beta signaling pathway. LOC102400339 (or MAD2L1) is a down-regulated gene associated with five DMRs, two of which are located within its promoter (one hypo-DMR, 6 CpGs, and one hyper-DMR, 7 CpGs) and three of which within gene body (two hypo-DMRs, 14 and 30 CpGs, and one hyper-DMRs, 11 CpGs). LOC102400339 encodes a component of the mitotic spindle involved in the assembly checkpoint. UPP1 is an up-regulated gene associated with one hypo-DMR spanning its promoter and gene body (3 CpG sites), which encodes a uridine phosphorylase with a role in the degradation and salvage of pyrimidine ribonucleosides.
LOC112581380, LOC112587105, LOC102399222, LOC112582971, are still uncharacterized loci and no reliable information is available.
Discussion
The present study aimed to assess the influence of a green forage diet on genomic DNA methylation in the ruminal wall of Italian Mediterranean dairy buffaloes through a genome-wide approach, represented by Reduced Representation Bisulfite Sequencing (RRBS). The latter is a powerful yet cost-efficient method for studying DNA methylation on a genomic scale combining restriction enzymes and bisulfite sequencing to obtain information on DNA methylation in regions of the genome with a high CpG content24. Additionally, we explored the relationship between changes in DNA methylation and gene expression upon administration of a green forage diet, by combining DNA methylation data, here reported, with transcriptomic profiles, published in our previous study obtained from rumen samples of the same animals23. Although several studies conducted in ruminants highlighted that different diet regimens impact DNA methylation21,25, to the best of our knowledge, we underlined for the first time that green forage diet has an impact at the genome-wide level on DNA methylation signature in buffaloes ruminal wall, and to some extent, this has an effect on gene expression.
We found that green forage modulates globally DNA methylation, without a specific genomic preference, with more than 6500 DMRs uniformly distributed along chromosomes. The slight prevalence of hypo-DMRs, with respect to hyper-DMRs, allows to hypothesize a trend of green forage diet to induce DNA hypomethylation. DNA methylation changed in both genes and intergenic regions. However, DMRs localized in genes (in promoters, exons and introns) are mainly hypomethylated, whereas intergenic DMRs are mostly hypermethylated, which suggests a preferential functional effect of green forage on DNA methylation, based on DMR specific genomic location. Among the DMRs associated with genic regions, we highlighted a slight prevalence of DMRs at regulatory regions upstream genes, with respect to exon and introns, according to the enrichment of CpG sites within the promoters26 and their known functional role in these genomic regions.
Our study identified 4648 genes closely associated with one or more DMRs (DMGs), the majority of which are classified as protein-coding genes. These genes encode factors implicated in the regulation of biological processes important for buffalo health and welfare, such as those related to response to oxidative stress, circadian rhythms, fungal infections and inflammation, as well as factors involved in rumen functionality, such as muscle structure/function, and energy, lipid, and amino acid metabolism. Of note, green forage administration resulted to impact similar GO categories in buffalo rumen also when its effect on transcriptional program has been investigated23, reinforcing the idea of a stepwise process by which the administration of green forage might lead to changes in the epigenetic landscape, such as DNA methylation, which, in turn, modulates gene expression. Moreover, our data suggested that green forage administration increases rumen activity, and this evidence allows to hypothesize that this diet might have a role in the production of nutraceutical molecules, the levels of which might be enhanced also in milk. Accordingly, it was previously reported that 30% green forage enhances the concentration of beneficial biomolecules, including betaine and carnitine, in milk and other dairy products in buffaloes [17, 26].
The integration of rumen DNA methylation profiles with transcriptomic data, previously obtained in the ruminal wall of the same animals23, did not show a global correlation between DNA methylation status in the DMRs and gene expression in Treated buffaloes. This result is not entirely unexpected, considering that methylation is not the only epigenetic factor that regulates gene expression, which is known to be influenced also by other epigenetic features, such as histone methylation and/or acetylation, and microRNAs contribution27,28. Furthermore, it should be pointed out that RRBS analysis allows the study of the methylation status of DNA regions enriched in CpG sites, which are not the only regions directly responsible for transcriptional regulation29. We may suppose that genes for which the expression do not correlates with the methylation status are regulated by DNA methylation of other DNA regions that have not been investigated with this study. On the other hand, when we integrated the previously identified differential expression data23 with changes of DNA methylation in genes associated with DMRs upon green forage administration, we identified a discrete number of genes the expression of which varies as a function of methylation in their promoter and/or gene body, this latter DNA regions known to be functionally implicated in the regulation of gene expression, in dependence on their methylation status. Indeed, in mammals, DNA methylation of gene promoter is usually associated with transcriptional repression, whereas promoter hypomethylation leads to enhanced gene expression. On the contrary, DNA methylation of gene body, generally, associates with gene transcription30. However, it should be noted that not all previously identified green forage-modulated genes are coherently associated with changed DNA methylation status in their promoter or gene body and vice-versa. We might speculate that, for these genes, DNA methylation is not the leading event that regulates their transcription. In addition, it is known that even the methylation of a single CpG within a binding site of transcription factors, that could have escaped our analysis due to the technical limitation, theoretically, could influence gene regulation31,32. Likewise, regarding genes associated with changed DNA methylation but not with differential expression, one might speculate that the observed methylation changes in their DMRs are not sufficient to trigger a differential expression, as observed also in other biological contexts22.
Our study underlined that the expression of several genes potentially important for buffalo welfare and rumen functionality, is modulated by green forage administration, probably through changes in DNA methylation. In this frame, we decided to examine the role of a sub-set of them more extensively. SIM2 is one of the most significantly up-regulated genes upon green forage administration in our datasets, associated with two hypomethylated DMRs in its promoter, according to its over-expression. This gene encodes the X1 isoform of the single-minded homolog 2 protein, involved in the antimicrobial response in the gastro-intestinal tract through the regulation of the expression of antimicrobial peptides and defensins33. We might suppose that SIM2 enhanced expression induced by green forage diet strengthens the antimicrobial response in the buffalo ruminal wall. Conversely, WARS gene is down-regulated and associated with a hypermethylated DMR located in the promoter. WARS1 encodes a cytoplasmatic aminoacyl-tRNA synthetase known to have roles in innate immune system and inflammatory responses34,35,36. In a Genome-Wide Association Study (GWAS), genetic variants of WARS2, encoding the mitochondrial form of aminoacyl-tRNA synthetases, were associated with fatty acid composition in beef37. DUOX1, is up-regulated in Treated animals and coherently associated with increased DNA methylation in its gene body. This gene encodes a member of the NADPH oxidase family involved in the production of hydrogen peroxide. Notably, increased expression of DUOX1 in goat rumen was negatively correlated with the microbial metabolism of starch and sucrose, evidence that suggests a role of this factor in rumen immunity and the control of microbial colonization38. Accordingly, DUOX1-mediated ROS production has been proposed to have an antimicrobial effect39. As for DUOX1, also KCNK10 gene is up-regulated upon green forage administration and coherently associated with two hypermethylated DMRs in the gene body. This gene encodes a potassium transport channel activated by arachidonic acid and other naturally occurring unsaturated free fatty acids40, with a function for muscle contraction and differentiation. This evidence suggests a role of green forage in stimulating the ruminal muscular function. The TFAP2E is an up-regulated gene consistently associated with a hypermethylated DMR located in the gene body. The product of this gene is a transcription factor that regulates numerous genes involved in multiple biological processes, such as development, cell cycle and differentiation, through the binding of their promoters and enhancers41. In a recent GWAS study, genetic variants of TFAP2E were associated with ketosis in cattle42. Furthermore, TFAP2E has been indicated among the candidate genes involved in the determination of the genetic trait of fat yield in buffalo based on single-step genomic best linear unbiased prediction and random regression models43. TNXB gene, which encodes a glycoprotein of extracellular matrix with a role in fatty acid metabolism, is down-regulated upon green forage administration and associated with a hypomethylated DMR in the gene body. TNXB down-regulation was observed in Brazilian bull thoracic muscle in association with high levels of fatty acids in the muscle44. It is reasonable to hypothesize that the down-regulation of this gene, observed following a green forage diet, could have an impact on the quality of the meat. Recently, it was reported that different expression levels of TNXB are associated with variations in adipose tissue deposition in different sheep breeds45. These data are in agreement with the role of the extracellular matrix in the regulation of adipogenesis, through the formation of a structural support for adipocytes.
We found changes in both DNA methylation status and in the expression levels of genes encoding factors related to energy metabolism and, hence, rumen fermentability due to the composition of green forage, enriched in simple sugars. Indeed, different genes the product of which participates in the regulation of systemic glucose homeostasis, such as BMP6 and UPP1, have been highlighted in this study. The BMP6 gene, which encodes a key regulator of iron metabolism46 and glucose homeostasis47, is up-regulated in Treated animals and associated with four DMRs located in two annotated CpG islands in the gene body, three of which are hypomethylated and one is hypermethylated. Three of the four DMRs are included in an annotated CpG island that spans the exon 1 and the promoter of the gene. It is interesting to note that the up-regulation of BMP6 in the rumen epithelium had been reported in cattle fed a high grain content diet48, and in Japanese Black cattle during weaning, this latter a period in which rumen papillae start to increase their density and length49. UPP1 gene is up-regulated in treated animals and associated with two hypomethylated DMRs, one located in the promoter and the other one in the gene body. UPP1 encodes an enzyme with a role in the degradation and salvage of pyrimidine ribonucleosides that catalyzes the reversible phosphorylation of uridine to uracil and ribose-1-phosphate. The latter is converted into fructose-6-phosphate and glyceraldehyde-3-phosphate by the non-oxidative branch of the pentose phosphate pathway. According to this, uridine has been reported to be a meaningful source of energy in specific tissues (blood, lung, brain and kidney) serving as an alternative energy source50. Since it has already been reported that fresh forage and hay-based feeding systems resulted in a greater content of uridine in ruminants51, probably due to the activity of the Christensenellaceae_R-7_group52 we may speculate an UPP1 transcriptional up-regulation in response to the higher levels of uridine in Treated animals. UPP1 expression, seems to be also correlated with ruminant health status. Indeed, this gene is down-regulated in rumen tissue under acidotic conditions53, and its lower expression has been reported in peripheral blood of cattle showing high levels of stress54.
The inconsistency between gene expression and methylation status of gene body observed for BMP6, TENM3 and UPP1 genes is not completely unexpected. Although the methylation of the gene body seems to be a feature of active genes, indeed, its role is poorly understood; it was proposed that this phenomenon does not impede transcription, but whether it promotes gene expression is still not clear55. Of note, in breast cancer cells, hypomethylation of the gene body has been associated with gene down-regulation56. Moreover, additional functions have been proposed for this epigenetic feature, such as the regulation of the splicing57.
Conclusions
To the best of our knowledge, our study provides, for the first time, a broad overview of DNA methylation changes that occur in buffaloes following green forage feeding, and suggests that green forage modulates, at least in part, the ruminal wall transcription through changes in this epigenetic signature.
We found that green forage diet has a genome-wide impact on DNA methylation, leading mainly to DNA hypomethylation with respect to hypermethylation. Furthermore, the inspection of DNA methylation changes in the proximity of genes highlighted that these changes regard principally protein-coding genes, many of which encode factors implicated in the regulation of biological processes important for buffalo welfare and rumen functionality.
The integration of changes in rumen DNA methylation profiles with the previously obtained differential gene expression data upon green forage administration deepened our understanding of the response of metabolism, activity, and physiology of rumen to this diet regimen, and pinpoints DNA methylation as a key player in the diet-dependent transcriptional modulation, adding a further piece of knowledge on the potential effects of the diet on molecular mechanisms linked to the animal welfare and production of nutraceutical molecules.
Our results will open new perspectives for the setting of innovative precision feeding-based strategies, taking into account both animal welfare and the production of animal derived food rich in health promoting biomolecules.
Materials and methods
Animals, dietary treatment, and ruminal tissue collection
Animals, dietary treatment, and ruminal tissue collection
All experimental procedures were performed according to the European Directive 2010/63/EU and the Italian Legislative Decree No. 26 dated 4th March 2014 and received institutional approval from the Ethical Animal Care and Use Committee of the University of Naples “Federico II” (Protocol No. 25532 − 2022). The animals involved in this study were the same used for the transcriptomic analyses previously published [23]. The study was carried out over 60 days at a commercial buffalo dairy farm located in southern Italy, by using Italian Mediterranean dairy buffaloes (n 14), for which the slaughter was already scheduled. Buffaloes were reared in pens with a concrete floor and milked twice daily and had undergone an adaptation period of 14 days before beginning the experiment. Animals were randomly assigned to two homogeneous groups (Control and Treated) according to parity (5.6 ± 0.6 vs. 5.7 ± 0.5 in Treated vs. Control group), age (8.2 ± 0.9 vs. 9.2 ± 0.4 in Treated vs. Control group), days in milk (191 ± 37.5 vs. 208 ± 43.0 in Treated vs. Control group), and average milk production (6.5 ± 0.5 vs. 6.0 ± 0.8 kg/day in Treated vs. Control group), whose values were recorded in the ten days before the beginning of the trial. The diet of Control buffaloes was a Total Mixed Ration (TMR) whilst Treated buffaloes received TMR + green forage which comprised ryegrass (approximately 30% of the diet). The forage was re-blossoming stage ryegrass cut twice a day to avoid any fermentation, without storage and administered to animals through mixing wagon. The forage: concentrate ratio was 56:44 and 69:31 respectively in Control and Treated animals. The two diets were created isonitrogenous and isoenergetic and differed only in the inclusion of ryegrass in treated animals (Table 5). Animals were fed in the morning and in the evening. At the conclusion of the trial the animals, which were at the end of their reproductive and productive careers, were euthanized all together by penetrating captive bolt, according to the AVMA Guidelines for the Humane Slaughter of Animals, following procedures approved by the Ethical Animal Care and Use Committee of the University of Naples “Federico II”. After the slaughter, samples of rumen were collected for molecular studies. Each sample was labeled and stored at -80 °C until analysis.
Isolation of genomic DNA from ruminal wall
Genomic DNA was extracted from 100 mg of frozen ruminal wall tissue of 14 dairy buffaloes, of which 6 received TMR (Control group) and 8 were fed a TMR + green forage (Treated group). Genomic DNA was extracted using Wizard Genomic DNA Purification Kit (Promega, Madison, WI, USA), according to the manufacturer’s instructions. As a further purification step, 10 µg of extracted DNA were subjected to an additional proteinase K digestion in 200 mM Tris HCl pH 7.5, 0.5 mM EDTA, 0.5% SDS and 15 µg proteinase K (Roche Basel, Switzerland) for 2 h at 55 °C. The quality and the concentration of DNA samples were determined using NanoDrop (NanoDrop Technologies, Wilmington, DE, USA) and agarose gel electrophoresis. For all samples, the ratios of the absorbance at 260 nm and 280 nm (A260/280) and at 260 nm and 230 nm (A260/230) were between 1.8 and 2.0 and 1.9–2.2, respectively. Five µg of each DNA sample was sent to IGA Technology Services (Udine, Italy) for library generation and sequencing.
RRBS library Preparation and sequencing
After quantification of DNA samples with Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA), 100 ng of high-quality DNA from each sample was used as input for library generation using the Ovation RRBS Methyl-Seq System (NuGEN, Redwood City, CA), following the manufacturer’s instructions. Briefly, DNA was digested with the methylation insensitive restriction enzyme MspI, which recognizes CCGG, ligated to sequencing adaptors, repaired, converted with bisulfite, and amplified by PCR (7 cycles of PCR). Unmethylated lambda phage DNA was spiked-in to estimate the bisulfite conversion rate. The concentration of the final libraries was measured with both Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA) and Caliper GX DNA Assay (PerkinElmer, Waltham, MA). Libraries were then sequenced in paired-end mode of 100 bp, 20 M reads per sample, using a NovaSeq6000 (Illumina, San Diego, CA). The processing of raw data for both format conversion and de-multiplexing was performed by Bcl2Fastq 2.20 version of the Illumina pipeline.
RRBS reads processing and mapping
Bioinformatic analysis were carried out by Sequentia Biotech (Barcelona, Spain). Raw reads were first processed with UMI-tools58 to extract the UMI sequences. After that, trimming was performed to remove low quality bases (minimum Phred = 15), adapter sequences and short reads (minimum length = 35 bp) using Trimmomatic59. Filtered reads were then aligned to the buffalo reference genome sequence (NCBI: GCF_003121395.1, “UOA_WB_1”) using Bismark60. After mapping, UMI-tools dedup was used to remove duplicates based on UMI information. CpGs methylation levels were then extracted from the deduplicated BAM files using MethylDackel [https://github.com/dpryan79/MethylDackel].
Identification of differentially methylated regions (DMRs)
The CpG methylation values obtained from MethylDackel were used as input for the statistical analyses, which were performed with the R programming language. The package bsseq 61 was used to obtain normalized smoothed methylation values with the following options: ns = 20, h = 150, maxGap = 10,000,000. Only CpGs with a coverage of at least 5 reads across all the samples were analyzed. Differentially methylated regions between the two groups were detected using the dmrFinder function from bsseq. Significant DMRs were those having an areaStat <-3.2 and > + 3.2 and a quantile-based cutoff of 5%.
Gene ontology analysis
The gene annotations available in NCBI (NCBI: GCF_003121395.1, “UOA_WB_1”) were used to associated the DMRs to the genomic features of B. bubalis. For operations between BED and GFF files, bedtools was used62. Gene ontology (GO) enrichment analysis of DMR-associated genes (DMGs) was performed using an in-house developed script based on the hypergeometric test of GO distribution, setting an FDR < 0.05 as the threshold for statistical significance. All annotated genes in the genome were used as background for the GO analysis.
Methylome and RNA-seq data integration
The FPKM values from the RNA-seq data produced from the same samples described in this study [23] were downloaded and converted to log2 FPKM. The genes associated to differentially methylated regions (DMG) were split between those having a DMR in the promoter region and those in the gene body. Pearson correlation values were then calculated between the percentage of methylation of the DMRs and the log2 FPKM of the corresponding genes. In addition, the list of DMG and DEGs, from [23], were intersected and a scatter plot was generated to show the relationship between the log2 fold change values from the RNA-seq and the difference in methylation. Only the significantly differentially expressed genes associated with a significantly differentially methylated region were highlighted in the plot. Plots and data analysis were carried out with R using ggplot2 (https://ggplot2.tidyverse.org/) for visualization and the corr.test function from the psych package (https://rdrr.io/cran/psych/), for the correlation analysis.
Data availability
The RRBS dataset supporting the conclusions of this article is available in the NCBI repository (BioProject: PRJNA1107499).
References
Deb, G. K., Nahar, T. N., Duran, P. G. & Presicce, G. A. Safe and sustainable traditional production: the water Buffalo in Asia. Front. Environ. Sci. 4 https://doi.org/10.3389/fenvs.2016.00038 (2016).
Dash, S. et al. Effect of heat stress on reproductive performances of dairy cattle and buffaloes: A review. Vet. World. 9, 235–244. https://doi.org/10.14202/vetworld.2016.235-244 (2016).
El Debaky, H. A. et al. Potential of water Buffalo in world agriculture: challenges and opportunities. Appl. Anim. Sci. 35, 255–268. https://doi.org/10.15232/aas.2018-01810 (2019). Review.
Hamid, M. A., Siddiky, M. N. A., Rahman, M. A. & Hossain, K. M. Scopes and opportunities of Buffalo farming in Bangladesh: a review. SAARC J. Agric. 14, 63–77 (2016).
Salzano, A. et al. Breed and feeding system impact the bioactive Anti-Inflammatory properties of bovine milk. Int. J. Mol. Sci. 23 https://doi.org/10.3390/ijms231911088 (2022).
Salzano, A. et al. Green feed increases antioxidant and antineoplastic activity of Buffalo milk: A globally significant livestock. Food Chem. 344, 128669. https://doi.org/10.1016/j.foodchem.2020.128669 (2021).
Matthews, C. et al. The rumen microbiome: a crucial consideration when optimising milk and meat production and nitrogen utilisation efficiency. Gut Microbes. 10, 115–132. https://doi.org/10.1080/19490976.2018.1505176 (2019).
Schwab, C. G. & Broderick, G. A. A 100-Year review: protein and amino acid nutrition in dairy cows. J. Dairy. Sci. 100, 10094–10112. https://doi.org/10.3168/jds.2017-13320 (2017).
Flint, H. J. & Bayer, E. A. Plant cell wall breakdown by anaerobic microorganisms from the mammalian digestive tract. Ann. N Y Acad. Sci. 1125, 280–288. https://doi.org/10.1196/annals.1419.022 (2008).
Sun, H. Z. et al. Multi-omics reveals functional genomic and metabolic mechanisms of milk production and quality in dairy cows. Bioinformatics 36, 2530–2537. https://doi.org/10.1093/bioinformatics/btz951 (2020).
Baaske, L., Gabel, G. & Dengler, F. Ruminal epithelium: a checkpoint for cattle health. J. Dairy. Res. 87, 322–329. https://doi.org/10.1017/S0022029920000369 (2020).
Chaudry, A. S. Forage based animal production systems and sustainability, an invited keynote. Revista Brasileira De Zootecnia. 37, 78–84. https://doi.org/10.1590/S1516-35982008001300010 (2008).
Elgersma, A., Soegaard, K. & Jensen, S. K. Fatty acids, alpha-tocopherol, beta-carotene, and lutein contents in forage legumes, Forbs, and a grass-clover mixture. J. Agric. Food Chem. 61, 11913–11920. https://doi.org/10.1021/jf403195v (2013).
Eden, A., Gaudet, F., Waghmare, A. & Jaenisch, R. Chromosomal instability and tumors promoted by DNA hypomethylation. Science 300, 455. https://doi.org/10.1126/science.1083557 (2003).
Li, E., Beard, C. & Jaenisch, R. Role for DNA methylation in genomic imprinting. Nature 366, 362–365. https://doi.org/10.1038/366362a0 (1993).
Li, E., Bestor, T. H. & Jaenisch, R. Targeted mutation of the DNA methyltransferase gene results in embryonic lethality. Cell 69, 915–926. https://doi.org/10.1016/0092-8674(92)90611-f (1992).
Lorincz, M. C., Dickerson, D. R., Schmitt, M. & Groudine, M. Intragenic DNA methylation alters chromatin structure and elongation efficiency in mammalian cells. Nat. Struct. Mol. Biol. 11, 1068–1075. https://doi.org/10.1038/nsmb840 (2004).
Cai, D. et al. Betaine supplementation in maternal diet modulates the epigenetic regulation of hepatic gluconeogenic genes in neonatal piglets. PLoS One. 9, e105504. https://doi.org/10.1371/journal.pone.0105504 (2014).
Hu, Y. et al. In Ovo injection of betaine affects hepatic cholesterol metabolism through epigenetic gene regulation in newly hatched chicks. PLoS One. 10, e0122643. https://doi.org/10.1371/journal.pone.0122643 (2015).
Murdoch, B. M., Murdoch, G. K., Greenwood, S. & McKay, S. Nutritional influence on epigenetic marks and effect on livestock production. Front. Genet. 7, 182. https://doi.org/10.3389/fgene.2016.00182 (2016).
Haluskova, J., Holeckova, B. & Stanicova, J. DNA methylation studies in cattle. J. Appl. Genet. 62, 121–136. https://doi.org/10.1007/s13353-020-00604-1 (2021).
Li, Y. et al. DNA methylation, MicroRNA expression profiles and their relationships with transcriptome in grass-fed and grain-fed Angus cattle rumen tissue. PLoS One. 14, e0214559. https://doi.org/10.1371/journal.pone.0214559 (2019).
Salzano, A. et al. Transcriptomic profiles of the ruminal wall in Italian mediterranean dairy buffaloes fed green forage. BMC Genom. 24, 133. https://doi.org/10.1186/s12864-023-09215-6 (2023).
Xi, Y. et al. RRBSMAP: a fast, accurate and user-friendly alignment tool for reduced representation bisulfite sequencing. Bioinformatics 28, 430–432. https://doi.org/10.1093/bioinformatics/btr668 (2012).
Zhang, N. Epigenetic modulation of DNA methylation by nutrition and its mechanisms in animals. Anim. Nutr. 1, 144–151. https://doi.org/10.1016/j.aninu.2015.09.002 (2015).
Saxonov, S., Berg, P. & Brutlag, D. L. A genome-wide analysis of CpG dinucleotides in the human genome distinguishes two distinct classes of promoters. Proc. Natl. Acad. Sci. U S A. 103, 1412–1417. https://doi.org/10.1073/pnas.0510310103 (2006).
Fang, L. et al. Functional annotation of the cattle genome through systematic discovery and characterization of chromatin States and butyrate-induced variations. BMC Biol. 17 https://doi.org/10.1186/s12915-019-0687-8 (2019).
Zhang, L., Lu, Q. & Chang, C. Epigenetics in health and disease. Adv. Exp. Med. Biol. 1253, 3–55. https://doi.org/10.1007/978-981-15-3449-2_1 (2020).
Yong, W. S., Hsu, F. M. & Chen, P. Y. Profiling genome-wide DNA methylation. Epigenetics Chromatin. 9, 26. https://doi.org/10.1186/s13072-016-0075-3 (2016).
Moore, L. D., Le, T. & Fan, G. DNA methylation and its basic function. Neuropsychopharmacology 38, 23–38. https://doi.org/10.1038/npp.2012.112 (2013).
Furst, R. W., Kliem, H., Meyer, H. H. & Ulbrich, S. E. A differentially methylated single CpG-site is correlated with Estrogen receptor alpha transcription. J. Steroid Biochem. Mol. Biol. 130, 96–104. https://doi.org/10.1016/j.jsbmb.2012.01.009 (2012).
Tsuboi, K. et al. Single CpG site methylation controls Estrogen receptor gene transcription and correlates with hormone therapy resistance. J. Steroid Biochem. Mol. Biol. 171, 209–217. https://doi.org/10.1016/j.jsbmb.2017.04.001 (2017).
Chen, K. J., Lizaso, A. & Lee, Y. H. SIM2 maintains innate host defense of the small intestine. Am. J. Physiol. Gastrointest. Liver Physiol. 307, G1044–1056. https://doi.org/10.1152/ajpgi.00241.2014 (2014).
Ahn, Y. H. et al. Secreted tryptophanyl-tRNA synthetase as a primary defence system against infection. Nat. Microbiol. 2, 16191. https://doi.org/10.1038/nmicrobiol.2016.191 (2016).
Lee, H. C. et al. Released Tryptophanyl-tRNA synthetase stimulates innate immune responses against viral infection. J. Virol. 93 https://doi.org/10.1128/JVI.01291-18 (2019).
Nguyen, T. T. T. et al. Tryptophan-dependent and -independent secretions of tryptophanyl- tRNA synthetase mediate innate inflammatory responses. Cell. Rep. 42, 111905. https://doi.org/10.1016/j.celrep.2022.111905 (2023).
Cesar, A. S. et al. Genome-wide association study for intramuscular fat deposition and composition in Nellore cattle. BMC Genet. 15, 39. https://doi.org/10.1186/1471-2156-15-39 (2014).
Pan, X. et al. Dynamics of rumen gene expression, Microbiome colonization, and their interplay in goats. BMC Genom. 22, 288. https://doi.org/10.1186/s12864-021-07595-1 (2021).
Rada, B. & Leto, T. L. Oxidative innate immune defenses by Nox/Duox family NADPH oxidases. Contrib. Microbiol. 15, 164–187. https://doi.org/10.1159/000136357 (2008).
Bang, H., Kim, Y. & Kim, D. TREK-2, a new member of the mechanosensitive tandem-pore K + channel family. J. Biol. Chem. 275, 17412–17419. https://doi.org/10.1074/jbc.M000445200 (2000).
Tummala, R., Romano, R. A., Fuchs, E. & Sinha, S. Molecular cloning and characterization of AP-2 epsilon, a fifth member of the AP-2 family. Gene 321, 93–102. https://doi.org/10.1016/s0378-1119(03)00840-0 (2003).
Soares, R. A. N., Vargas, G., Duffield, T., Schenkel, F. & Squires, E. J. Genome-wide association study and functional analyses for clinical and subclinical ketosis in Holstein cattle. J. Dairy. Sci. 104, 10076–10089. https://doi.org/10.3168/jds.2020-20101 (2021).
Lazaro, S. F. et al. Genomic studies of milk-related traits in water Buffalo (Bubalus bubalis) based on single-step genomic best linear unbiased prediction and random regression models. J. Dairy. Sci. 104, 5768–5793. https://doi.org/10.3168/jds.2020-19534 (2021).
Berton, M. P. et al. Gene expression profile of intramuscular muscle in Nellore cattle with extreme values of fatty acid. BMC Genom. 17 https://doi.org/10.1186/s12864-016-3232-y (2016).
Li, B. et al. Transcriptome analysis of adipose tissues from two fat-tailed sheep breeds reveals key genes involved in fat deposition. BMC Genom. 19, 338. https://doi.org/10.1186/s12864-018-4747-1 (2018).
Andriopoulos, B. Jr. et al. BMP6 is a key endogenous regulator of Hepcidin expression and iron metabolism. Nat. Genet. 41, 482–487. https://doi.org/10.1038/ng.335 (2009).
Pauk, M. et al. A novel role of bone morphogenetic protein 6 (BMP6) in glucose homeostasis. Acta Diabetol. 56, 365–371. https://doi.org/10.1007/s00592-018-1265-1 (2019).
Zhao, K., Chen, Y. H., Penner, G. B., Oba, M. & Guan, L. L. Transcriptome analysis of ruminal epithelia revealed potential regulatory mechanisms involved in host adaptation to gradual high fermentable dietary transition in beef cattle. BMC Genom. 18, 976. https://doi.org/10.1186/s12864-017-4317-y (2017).
Nishihara, K. et al. Comparative transcriptome analysis of rumen papillae in suckling and weaned Japanese black calves using RNA sequencing. J. Anim. Sci. 96, 2226–2237. https://doi.org/10.1093/jas/skx016 (2018).
Skinner, O. S. et al. Salvage of ribose from uridine or RNA supports Glycolysis in nutrient-limited conditions. Nat. Metab. 5, 765–776. https://doi.org/10.1038/s42255-023-00774-2 (2023).
Rocchetti, G., Gallo, A., Nocetti, M., Lucini, L. & Masoero, F. Milk metabolomics based on ultra-high-performance liquid chromatography coupled with quadrupole time-of-flight mass spectrometry to discriminate different cows feeding regimens. Food Res. Int. 134, 109279. https://doi.org/10.1016/j.foodres.2020.109279 (2020).
Ma, X. et al. Multi-omics revealed the effects of dietary energy levels on the rumen microbiota and metabolites in Yaks under house-feeding conditions. Front. Microbiol. 14, 1309535. https://doi.org/10.3389/fmicb.2023.1309535 (2023).
Enemark, J. M. The monitoring, prevention and treatment of sub-acute ruminal acidosis (SARA): a review. Vet. J. 176, 32–43. https://doi.org/10.1016/j.tvjl.2007.12.021 (2008).
Del Corvo, M. et al. Genome-Wide DNA methylation and gene expression profiles in cows subjected to different stress level as assessed by cortisol in milk. Genes (Basel). 11 https://doi.org/10.3390/genes11080850 (2020).
Jones, P. A. Functions of DNA methylation: Islands, start sites, gene bodies and beyond. Nat. Rev. Genet. 13, 484–492. https://doi.org/10.1038/nrg3230 (2012).
Hon, G. C. et al. Global DNA hypomethylation coupled to repressive chromatin domain formation and gene Silencing in breast cancer. Genome Res. 22, 246–258. https://doi.org/10.1101/gr.125872.111 (2012).
Laurent, L. et al. Dynamic changes in the human methylome during differentiation. Genome Res. 20, 320–331. https://doi.org/10.1101/gr.101907.109 (2010).
Smith, T., Heger, A. & Sudbery, I. UMI-tools: modeling sequencing errors in unique molecular identifiers to improve quantification accuracy. Genome Res. 27, 491–499. https://doi.org/10.1101/gr.209601.116 (2017).
Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for illumina sequence data. Bioinformatics 30, 2114–2120. https://doi.org/10.1093/bioinformatics/btu170 (2014).
Krueger, F. & Andrews, S. R. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics 27, 1571–1572. https://doi.org/10.1093/bioinformatics/btr167 (2011).
Hansen, K. D., Langmead, B. & Irizarry, R. A. BSmooth: from whole genome bisulfite sequencing reads to differentially methylated regions. Genome Biol. 13, R83. https://doi.org/10.1186/gb-2012-13-10-r83 (2012).
Quinlan, A. R. & Hall, I. M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841–842. https://doi.org/10.1093/bioinformatics/btq033 (2010).
Kilkenny, C., Browne, W. J., Cuthill, I. C., Emerson, M. & Altman, D. G. Improving bioscience research reporting: the ARRIVE guidelines for reporting animal research. PLoS Biol. 8, e1000412. https://doi.org/10.1371/journal.pbio.1000412 (2010).
Acknowledgements
The authors are grateful to M. D’Esposito, prematurely passed away, for giving relevant advice in the conception of the study.
Funding
This research was financially supported by Italian Ministry for Economic Development (grant number F/200016/01–03/X45) through a research project entitled “CAPSULE - Ottimizzazione delle teCniche di Allevamento e dei Processi produttivi del Settore lattiero-caseario bUfalino e del vino per la produzione di aLimEnti funzionali, and Agritech National Research Center and European Union Next-GenerationEU (PIANO NAZIONALE DI RIPRESA E RESILIENZA (PNRR) - MISSIONE 4 COMPONENTE 2, INVESTIMENTO 1.4 - D.D. 1032 17/06/2022, CN00000022). This manuscript reflects only the authors’ views and opinions, neither the European Union nor the European Commission can be considered responsible for them.
Author information
Authors and Affiliations
Contributions
G.C. conceived the study; F.D.R., S.F., A.S. and G.C. planned and designed the experiments; G.B. managed the animals and organized the collection of the samples; S.F. and A.S. performed the experiments; R.A.C. and S.F. performed bioinformatics analysis; G.C., F.D.R., A.S., and S.F. wrote the manuscript; G.B. and R.A.C. gave conceptual advice and edited the manuscript; all of the authors have reviewed and approved the final manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Ethics approval and consent to participate
The study was performed according to the ARRIVE guidelines for reporting animal research63. All experimental procedures were carried out in accordance with the European Directive 2010/63/EU and the Italian Legislative Decree No. 26 dated 4th March 2014. The Ethical Animal Care and Use Committee of the University of Naples “Federico II” approved the research (Protocol No. 25532 − 2022).
The informed consent for the use of animals was provided by the cooperating commercial buffalo dairy farm located in southern Italy.
Additional information
Publisher’s note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Electronic supplementary material
Below is the link to the electronic supplementary material.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
About this article
Cite this article
Fioriniello, S., Salzano, A., Bifulco, G. et al. Green forage impacts on the DNA methylation in the ruminal wall of Italian mediterranean dairy buffaloes. Sci Rep 15, 8074 (2025). https://doi.org/10.1038/s41598-025-91969-y
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41598-025-91969-y
