DNA Imprinting and Differentially Expressed Genes in Longissimus Thoracis Muscle of Bos indicus Submitted to Early Weaning Management

1. Introduction

In cattle, early weaning consists of separating the calves from cows at less than 150 days of age [1] and is aimed at reducing the stress of suckling and nutritional requirements of the cow, thus enabling early body condition recovery and return to estrus [2]. However, attention must be paid to early weaned calves in order to avoid malnutrition [3,4]. Furthermore, energy supplementation of early weaned calves can be used to enhance characteristics of interest in animal production [5,6]. In breeds of European origin, early weaning has been shown to exert a major effect on the physiology of calves destined for slaughter when compared to those conventionally weaned at 210 days of age [6,7].

Studies have shown that early weaned calves supplemented with energy sources can achieve superior performance in systems that employ a slaughter age between 18 and 24 months, such as higher body weight, greater subcutaneous fat thickness, greater Longissimus thoracis muscle area, and higher intramuscular fat levels [5,6,8]. Among the key mechanisms underlying these metabolic changes are epigenetic modifications that can affect gene modulation and expression in response to environmental and nutritional stimuli [9] and, consequently, phenotypes.

DNA methylation, which is an important epigenetic modification, can alter gene expression through the formation of a physical barrier that inhibits the binding of general transcription factors or affects the state of chromatin condensation; in addition, DNA methylation is involved in the recruitment of methyl-binding proteins that form repressive complexes [10,11]. The diet is known to play a fundamental role in the determination of the DNA methylation pattern by providing substrates necessary for methylation and cofactors that modulate the enzymatic activity of DNA methyltransferase, as well as in the activity of enzymes that regulate one-carbon metabolism [12]. S-adenosylmethionine (SAM), which is synthesized in the methionine cycle from various dietary precursors, is a potent aminopropyl and methyl donor within the one-carbon cycle and serves as the main substrate for the methylation of DNA, associated proteins, and RNA [13,14].

Epigenetic programming is sensitive to nutritional factors, especially during the prenatal period and early phases of life of animals [12,15]. Thus, nutritional alterations can induce “metabolic imprinting”, a mechanism whereby adequate nutrition in the first months after birth alters glucose metabolism and causes long-term beneficial effects on the health and performance of animals [16]. This mechanism proposes critical periods of development, known as “nutritional windows”, in which changes in the DNA methylation pattern and gene expression can permanently modify the animal’s metabolism [17].

At present, there are no omics studies that have associated epigenetics, particularly DNA methylation, with gene expression in Nellore cattle submitted to early weaning and energy supplementation. In this context, we examined the effects of distinct nutritional management protocols on DNA methylation patterns in the Longissimus thoracis muscle of Nelore calves during the perinatal phase, and how these epigenetic modifications may have influenced the gene expression profile.

2. Results

2.1. Reduced-Representation Bisulfite Sequencing

A mean number of 30 million 150-bp reads were generated for the EW group and 28 million 150-bp reads for the CW group. After quality control and adapter removal, the raw number of reads decreased by approximately 1.6% in both groups, corresponding to the number of clean reads. The overall read mapping rate was 76.45% and 76.55% in the EW and CW groups, respectively (Table 1).

2.2. DNA Methylation Patterns in Muscle

A mean percentage of 3.84% mCpG, 0.54% mCHG and 0.47% mCHH was detected in the genome of EW animals, while these percentages were 3.79%, 0.51% and 0.44%, respectively, in the CW group (Table 2).

We found 481 DMRs showing differences in methylation between the two weaning groups at 205 days (qvalue < 0.05 and methylation difference >= 20%) (Supp. Table S1). Balanced ratios of hypo- and hypermethylated regions were observed in the two groups (CW vs. EW). However, there was a slightly higher level of methylation in early weaned animals, with approximately 52% (250) of hypermethylated regions and 48% (231) of hypomethylated regions (Figure 1).

2.3. Distribution of Differentially Methylated Regions across the Genome

Regarding functional location, DMRs were most frequently present in intergenic regions (65%), followed by introns (14%), exons (13%), and promoters/TSSes (7%) (Figure 2a). Approximately 23% of DMRs were located in CpG islands, while 26% were present in shores (Figure 2b; Supp. Table S1).

There were 172 DMRs overlapping 53 genes, some of them in more than one window and at different functional locations (Supp. Table S2). Of these prospected regions, 80 DMRs were hypermethylated in the EW group, while the other 92 were hypomethylated. These DMRs overlapped 53 different genes, which were classified as DMGs. Regarding functional location, 16 DMRs were located in gene promoters, while three others were present in TSSes (Table 3).

2.4. Functional Enrichment Analysis

Functional enrichment analysis considering the 53 DMGs/DEGs between the EW and CW groups highlighted the following biological processes related to 3’-5’-cyclic adenosine monophosphate (cAMP) signaling and positive lipid regulation: cAMP metabolic process (GO:0046058), cAMP biosynthetic process (GO:0006171), and positive regulation of lipid metabolic process (GO:0045834). Regarding metabolic pathways (KEGG), growth-related pathways can be highlighted, such as growth hormone synthesis, secretion and action and GnRH signaling and others (Table 4).

2.5. Relationship between Methylation and Gene Expression

In general, a trend can be seen towards a relationship between gene expression and levels of intragenic DNA methylation or methylation in promoter regions of the respective genes (Figure 3, Supp. Table S3). The estimated regression coefficient was -0.00316, with a pvalue of 0.0002. The significant slope indicates that many genes harbor hypermethylated regions associated with a decline in expression and hypomethylated regions associated with high gene expression. However, different associations have been reported in the literature in which highly methylated gene bodies tend to exhibit high transcription activity [18]. The ADCY3 (adenylate cyclase 3), STUB1 (STIP1 homology and U-box containing protein 1), and ZNF557 (zinc finger protein 557) genes contained low proportions of 5mC sites in the promotor regions of genes upregulated in the EW group. The same was observed in the CW group for the TECR (trans-2,3-enoyl-CoA reductase), GMDS (GDP-mannose 4,6-dehydratase), SSPO (SCO-spondin), SULF2 (sulfatase 2), MACROD1 (mono-ADP ribosylhydrolase 1), ARID3B (AT-rich interaction domain 3B), and ARHGEF10L (rho guanine nucleotide exchange factor 10 like) genes (Supp. Table S3).

Since not all genes involved in a pathway have their regulation mediated by DNA methylation but can be affected by genes regulated by these mechanisms, the DEG data are important for elucidating the activity of metabolic pathways. These results are therefore essential for studies on epigenetic regulation by DNA methylation and are presented in Supp. Table S4.

3. Discussion

In the present study, the comparison of groups of Nellore calves submitted to EW and energy supplementation and calves submitted to CW revealed an important correlation between methylated DNA regions and gene expression, with many of the genes being associated with fat metabolism. It is known that the epigenome is sensitive to nutritional factors [19]. Comparison of the EW vs CW groups showed a low percent variation in methylation. This finding can be explained by the conserved methylation profile of cattle [19], as well as by the fact that the animals were exposed to the same management conditions from birth to 120 days of age. However, the exposure of EW animals to a period of “nutritional window” (121 to 205 days) under energy supplementation exerted an epigenetic effect, with the identification of 485 DMRs. This finding might be related to the presence of so-called epigenetic modifiers – molecules that are important methyl donors and/or cofactors in DNA methylation reactions such as methionine, choline, folate, and vitamin B12 [20]. Thus, the hypermethylation of genes during the postnatal metabolic programming that occurred in the EW group may be the result of nutritional factors. Since these animals received a high-protein diet (20%), the availability of the compounds mentioned above, which are recognized as methyl donors for DNA methylation, was greater. These compounds enter one-carbon metabolism at different sites but are ultimately converted to the methyl group donor SAM, which provides one methyl group for DNA methylation [21]. In one-carbon metabolism, most methyl groups are derived from choline (60%), followed by methionine (20%) and folate (10-20%) [22]. Methionine is an essential amino acid that is regenerated from homocysteine. In addition to its importance for protein synthesis, when converted to SAM, methionine participates in more than 100 methylation reactions in the body [23]. Thus, variations in dietary methionine may affect DNA methylation [24]. However, some studies suggest that there is no simple correlation between the concentration of methyl donors and DNA methylation and that complex mechanisms of competitive inhibition of DNA methyltransferases may contribute to DNA imprinting [12,25].

Bos indicus breeds are known for their lower meat tenderness and intramuscular fat deposition when compared to Bos taurus [26]. Thus, the identification of gene pathways that could improve fat deposition and meat tenderness in Bos indicus is an important strategy for this animal production sector [26]. Within this context, studies of European cattle breeds suggested that early weaning management leads to changes in DNA methylation patterns and in the expression of genes involved in metabolism [6,8]. In the study by [6], metabolic programming due to early weaning resulted in significant changes in carcass fat deposition, suggesting that early postnatal metabolic imprinting events can be explored as a management tool to improve the value and quality of meat based on increased marbling. In our study, EW animals exhibited 53 DMGs/DEGs. These results are consistent with the effects of early postnatal metabolic programming on the expression of genes levels involved in numerous biological processes, including those related to the positive regulation of lipid metabolic processes.

Among the main genes exhibiting significant epigenetic changes and alteration in expression due to early weaning, we highlight the ADCY3 and AKT1 genes. The ADCY3 gene is involved in the regulation of different processes, including embryogenesis, hormone secretion, glycogen degradation, smooth muscle relaxation, cardiac muscle contraction, and olfaction [27,28,29]. In general, adenylyl cyclases (ADCYs) are enzymes that catalyze the synthesis of cAMP from adenosine triphosphate (ATP) [30]. The present results showed that the cAMP signaling pathway was enriched in the EW group. Cyclic AMP is a second messenger involved in intracellular signal transduction that is associated with the function of kinases in different biochemical processes, including the regulation of carbohydrate and lipid metabolism [30,31]. Increases in cAMP ultimately activate a number of enzymes and transcription factors that regulate satiety and metabolism [32].

In this scenario, we evaluated some pathways linked to energy metabolism and regulated by overexpression of ADCY3 regarding genes differentially expressed in the EW group. The altered expression of ADCY3 is related to an increase in protein kinase cAMP-activated catalytic subunit alpha (PRKACA). This kinase signals changes in adipogenesis, thermogenesis, fatty acid oxidation, and insulin resistance pathways [28]. However, the evaluation of DEGs in this study did not provide solid data on changes in adipogenesis, thermogenesis, or fatty acid oxidation; however, genes involved in de novo lipogenesis were upregulated.

Adipogenesis consists of the differentiation of mesenchymal stem cells into preadipocytes and the proliferation of these cells, the differentiation of preadipocytes into adipocytes, and the conversion to lipid-assimilating cells [33,34,35]. The mechanisms that control adipogenesis in skeletal muscle are largely regulated by transcription factors of the enhancer binding protein (CEBP) and peroxisome proliferator-activated receptor γ (PPARg) family [35,36,37]. The CEBP and PPARg genes were found to be overexpressed but not differentially methylated in the EW group. However, various other genes that are important for adipogenesis did not show consistent changes in expression, including bone morphogenetic protein 4 (BMP4), PR domain containing protein 16 (PRDM16), fibroblast growth factor 10 (FGF10), and fatty acid binding protein 4 adipocyte (FABP4). For example, PRDM16 acts as a switch between the myogenic lineage and brown adipocytes [38,39]. Cells expressing PRDM16 do not undergo differentiation to the myogenic lineage. Thus, although the EW group exhibited a significant change in ADCY3, the alteration in the expression of genes involved in adipogenesis was less consistent. In fact, it has been proposed that adipogenesis starts in the second half of gestation in ruminants and reaches is highest activity during the perinatal period [35]. In the present study, the postnatal period was the nutritional window, i.e., after the main phase of adipocyte differentiation.

Beta-oxidation and thermogenesis are other processes that can be influenced by excess ADCY3 and PRKACA [28]. Although the expression of these regulators was altered in the EW group, no significant changes were observed in the expression of the hormone-sensitive lipase (HSL or LIPE) or carnitine palmitoyltransferase (CPT1 and CPT2) gene, indicating the lack of a difference in beta-oxidation between groups. Cyclic AMP-responsive element-binding protein 1 (CREB1) is regulated by PRKACA during thermogenesis. Likewise, neither the expression of CREB1 nor of other genes involved in thermogenesis was significantly altered in the EW group.

The insulin resistance pathway has also been suggested to respond to cAMP signals [28]. According to [6], an increase in gluconeogenic substrates for ruminal fermentation and in the concentration of insulin in response to circulating glucose leads to greater cellular uptake and utilization of glucose for fat deposition in intramuscular adipocytes proliferated by imprinting, increasing marbling during the growth phase. The gene expression data obtained for the EW group showed increased expression of protein kinase AMP-activated catalytic subunit alpha 1 (PRKAA1 or AMPK), an important mediator of the insulin resistance pathway, along with phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit beta (PI3KCB) and AKT serine/threonine kinase 1 (AKT1) [28]. AKT1 is responsible for the regulation of glucose uptake by mediating the insulin-induced translocation of solute carrier family 2 member 4 or glucose transporter type 4 (SLC2A4 or GLUT4) to the cell surface [40]. AKT1 was hypermethylated in the promoter and intron regions in the EW group; however, no significant alteration in gene expression was observed. Furthermore, the analysis of DEGs in this study revealed no alteration.

AKT signaling has been reported to prevent the degradation of sterol regulatory element binding transcription factor 1 (SREBP1) and thus to promote de novo lipid synthesis [41]. This gene belongs to a family of transcription factors [41] and encodes an important transcriptional regulator of genes involved in lipid synthesis and metabolism, including fatty acid synthase (FASN) and ATP-citrate lyase (ACLY) [42,43]. SREBP1 was overexpressed, demonstrating that AKT1 hypermethylation did not compromise gene expression. Other factors may be involved in the positive regulation of de novo lipid synthesis in the EW group.

Indeed, the DEG data demonstrated overexpression of genes involved in the de novo lipid synthesis pathway, including acetyl-CoA synthetase (ACSS2), ACLY, acetyl-CoA carboxylase (ACACA), stearoyl-CoA desaturase (SCD), and FASN, as well as the transcription factor SREBF-1, as also reported by [44]. An intriguing question raised here is how the postnatal metabolic programming performed in this study has stimulated lipogenesis. According to the review by [35], fat synthesis requires the incorporation of triglycerides into the adipose tissue of the animal after absorption of dietary fatty acids or de novo synthesis of other fatty acids. Acetate – a volatile fatty acid produced during ruminal fermentation – is the main precursor molecule for the synthesis of fatty acids in ruminants [45]. After ruminal absorption and already in the adipocyte, acetyl-CoA synthetase (ACSS2) converts acetate into acetyl-CoA [46] which, in the presence of reduced nicotinamide adenine dinucleotide phosphate (NADPH), is crucial for de novo fatty acid synthesis [35].

Interestingly, the cAMP signaling pathway enriched in this study regulated various other processes through ADCYs. One such process involves the sodium/potassium-transporting ATPase subunits gamma, alpha-1 and beta-2 (KEGG Bos taurus: 281773). These proteins, in turn, participate in digestion and absorption pathways (bta04973 and bta04974) in the intestinal epithelium, including acetate absorption in colonic epithelial cells. Although the expression of ATPases was not obtained for intestinal cells in this study, indirect evidence of increased acetate absorption and de novo lipid synthesis stems from the overexpression of FABP4. Expressed mainly in adipocytes, the product of the FABP4 gene is a cytoplasmic protein that binds to long-chain fatty acids and other hydrophobic ligands and plays a role in uptake, transport, and hydrolysis [47,48]. Increased expression of this gene has been associated with greater intramuscular fat deposition, carcass weight, and fatty acid composition in beef cattle [44,49,50].

Taken together, the data reinforce the proposals of greater fat deposition in cattle submitted to early weaning [6,8,44]. Furthermore, emerging evidence indicates that postnatal metabolic programming through epigenetic modification mainly exerts an effect on the ADCY3 receptor which, in turn, leads to numerous changes in pathways linked to cAMP signaling.

Authors should discuss the results and how they can be interpreted from the perspective of previous studies and of the working hypotheses. The findings and their implications should be discussed in the broadest context possible. Future research directions may also be highlighted.

4. Materials and Methods

4.1. Animals and Collection of Muscle Tissue

The field experiment was conducted in the town of Cáceres, State of Mato Grosso, Brazil. The climate of the region is characterized as tropical, with two well-defined seasons (dry in winter and wet in summer), an average annual temperature of 25°C, and average precipitation of 1,396 mm/year [51]. Forty purebred Nellore cows that met the following criteria were selected from a herd of 3,500 animals: same age, same number of calvings, and having given birth to male calves. In addition, cows with a maximum difference between calvings of 15 days were chosen.

These cows and their calves were kept in a rotational grazing system of Brachiaria decumbens paddocks, with minerals being available ad libitum. At 120 days after the average calving date, 20 calves were weaned and started to receive a concentrate-based supplementation in a new paddock consisting of Brachiaria decumbens where they remained until 205 days of age. These calves formed the early weaning (EW) group. The other 20 calves remained with the cows until 205 days of age, when they were weaned conventionally (CW group) (Table 5). Thus, except for the period between 120 and 205 days corresponding to the different treatments, the animals shared a common environment and received the same management and diet.

For molecular analyses, biopsies of Longissimus thoracis muscle performed between the 12th and 13th rib were obtained from 40 calves (20 EW and 20 CW) at 205 days of age. The collected biopsies were kept in liquid nitrogen and transported to the Laboratory of Genetics and Animal Breeding (LGMA) of the School of Veterinary Medicine and Animal Science (FMVZ), São Paulo State University (Unesp), Botucatu. These samples were stored frozen at -80oC until the time of RNA and DNA extraction.

4.2. Genomic DNA Extraction, RRBS, Quality Control and Alignment to Reference Genome

For DNA methylation analysis, total genomic DNA was extracted from the muscle biopsies of five individuals from each group (EW and CW) for reduced-representation bisulfite sequencing (RRBS), which focuses on genomic regions rich in CpG. For this purpose, 10 genomic libraries were prepared and sequenced (five EW vs five CW) using 30 U of the MspI restriction enzyme to digest CpG-rich regions, followed by bisulfite treatment using the EZ DNA Methylation-Lightning™ Kit (Zymo Research, USA). Sequencing was performed using 150 paired-end reads on the Illumina NovaSeq 6000 platform (Illumina, USA).

The reads obtained were identified using standard base calling software (Illumina, USA) and the quality of the raw fastq files was analyzed with FastQC 0.11.9 [52]. Adapters was trimmed and low-quality reads were removed using the TrimGalore 0.6.4 software [53]. Alignment to the reference genome (Bos taurus ARS-UCD1.2) was performed using the Bismark 0.22.3 software [54]. Total methylated (C) and unmethylated (T) cytosine reads were called using the MethylDackel v0.5.0 tool (https://github.com/dpryan79/MethylDackel).

4.3. Identification and Characterization of Differentially Methylated Regions and Genes

First, to better understand the context of cytosine methylation, metrics that distinguish cytosines directly linked to guanines (CpG), indirectly linked to guanines (CHG), or not linked to guanines (CHH) were obtained with the MethylDackel v0.5.0 tool, in which H is the IUPAC ambiguity code for any nucleotide other than G. In addition, coverage metrics, correlation of methylation levels between samples, and chromosomal distribution of methylated cytosines were generated with the R package methylKit [55].

To identify differential DNA methylation patterns, the methylKit v1.24.0 package [55] was used for previous coverage filtering (>10x) and normalization of T/C ratios across samples. In addition, for each sample, 1000 bp tiles were formed every 500 bp (sliding tiles), which allowed to calculate the T/C ratios per delimited and overlapping genomic regions. Based on Fisher’s exact test, it was then possible to identify positions containing differentially methylated cytosines (DMCs) between groups, as well as genomic regions containing average proportions of DMCs between groups (DMRs). For more robust identification of DMCs and DMRs, in addition to the p-value adjusted for multiple comparisons (qvalue < 0.05), the minimum difference of methylated cytosines of 20% between groups was also used in both cases (DMCs and DMRs).

For description of the basic DMR annotation, the plotTargetAnnotation function was applied through the interface of the genomation [56] and methylKit [55] packages of R [57], in which DMRs can overlap with only promoter regions, intergenic regions, exons, and introns. Overlaps with CpG islands and shores are also observed.

For the subsequent gene expression and methylation combined analyses, which require a greater amount of genomic information in terms of DMR annotation, were performed a customized annotation script using NCBI RefSeq and CpG_island files obtained from the UCSC Genome Browser (https://genome.ucsc.edu/cgi-bin/hgTables) based on the URS-ACD1.2 assembly. In this case, in addition to the information obtained above, the position of introns and exons was discriminated, as well as the overlap between CpG islands and shores in gene and promoter regions. To definition of the promoter region was considered 2000 bp downstream and 500 bp upstream of the transcription start site (TSS).

4.4. Relative Gene Expression

To access to the expression level of the genes studied here, a transcriptome dataset of Longissimus thoracis muscle obtained for the same samples as used in this experiment was employed, which are available at https://drive.google.com/drive/folders/113q1QKulYBU22EC7ZwtbWk-2ORFIEzek?usp=share_link. The data obtained from the n x m matrix containing the genes symbol (n) and the number of samples (m) were normalized using counts per million (CPM) to obtain the relative expression between treatments (CW and EW), given as the logarithmic function of fold change (logFC). The significance values of the difference between mean treatment counts were obtained with a generalized model assuming a negative binomial distribution, followed by application of the likelihood ratio test assuming an F-distribution of pvalues. For this purpose, the dispersion parameter was estimated in a standard manner and TNM normalized data were obtained considering the relative size of the libraries. The edgeR v. 4.2.0 package [58] in R v. 4.3.3 was used for these differential expression analyses [57].

4.5. Relationship between Differentially Methylated Genes and Genes Expression

To investigate the relationship between differentially methylated genes (DMGs) and their expression levels, graphical representation and statistical associations between methylation percentage and gene expression level were performed. In general, the positions and relative methylation levels (hyper- or hypomethylated) of genomic regions to which these DMRs were mapped were analyzed in relation to the level of gene expression (over- or underexpressed). Next, the correlation between the gene expression level given as logFC of DMGs and the methylation percentage of the DMRs was analyzed using Pearson’s correlation coefficient. Finally, a trend line was fitted by simple linear regression in which the level of gene expression given as logFC was the response variable and the methylation level expressed as percentage was the explanatory variable. For this, only genes that were DMRs and DMGs was considerated.

4.6. Functional Enrichment of Differentially Methylated Genes

Enrichment analysis of gene ontology biological processes (GO_BP) and metabolic pathways (KEGG) was performed using the EnrichR software [59], considering both DMGs and differentially expressed genes (DEGs) obtained by differential gene expression analysis [44]. GO terms and KEGG pathways with a pvalue < 0.01 and qvalue < 0.1 were considered to be significantly enriched.

5. Conclusions

The data obtained in this study suggest that early weaning of Nellore cattle alters the methylation profile of some genes linked to metabolism. The ADCY3 gene is one of the main genes that were differentially methylated and differentially expressed. Functional enrichment analysis suggested that the cAMP signaling pathway may play an important role in the metabolic changes observed in early weaned cattle. Cyclic AMP signaling is involved in adipogenesis and lipogenesis processes and the DEG data found here demonstrated strong alteration in the expression of genes involved in de novo lipogenesis. The differential methylation of metabolism genes due to early weaning supplementation found in this study suggests that ADCY3 may be one of keys gene for understanding metabolic programming associated with fat deposition pathways in cattle and indicate that epigenetic mechanisms such as DNA methylation may be involved in the energetic metabolism adaptation caused by diet and environmental modification in Nellore calves early weaned.