The impact of different feeds on DNA methylation, glycolysis/gluconeogenesis signaling pathway, and gene expression of sheep muscle

Animal experiments and sampling

To reduce the effect of genetic background, four pairs of 3-month-old Chinese Sunit lamb (Ovis aries) female twins weighing 24 ± 2.3 kg were used in the pairing experiments. They were randomly divided into two groups and fed for 110 days either with HFLP or LFHP pellets. The twin lambs were selected from Siziwangqi, Jining, China. All the lambs were fed during the test period. Immunization, deworming, bathing, disease prevention and control were carried out as usual, and water was given ad libitum. The two groups were kept in a small corral. After a fasting period, the lambs were slaughtered in a local abattoir, and were euthanised through captive bolt stunning, followed by exsanguination. After the slaughter, two pieces of LD muscle between the 12th and 13th ribs were sampled from each individual and preserved in a nitrogen canister.

WGBS and identification of DMRs

The Illiumina HiSeqTM2500 platform was used to sequence libraries (Biomarker Technologies, Beijing, China). A base detection was done to convert the peak signal to sequence data, and the original reading was quality filtered to obtain a pure reading. The reading of the 3′ junction subsequence was first trimmed, and the values that had more than 10% unknown base (N) and low quality (more than 50% PHRED score ≤5 base) were discarded. Also, the Q30 and CG contents were estimated at the same time. The clean reads were aligned to the sheep reference genome (Oar_v4.0) (https://www.sheephapmap.org/), and the bisulfite mapping of methylation sites was performed using Bismark software. The sequencing depth and coverage were assessed after the duplicate reads were aligned with the same region of the genome. The bisulfite conversion rate is the proportion of methylated clean reads in the genome compared to the total number of clean reads. To confirm C-site methylation by screening conditions for ≥4× and false discovery rate (FDR) <0.05, the binomial distribution test for each C site was used (Fan et al., 2020).To identify DMR, the model of the estimated methylation level was referenced. All C sites with read coverage >10× were used for DMR analysis with MOABS (Sun et al., 2014). There were at least three differential methylation sites, and the minimum methylation level difference was 0.1 (CG type was 0.2) using Fisher’s exact test, P < 0.1.

Functional enrichment analysis

To annotate gene functions, DMGs were compared with the functional databases (GO, KEGG, and COG) by BLAST. The WALLSEius non-center hypergeometric distribution of the GOseq R package (Young et al., 2010) was used for DMGs GO enrichment analysis. In KEGG pathway analysis, the KOBAS software was used to determine the relevance of DMG enrichment (Mao et al., 2005). An enrichment analysis of the differential genes was performed with ClueGO (Bindea et al., 2009), which is a plugin of Cytoscape (Shannon et al., 2003). The GO annotation was used to classify the DMGs and DEGs into cellular components, molecular functions, and biological processes.

DNA methylation mapping and patterns

The entire DNA methylation analysis of the LD muscle of the twin lambs was carried out by WGBS with 10× genome coverage and >99% conversion rate. Approximate values of 293.17 M clean reads were uniquely mapped to produce the twin lamb HFLP methylome, and 290.10 M clean reads for the production of the twin lamb LFHP methylome. The quality of sequencing data is detailed (Table 1). Approximately a total of 3.0% of all genomic sites were methylated in every group (Table 1). Methylation in the twin lambs was detected to appear in three contexts of the sequence, namely, CG, CHG, and CHH (where H = A/C/T). However, their methylation levels were different. The entire genome-wide mC levels of 90.21% (CG), 2.12% (CHG), and 7.67% (CHH) within the HFLP specimen (Fig. 1A), and 88.56% (CG), 2.50% (CHG), and 8.93% (CHH) within the LFHP specimen (Fig. 1B). These contexts of sequence were almost the same in every group.

Characterization of DMRs and gene annotation

The DMR were determined within the two samples, and were annotated into the gene functional regions based on varied methylation contexts. A total of 1,945 CG DMRs, 127 CHH DMRs, and three CHG DMRs were detected. Most of these CG DMRs were detected in the distal intergenic region, with only 19 and 30 found in the 5′UTR and 3′UTR, respectively. CHG DMRs were all located in other introns, and CHH DMRs were mainly located in the distal intergenic regions. Of the 1,945 CG DMRs, a total of 1,022 were hypermethylated and 923 were hypomethylated. Out of the 127 CHH DMRs, a total of 105 were hypomethylated and 22 were hypermethylated. However, all the three CHG DMRs were hypomethylated. Except for those in the distal intergenic regions, the ratio of DMRs in introns was the greatest for all methylation types. A total of 1,471 CG DMGs, 106 CHH DMGs, and three CHG DMGs were annotated (Fig. 2).

Function enrichment analysis of DMGs

Analysis of the COG, GO and KEGG pathway databases was carried out, and 1,471 DMGs were found in the DMRs. The COG analyses showed that DMGs were enriched in general function prediction only, mostly for CG. The GO analysis for the CG context showed that DMGs were substantially enriched in the classification of negative regulation of transcription from RNA polymerase II promoter, cell junction, and calcium ion binding. The 20 most substantially uniquely enriched muscle development associated GO terms for DMGs between HFLP and LFHP LD samples are displayed. The KEGG analysis for the CG context also discovered that the DMGs were greatly enriched in pathways in cancer, axon guidance, MAPK, and Wnt signaling pathways. The results revealed that DMGs, which are affected by DNA methylation, can influence muscle development, glucose and fatty acid metabolism (Table S1, Fig. S1).

According to GO function classification, DMGs were divided into molecular function, cellular components, and biological processes (BP). The cAMP was enriched in BP by five genes (LOC101106702, GNAI1, CFTR, PIK3CG, and SOX9), and the cellular processes, single biological processes, metabolic processes, and biological regulation. Increasing LOC101106702 promoted positive regulation of BP cAMP-dependent protein kinase activity. GNAI1 can regulate cAMP-mediated signaling pathway. CFTR and PIK3CG are responsible for cell response to cAMP, and SOX9 can induce cAMP-mediated signaling pathway (Fig. 3A).

The KEGG pathway enrichment analysis of DMGs involved a total of 49 functions. It was found that genes with differential DNA methylation were enriched when human diseases, environmental information processing, organismal systems, genetic information processing, and cellular processes were included in the KEGG enrichment analysis. In the process of environmental information processing, the cAMP pathway was significantly enriched, including 21 genes (MAPK10, DRD5, DRD1, MLLT4, SOX9, LOC101106702, AKT3, RRAS2, HHIP, SSTR1, RELA, GRIN3A, FSHR, PPP1R12A, LOC101104054, ATP2B1, GRIN2B, GNAI1, PIK3CG, CFTR, and VAV3). The glycolysis/gluconeogenesis pathway was significantly enriched, including four genes (HKDC1, PDHA2, PGM2, and PDHA1). The FOXO pathway was significantly enriched, including 13 genes (MAPK10, HOMER1, FOXO1, AKT3, FOXG1, BCL6, LOC101103187, EGFR, CHUK, IRS4, GADD45A, PIK3CG, and S1PR1) (Fig. 3B).

Gene expression of LD muscle under two feeds

The Illumina Hi-Seq 2500 was used to sequence 8 LD muscle samples (four repetitions per set), and 22.5 to 28.5 million 125-bp clean paired-end reads were generated. The statistical power of this experimental design, calculated in RNASeqPower (Hart et al., 2013) is 0.9917946 (rnapower (depth = 23.68, cv = 0.4, effect = 4, n = 4, alpha = 0.05)). A reference transcriptome using obtained RNA-seq reads was reconstructed by the Trinity software. There were 347,335 transcripts in the reference transcriptome, with an average length of 2,166.29 bp. The edgeR package in R was used to analyze the DEGs between the twin lambs under different diets. With P = 0.0001 and a false discovery rate (FDR) of 0.05, 487 differentially expressed transcripts belonging to 368 DEGs were discovered, including 67 novel genes.

A KEGG pathway enrichment analysis on the DEGs in the muscles of twin lambs under different diets (Fig. S2) was performed, and the PDK4 gene was determined to play a vital role in the regulation of glucose and fatty acid metabolism. The key gene, PDK4 (Log CPM = 11.584, P = 0.000, and FDR = 0.016), was found to exhibit the greatest expression levels in LD muscle samples under the feed conditions.

Joint analysis of DMGs and DEGs in LD muscle

It was noticed that 1,471 DMGs were found between the two groups (LFHP and HFLP). A total of 11 overlapping genes (AKAP7, ANKRD17, FKBP5, FOXO1, MIER1, NFYC, PITX2, PODN, TOR1AIP1, TTN, and UTRN) were detected between DEGs and DMGs (Fig. 4). Six of the 11 genes were hypermethylated (AKAP7, FKBP5, MIER1, PITX2, PODN, and UTRN), while the other five (ANKRD17, FOXO1, NFYC, TOR1AIP1, and TTN) were hypermethylated in the LFHP group compared with the HFLP group. Among the 11 genes, seven genes (AKAP7, MIER1, NFYC, PITX2, PODN, TTN, and UTRN) were expressed at higher levels in the LFHP diet than in the HFLP diet. Also, four genes (ANKRD17, FKBP5, FOXO1, and TOR1AIP1) had a higher expression level in the HFLP diet than those in the LFHP diet. Expression levels in the HFLP diet group were significantly higher than in the LFHP diet group (P = 0.02). The expression levels of FKBP5 and FOXO1 genes were significantly higher in the HFLP diet than in the LFHP diet (P = 0.03, P = 0.01). Of these genes, FKBP5 (Log CPM = 8.101, P = 0.000, and FDR = 0.000), which plays a role in the intracellular trafficking of heterooligomeric forms of steroid hormone receptors maintaining the complex in the cytoplasm when unliganded, had high expression levels in the two diets. FOXO1 (Log CPM = 4.461, P = 0.000, and FDR = 0.025) is involved in the regulation of gluconeogenesis by regulation of transcription from the RNA polymerase II promoter.

cAMP pathway analysis

A KEGG pathway enrichment analysis was executed using ClueGO to determine the functions of identified DMGs. A total of 274 enriched pathways were detected, including the cAMP signaling pathway (Table S1). A total of 21 DMGs genes were detected in the cAMP signaling pathway. From the data, 12 genes (MAPK10, DRD1, MLLT4, LOC101106702, RRAS2, HHIP, RELA, PPP1R12A, LOC101104054, GNAI1, PIK3CG, and VAV3) were hypomethylated, eight of the genes (DRD5, SOX9, AKT3, GRIN3A, FSHR, ATP2B1, GRIN2B, and CFTR) were hypermethylated in the LFHP group compared with the HFLP group. These genes (MAPK10, DRD5, DRD1, MLLT4, SOX9, and AKT3) are mostly related to many functions, such as glucose metabolism, lipid metabolism, production of insulin and leptin signaling, and insulin secretion by regulating the ERK/P38 signaling pathway. The cAMP activity in twin lambs fed with the HFLP diet was higher than that in twin lambs fed with the LFHP diet. In the HFLP diet, nine genes had a higher expression than those in the LFHP diet. However, two genes in the LFHP diet had a higher expression than in the HFLP group (Fig. 5). The expression of the CFTR gene in HFLP was significantly higher than that in the LFHP diet (P = 0.1), and the expression level of RRAS2 was detected to have an extremely significant difference in the HFLP diet than in the LFHP diet (P = 0.01). Based on the cAMP signaling pathway activity, the RRAS2 was extremely increased by the HFLP diet. The RRAS2, which is one of the Ras families, is capable of controlling insulin and glucose metabolism. Finally, the CFTR, which was significantly increased by the HFLP diet, is responsible for the movement of ions such as chloride and phosphorite ions across the cell, and could play a role by impairing insulin secretion (Fig. 5).

Glycolysis/gluconeogenesis pathway

Further analysis of the functional classification of DMGs using ClueGO software revealed that the glycolysis/gluconeogenesis pathway and the FOXO pathway were significantly enriched, and the FOXO pathway was the key pathway regulating gluconeogenesis. Of the 17 genes in both pathways, eight genes were hypermethylated (HOMER1, FOXO1, AKT3, FOXG1, EGFR, CHUK, S1PR1, and PDHA1) and eight genes were hypomethylated (MAPK10, BCL6, IRS4, GADD45A, PIK3CG, HKDC1, PDHA2, and PGM2) in the LFHP group compared with the HFLP group. Eight genes were highly expressed in the HFLP group than in the LFHP group (HOMER1, FOXO1, AKT3, EGFR, GADD45A, PIK3CG, S1PR1, and PDHA1), while three genes were highly expressed in the LFHP group than in the HFLP group (BCL6, CHUK, and PGM2). The HFLP group had considerably higher expression of HOMER1 and FOXO1 than the LFHP group (P = 0.01). There is a significant correlation between the upregulated gene expression and hypomethylation of HOMER1 and FOXO1 gene in HFLP group (Fig. 6).