DNA Methylation May Affect Beef Tenderness Through Signal Transduction in Bos Indicus

doi:10.21203/rs.3.rs-1415533/v1

Download PDF

Research Article

DNA Methylation May Affect Beef Tenderness Through Signal Transduction in Bos Indicus

https://doi.org/10.21203/rs.3.rs-1415533/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Background

Beef tenderness is a complex trait of economic importance for the beef industry. Understanding the epigenetic mechanisms underlying this trait may help improve the accuracy of breeding programs. However, little is known about epigenetic effects on Bos taurus muscle and their implications in tenderness, and no studies have been conducted in Bos indicus.

Results

Comparing methylation profile of Bos indicus skeletal muscle with contrasting beef tenderness at 7 days after slaughter, we identified differentially methylated cytosines and regions associated with this trait. Interestingly, muscle that became tender beef had higher levels of hypermethylation compared to the tough group. Enrichment analysis of predicted target genes suggested that differences in methylation between tender and tough beef may affect signal transduction pathways, among which G protein signaling was a key pathway. In addition, different methylation levels were found associated with expression levels of GNAS, PDE4B, EPCAM and EBF3 genes. The differentially methylated elements correlated with EBF3 and GNAS genes overlapped CpG islands and regulatory elements. GNAS, a complex imprinted gene, has a key role on G protein signaling pathways. Moreover, both G protein signaling pathway and the EBF3 gene regulate muscle homeostasis, relaxation, and muscle cell-specificity.

Conclusion

We present differentially methylated loci that may be of interest to decipher the epigenetic mechanisms affecting tenderness. Supported by the previous knowledge about regulatory elements and gene function, the methylation data suggests EBF3 and GNAS as potential candidate genes and G protein signaling as potential candidate pathway associated with beef tenderness via methylation.

Nelore

RRBS

cattle

GNAS

EBF3

shear force

epigenome

muscle

Tenderness is one of the most appreciated beef quality traits by consumers (1), which makes this a desired trait by producers. Knowledge of the biological processes controlling this trait is crucial to improve beef quality through better management decisions and animal breeding. Beef tenderness is a complex trait with intra and inter-breed variability (2). It is well known that most Bos taurus breeds present a more tender beef product when compared to beef sourced from Bos indicus (3). Thus, the mechanisms that differ between these different subspecies impacting tenderness merit investigation. As Nelore is the predominant Bos indicus breed in Brazil, and Brazil is one of the biggest exporters of bovine beef in the world (4), knowledge of specific mechanisms affecting tenderness in Bos indicus may have great impact in the international beef market.

Omics studies exploring the biological mechanisms underlying tenderness in beef cattle have been developed for both Bos taurus (5) and Bos indicus (6–11). However, the knowledge regarding the role of epigenetics on tenderness is still limited, mainly for Bos indicus.

Epigenetics refers to modifications that occur in gene expression without altering the DNA sequence, i.e., modifications that alter the phenotype without modifying the genotype (12, 13). Methylation is one of the most common and more studied epigenetic marks in mammals and consists of the addition of a methyl group in a DNA cytosine (12, 14–16). Generally, there is an inverse correlation between DNA methylation (DNAm) and gene expression level (17). Higher levels of methylation in the promoter region are related to gene silencing (14), whereas promoters of active genes show lower intensity of methylation (15). The effect of CpGs methylation on gene expression is modulated by complex mechanisms and depends on its genomic location, being context-dependent (18). Studying human cells, Varley and colleagues (18) found that methylation of CpGs close to transcription start site (TSS) generally has repressing effect on gene expression, however, methylated CpGs in gene body and far from TSS may have either positive or negative correlation with expression, depending of its genome context. For instance, if the gene body methylated CpG is in CpG islands (CpGi), either positive or negative correlation was found while those not located in CpGi were positively correlated with gene expression.

Few studies in bovine muscle have investigated the methylation profile, as reviewed by Wang and Ibeagha-Awemu (19), and there is only one relating it to beef tenderness (20). Among these few studies, Huang et al. (21) compared the methylation profile and its relationship with mRNA and miRNA expression levels in fetuses and adult muscle of Chinese Qinchuan cattle. Lately, the methylome of muscle was studied in Angus (20) and in a crossed taurus vs. indicus Brangus-Angus population (22). The unique study analyzing muscle methylation and tenderness was recently published by Zhao et al. (20) that identified several differentially methylated regions between animals with higher and lower tenderness in the Bos taurus population. To our knowledge, there are no studies available investigating the relationship between methylation and tenderness in Bos indicus animals. Therefore, since Bos taurus and Bos indicus may present difference in their phenotype and transcriptomes, we also expect differences in methylation between them, thus turning methylation studies in Bos indicus a need for the comprehension of molecular basis guiding tenderness. The objective of this study was to investigate the association of DNA methylation profile in Nelore (Bos indicus) skeletal muscle by comparing contrasting animals for tenderness. In addition, we overlapped our results with previous QTL studies in the same population, in order to find methylation that are candidates to be affecting the trait. Further, we screened whether these regions are in known regulatory elements using the recently published database from the FAANG pilot project (23) and verified the relationship between methylation and gene expression in muscle.

2.1. Animals and Phenotypic Data.

The Tough group had higher EBVSF, with average of 0.64, while the group of animals presenting higher beef tenderness had lower EBVSF values (Tender group), with average of -0.82. Groups' sample composition, phenotypes and EBVSF values are presented in Additional File 2: Table S1. The averages of SF measurements and EBVSFs were different between the two groups with a t-Student p = 1.227e-06 and p = 1.483e-07, respectively (Additional file 1: Fig. S1).

2.2. RRBS DNA sequencing and Nelore muscle methylation profile.

The sequence depth ranged from 29 to 55 million reads per sample with an average of 20 million reads per sample aligned to the ARS-UCD1.2/bosTau9 genome and average mapping efficiency of 47.34%. A full summary of read mapping statistics is presented in Additional file 2: Table S2. Approximately, 275 million cytosines were analyzed, ranging from 203 to 353 million per sample. We observed that an average of 23.02% of all the analyzed cytosines (including non-CpGs) captured in the RRBS were cytosines followed by guanine, i.e. CpG, and 9.03% were methylated CpGs (Additional file 2: Table S3). From these CpGs, 635,469 were identified in the 12 samples and were used for the differential methylation analysis. Most of these CpGs overlapped intergenic regions (54%), followed by promoters (29%) and introns (11%), with exons overlapping the lowest number of CpGs analyzed (6%) (Additional file 1: Figure S2A). Around 60% of the cytosines analyzed were in CpG islands while those overlapping CpG shores represented only 15% (Additional file 1: Figure S2B). We tested 635,469 CpGs and 193,249 regions that were common in the 12 samples. We found a low variability in global methylation profile across samples, as shown by the remarkably similar correlations, which ranged from 0.94 to 0.96 (Additional file 1: Figure S3).

2.3. Most of the differentially methylated CpG and regions were hypermethylated in animals with higher tenderness.

Although we observed a low global variation in the methylation profile, we identified 123 DMCs and 42 DMRs showing methylation differences between the two groups of extreme GEBVs for SF at 7 days of aging (q-value < 0.05 and methylation difference > 25%). Interestingly, most of the DMCs and DMRs had higher levels of methylation in the tender group (n = 87 and n = 36, respectively) than in the tough (n = 37 and n = 6, respectively). All the DMCs and DMRs are described in Additional file 2: Table S4 and Additional file 2: Table S5, respectively.

2.4. Genome distribution of differentially methylated cytosines and regions.

Thirty-two DMRs (76.19%) contained at least one DMC. We found DMCs in almost all chromosomes, except for BTA27 and BTA28, while DMRs were found on 21 different chromosomes (Fig. 1). Regarding the functional location, DMCs had a higher frequency in intergenic regions (66.67%) followed by the gene body: intron (21.74%), exon (9.76%) and promoters (2.44%) (Additional file 1: Figure S4A). On the other hand, DMRs were more present in introns (45.25%), followed by intergenic regions (35.71%), promoters (11.9%) and exons (7.14%) (Additional file 1: Figure S4B). As Additional file 1: Figure S5 shows, some DMCs and DMRs were in more than one functional location. Thus, only one was considered to calculate the percentage, with preference in the following order promoter > exon > intron, as described by Methylkit version 1.8.1 (24). In addition, most of the DMCs and DMRs were found outside of CpG island and CpG shore regions (Additional file 1: Figures S4C and S4D, respectively).

2.5. Enrichment analysis of target genes suggests different methylation levels controlling genes acting on signal transduction.

Based on the annotations generated by the FAANG Consortium (23), enriched target genes for DMCs and DMRs were identified and are available in Additional file 2: Table S6. In total, 23 biological processes GO terms were enriched and clustered in 11 groups using the kappa score. Several signal transduction pathways were enriched (Fig. 2) in four different groups represented by positive regulation of cytokine-mediated signaling pathway (GO:0001961), regulation of Rho protein signal transduction (GO:0035023), phospholipase C-activating G protein-coupled receptor signaling pathway (GO:0007200) and adenylate cyclase-activating G protein-coupled receptor signaling pathway (GO:0007189). One group represented by inward rectifier potassium channel activity clustered five biological processes related to ions channel activity (GO:0005242). We also found three biological processes related to neurons grouped and represented by myelination (GO:0042552). Other five group had a unique biological process. Phagocytosis and inward rectifier potassium channel activity had around 18% of its associated genes enriched in our analysis, while the other pathways ranged between 4 and 9%. A complete description of ClueGo results, containing all gene ontologies, statistics, and genes clustered is available in Additional file 2: Table S7.

2.6. DMCs and DMRs overlapped QTLs for shear force.

To find DMCs and DMRs that may be also in QTLs related with tenderness, we looked for overlapping shear force QTLs identified in a previous Bos indicus study (6). Six DMCs overlapped QTLs associated with shear force at day 7 after slaughter, from which two also overlapped QTLs for SF at 24 hours after slaughter. All but one of these DMCs were hypermethylated in the tender muscle. Moreover, we found DMC8 and DMC83 overlapping QTLs for shear force from studies in other populations. Only one DMR (DMR18) overlapped a QTL, which was associated with shear force at day 7. Additional file 2: Tables S4 and S5 describe the QTLs overlapping DMCs and DMRs, respectively.

2.7. Manually curated list of candidate methylation regions and cytosines affecting tenderness.

The initial genome-wide DNAm screen, the results from enrichment analysis and the available information in the literature provided information for the development of a manually selected gene list (Additional file 2: Table S8). Thirty-two DMCs and nine DMRs were potentially acting as regulators of tenderness and differences in methylation level may result in phenotypic changes. Eight DMCs and one DMR overlapped QTLs for shear force and had 15 targets or closest genes with expression data. The enrichment analysis resulted in 27 genes that had expression data and were targets of 22 DMCs and 8 DMRs. Then considering the target genes or closest genes, we manually checked for function descriptions available in the literature related to muscle development or tenderness. This resulted in the addition of myosin VI (MYO6) and EBF transcription factor 3 (EBF3) genes.

2.8. Four genes had expression correlated with methylation.

We performed statistical association tests between the methylation percent of the selected DMCs/DMRs and the expression of its target genes. Eleven animals had paired RRBS methylation and RNA-Seq expression data. The RNA-Seq data showed an average of 9 billion reads with around 95% being uniquely mapped. RNA-Seq alignment results are described in Additional file 2: Table S9. We tested 65 combinations of DMC/DMRs and genes, which were 32 DMCs with 41 genes and 9 DMRs with 11 genes (Additional file 2: Table S8). DMC12, DMC 23, DMC89, DMC90, DMC98 and DMR40 showed a high correlation (p < 0.05) with one of their predicted target or closest genes as shown in Table 1 and Additional file 1: Figure S6. The four genes were epithelial cell adhesion molecule (EPCAM), phosphodiesterase 4B (PDE4B), GNAS complex locus (GNAS) and EBF3. Besides the DMC12 and DMC98, the other three DMCs and DMR40 were in CpG islands and regulatory elements (Additional file 2: Table S10). In addition, we tested the differential expression between extremes for these four genes and we found PDE4B with higher expression in Tough (p = 0.017) while GNAS was higher in Tender group (p = 0.030) (Additional file 1: Figure S7). EPCAM and EBF3 did not show differential expression (p > 0.05) between animals with divergent tenderness values.

Table 1

DMCs and DMRs that presented high correlation (Pearson; p-value < 0.05) between methylation level and gene expression of their predicted targets. The table includes whether the DMC or DMR overlaps a CpG island (CpGi), the regulatory element overlapped by the DMC or DMR in muscle and adipose tissues of *Bos taurus*, the target gene and the statistical association analysis results.
DM	CpGi	Regulatory element state			Target gene^§	r	p
DM	CpGi	Muscle	Adipose		Target gene^§	r	p
DMC12	other	ATAC island		medium enhancer with ATAC	EPCAM	-0.64	0.033
DMC98	other	Quiescent		Quiescent	PDE4B	-0.63	0.038
DMC23	CpGi	Bivalen/poised TSS		Strongly active promoter/transcript	GNAS	0.75	0.0085
DMC89	CpGi	strongly active promoters/transcripts		bivalent/poised TSS	EBF3	-0.72	0.013
DMC90	CpGi	strongly active promoters/transcripts		bivalent/poised TSS	EBF3	-0.75	0.0073
DMR40	CpGi	strongly active promoters/transcripts		bivalent/poised TSS	EBF3	-0.81	0.0023

2.8.1. GNAS

The GNAS gene was predicted to be the target of two DMCs and two DMRs, however, only the DMC23, which is located in a CpG island, showed correlations between methylation and expression (Table 1, Additional file 1: Figure S6A). The methylation percent in the DMC23 had a high positive correlation (0.75) with GNAS expression (Fig. 6A). To provide further evidence that DMC23 may act as a regulatory element, we checked its position for histone modification marks and ATAC-Seq using the data provided by Kern et al. (23). Figure 3 shows the enrichment of histone marks like H3K4me1, H3K27ac and H3K4me3, as well as an open chromatin state found by ATAC-seq in the DMC23 surrounding region (vertical black line in Fig. 3), suggesting that this region has features of an active promoter or enhancer. Indeed, in Bos taurus this region was classified as “strongly active promoter/TSS” in adipose tissue. However, in muscle, the presence of H3K27me3 characterizes an element repressing the expression which, together with the abovementioned epigenetic marks, may define this region as a “bivalent/poised TSS”. The importance of this region as a regulatory element was corroborated in other tissues, since all 15 tissues studied by Kern et al. (23) showed promoter features (Additional file 2: Table S10). Muscle was the only tissue which presented bivalent/poised TSS feature in this region, however seven tissues were classified as the adipose, a.i., as strongly active promoter/TSS, two as an active TSS while the remaining four tissues showed at gene features “promoter transcribed”.

2.8.2. PDE4B

Although DMC98 overlapped the gene JAK1 and it was not correlated to this gene expression but was negatively correlated (-0.63) with its predicted target PDE4B (Additional file 1: Figure S6B; Table 1). This gene was enriched in two pathways: adenylate cyclase-activating G protein-coupled receptor signaling pathway and voltage-gated cation channel activity. From the six selected DMCs and DMRs, this was the only one overlapping a region with no recognized regulatory element (Quiescent), showing very low signs of histone marks and ATAC in 14 tissues out of 15, according to Kern et al. (23) (Additional file 2: Table S10 and Additional file 1: Figure S8).

2.8.3. EBF3

EBF3 was the most correlated gene, having correlation with three elements (DMC89, DMC90 and DMR40) which are located in a CpG island. We found high negative pairwise correlations between EBF3 gene expression and either DMC89 (-0.72), DMC90 (-0.75) or DMR40 (-0.81) (Additional file 1: Figure S6C-E). Other two DMCs were tested with this gene, however no correlation was found. The Bos taurus data suggests that the region overlapped by the two DMCs and the DMR40 (vertical blue line in Fig. 4) has promoter features in Bos taurus (23) (Fig. 4, Additional file 2: Table S10). In muscle, this region was classified as strongly active promoters/transcripts which is characterized by the enrichment for two histones marks (H3K27ac and H3K4me3) and ATAC peaks. In adipose tissues, the region was part of a bivalent/poised TSS, characterized by the enrichment of all histones and ATAC marks, but mainly the H3K27me3, that has a repression function. Moreover, according to the Bos taurus data, EBF3 was expressed only in adipose and muscle tissue (Additional file 1: Figure S9), with higher levels in adipose (Fig. 4).

2.8.4. EPCAM

The DMC12, which was selected because it overlapped a QTL for shear force in the Nelore population, is in an intergenic region and more than 200,000 bp far away from its target EPCAM gene. We found a negative correlation (-0.64) between the methylation percent of DMC12 and the EPCAM expression (Table 1, Additional file 1: Figure S6E). The region overlapped by DMC12 also showed regulatory element features in several tissues of Bos taurus, including adipose and muscle (Additional file 1: Figure S10). Muscle had this region classified as an ATAC island and this state seems to be conserved in other five out 15 analyzed tissues (cecum, hypothalamus, rumen, abomasum, and spleen). In adipose tissue, this region was classified as medium enhancer based on ATAC peaks (Additional file 1: Figure S10).

3.1. A new layer of knowledge regarding tenderness: methylation.

The genetic mechanisms underlying beef tenderness have been extensively studied in the past years; however, little is known about the epigenetic marks affecting this important production trait. Epigenetic marks, such as DNAm, may affect the phenotype without altering the genome sequence, being an additional layer of information, which may improve comprehension of the biological processes underlying beef tenderness. Furthermore, epigenetic information complements the known QTLs and genes associated with the phenotype, helping to elucidate the intermediate mechanism between genotype and phenotype. To increase the epigenetic information related to tenderness, we screened CpGs and genomic regions looking for those showing differential methylation between muscle samples with extreme beef tenderness phenotypes. We found DMCs and DMRs overlapping QTLs associated with tenderness and overlapping regulatory elements, both previously identified (6, 23). Enrichment analysis showed that signal transduction and ions channel activity may be affected by different methylation status between tender and tough beef. Then, our results increase previous knowledge about mechanisms affecting beef tenderness in Bos indicus, since as far as we know, there are no previous studies in this subspecies.

3.2. DMCs and DMRs are more highly methylated in tender than in tougher beef.

The global methylation profile had low variation among animals. These findings may be, in part, due to the low genomic diversity found in the original population (25) and due to high conserved methylation profile in cattle (26). However, we found specific CpGs and regions with differences in methylation level among the two groups. Interestingly, muscles that will become higher tenderness beefs seem to be more methylated in DMCs and DMRs than the ones that will produce tough beefs. This finding corroborates Zhao et al. (20) who found genomic regions with high or medium levels of methylation had more hypermethylated DMRs in tender than in tough Angus muscle. It thus seems that both Bos indicus and Bos taurus show the same trend of hypermethylation in tender beef. This result may be complemented by recent reports that suggests CpG methylation may be a key mechanism in the determination of myofiber type and that slow-type myofibers, which are associated with higher tenderness, had more hypermethylated promoters affecting muscle-specific gene expressions than fast-type myofibers (27). Lu et al (28) suggested that increasing the number of slow-type myofiber in muscle of Chinese cattle (Luxi and Quinchuan) may improve beef tenderness. Unfortunately, in the present study the type of myofibers was not characterized.

3.3. Signal transduction and postmortem processing may be affected by methylation in muscle.

Identifying the genes that may be affected by the changes in methylation level may shed light on the mechanisms that link genotype to phenotype. In a recent study, Kern et al. (23) used several epigenetic marks and gene expression to find that most of the bovine regulatory elements do not target the nearest gene or even the gene that they overlap. These authors also reported that some regulatory elements targeted more than one gene, with an average of 1.7 genes per element in cattle. For this reason, we used their list of predicted target genes to identify those that may be affected by the DMCs and DMRs found here. Then we used the predicted target genes to perform functional enrichment analysis.

The biological processes enrichment suggests that signal transduction was affected by methylation, and this may be through the G protein signaling pathway. Four out of five groups of enriched pathways were directly related to G-protein-dependent signaling including adenylate cyclase-activating G protein-coupled receptor and phospholipase C-activating G protein-coupled receptor signaling pathways, regulation of Rho protein signal transduction and regulation of cytokine-mediated signaling pathway. Indeed, signal transduction through G protein-coupled is known to be a key mechanism for muscle homeostasis and growth (29). Signal transduction mediated by G-protein is responsible for converting the extracellular stimulus, reaching the membrane in downstream cascade of processes and resulting in reactions as muscle contraction (reviewed by 30). The hormones or other stimulus bind membrane G proteins-coupled receptor which in turn activates different subclasses of G proteins and its downstream pathways, like adenylate cyclase-activating G protein-coupled receptor and phospholipase C-activating G protein-coupled receptor signaling pathways (reviewed by 30), that were enriched in our analysis. Two genes grouped in the adenylate cyclase-activating pathways, GNAS and PDE4B, had their expression correlated with methylation level of the correspondent DMCs, adding further evidence to the hypothesis that G-protein signaling pathway is being affected by methylation and may affect beef tenderness. Rho protein signal transduction pathway can be activated by the G protein signaling and result in muscle contraction, this activation being mediated by Ca2+ (reviewed by 30–32). Takano et al. (33) found Rho family G protein acting on muscle differentiation in mice. G-protein signaling pathway was found acting on skeletal muscle atrophy through cytokine-mediated signaling pathways like tumor necrosis factor (29, 34). Both adenylate cyclase and phospholipase C-activating pathways regulate ion channel activities (reviewed by 30), which was also enriched in our analysis. Specifically, our analysis indicates that calcium and potassium channels are affected by methylation which corroborates Zhao et al. (20) findings in the muscle of Bos taurus. The direct participation of these channels in muscle development and function is well described (35). Tizioto et al. (2013a) found potassium and calcium channels enriched in a GWAS study for tenderness in a larger population of Nelore which included our animals. In addition, the enrichment of voltage-gated calcium channels and neuron-related biological processes in our analysis suggests that methylation may be affecting the muscle contraction signaling through the transfer of information from neuron to muscle in the neuromuscular junction (NMJ). The enriched calcium voltage-gated channel subunit alpha1 A (CACNA1A) gene encodes a subunit of the Cav2.1 channel which is located in the neuron terminal in the NMJ and has a key role in triggering muscular contraction in response to nerve impulse (36). The increase of Ca²⁺ ions influx in the cell by the Cav2.1 stimulates the release of neurotransmitters which will hit the receptor in the muscle cell and trigger downstream reactions to muscle contraction. Conversely, PDE4B that was also enriched in calcium voltage-gated channels has an opposite role by restricting the amount of Ca²⁺ entering the cell via these CaV1.2 channels in cardiomyocytes (37). Taking together these enriched pathways, we suggest that the transduction of signaling from extracellular stimulus to generation of second messengers are being affected and this may result in differences of muscle development and contraction.

Besides not related to signaling transduction, phagocytosis was one of the most associated genes enriched pathways and seems to be related with the postmortem process of the beef since muscle cells were found ingesting other cells and molecules from extracellular matrix by phagocytosis to degradation in the postmortem (38). Moreover, Tizioto et al. (6) found phagocytosis enriched for QTLs associated with beef tenderness on day 14 after slaughter.

3.4. DMC23 as a regulatory candidate site affecting tenderness through GNAS expression modulation

A key player of the G-signaling protein pathway is the G protein Gsα, encoded by the GNAS gene. After binding the extracellular stimulus, the membrane G protein-coupled receptors activate the Gsα which in turn stimulate the adenylyl cyclase pathway (39), enriched in our analysis. Using ATP, the enzyme adenylyl cyclase generates cyclic adenosine monophosphate (cAMP) which participates in several downstream pathways directly affecting skeletal muscle function including muscle contraction (40). Furthermore, bovine GNAS is a known imprinted gene with a complex locus that has two isoforms, according to the RefSeq database. The bovine NM_001271771.1 isoform, also known as NESP55, has an alternative splice site on exon 2 of the GNAS (41, 42) and has a maternal expression in bovine fetus (41–43) and adult (41) showing a paternal methylated promoter in B. taurus vs. B. indicus hybrid animals (42). The bovine NM_181021.3 isoform starts from exon 2 of GNAS locus and was shown to be paternally expressed in the conceptus of B. taurus vs. B. indicus hybrid animals (42). We found the DMC23 and DMR20 located in the intron 1 of GNAS locus, overlapping the NM_001271771.1 isoform, and located less than 1kb upstream from the TSS of the NM_181021.3 isoform. According to the pilot project from FAANG Consortium (23), these DMC23 and DMR20 overlap an active promoter region in muscle and adipose tissue and have GNAS as one of the target genes. It is important to keep in mind that our analysis is a pool of cells taken from skeletal muscle and this may include some adipose cells that are part of the intramuscular fat. Interestingly, the promoter was classified as poised in muscle, which is a promoter having epigenetic marks with both activation and repression effects allowing the gene to be activated or repressed as timely needed (44). This may explain the positive correlation between DMC23 methylation and the GNAS gene expression found in our samples, suggesting that the methylation may have an activation effect on the expression. A hypothesis is that the methylation may be disturbing the H3K27me3 binding, decreasing the heterochromatin formation which results in higher transcription of the gene. Moreover, the status of the region overlapped by the DMC23 as an active regulatory element seems to be conserved in other tissues (23). We present strong evidence that DMC23 overlaps a regulatory element that regulated the GNAS expression in muscle and adipose tissue. Therefore, DMC23 is a putative candidate site affecting tenderness due to the impact on gene expression of GNAS, which may affect beef tenderness through coupled receptor signaling pathways. Although little is known about the relationship of GNAS and tenderness, this gene was found downregulated in doubled-muscle bovine fetuses, known to result in higher shear force values (45) and enriched in the calcium signaling pathway (46). Polymorphisms and QTLs in this gene were associated with several production traits in bovine (43, 47). Moreover, abnormal methylation in the GNAS gene has been associated with disruption of several phenotypes, including diseases such as pseudohypoparathyroidism and obesity in humans (48, 49).

3.5. Differentially methylated CpG was negatively correlated with PDE4B expression.

Meanwhile GNAS gene product, Gsα, is known to activate the cAMP signaling, the phosphodiesterase enzymes (PDEs) seem to have opposite function by catalyzing the hydrolysis of cAMP decreasing its levels (50). In skeletal muscle of rats, the PDE4 was shown to be the major cAMP hydrolyzer in comparison to the other PDEs (51). This contradictory function of GNAS and PDE4 is in agreement with our results, since we found opposite expression levels of these two genes, and these expression levels seems to be regulated by DNA methylation, providing insights on the link of methylation and beef tenderness. We found higher transcript levels of PDE4B in tough beef while GNAS was more expressed in tender animals. This molecular pattern may suggest higher stimulation of G-protein-mediated cAMP signaling pathway in tender animals and this may reflect in the muscle functions including glycogenolysis, contractility, sarcoplasmic calcium dynamics (40) resulting in a different tenderness pattern. Interestingly, the expression of both genes was found correlated with methylation level of cytosines which were all hypermethylated in tender beef. Therefore, GNAS and PDE4B had opposite patterns of expression and methylation correlation between the beef with high and low tenderness suggesting that the signaling pathway that they partake, the cyclase-activating G protein-coupled, may have a role in the tenderness of Bos indicus beef and may be regulated by methylation.

3.6. Changes in methylation level ofEBF3 promoter as a potential candidate affecting beef tenderness.

The RRBS analysis identified EBF3 as a potential gene associated with tenderness through methylation. Several pieces of evidence among our findings and previous studies assure that DMC89, DMC90 and DMR40 may have a role in regulating EBF3 gene expression and tenderness. First, we found that these elements may regulate EBF3 gene expression due to the high negative correlation found with methylation. This is corroborated by Kern et al. (23) that found a correlation between EBF3 gene expression and H3K27ac signal overlapping the region where DMC89, DMC90 and DMR40 are located. Additionally, these elements are located in a promoter region and CpG island, which are known to regulate transcription and gene silencing (14). Second, EBF3 gene has specific functions in muscle. For instance, EBF3 gene encodes a transcription factor that acts together with MYOD gene in muscle relaxation and is a regulator of muscle cell-specific transcription (52). According to Kern et al. (23) data, EBF3 is a tissue specific gene being expressed only in adipose and muscle, compared to other four tissues. EBF3 is highly expressed in adipose of Bos taurus and with a lower expression in muscle (23). Intramuscular fat is known to affect skeletal muscle metabolism and beef tenderness (53, 54), therefore, EBF3 may be acting in muscle tenderness by both muscle and adipose pathways. Taken together, these findings make DMC89, DMC90, DMR40, and the EBF3 gene strong candidates for further studies of epigenetic events driving muscle development and tenderness.

3.7. DMC12 overlapped QTLs for beef tenderness.

Our methylation screen study resulted in regions and cytosines that are putative candidates underlying the phenotype beef tenderness in cattle, thus contributing to better understand the mechanisms regulating this trait. We highlight the DMC12 which lies within a QTL associated with tenderness in Nelore (6) and its methylation level was negatively correlated with the EPCAM gene expression. In Bos taurus, DMC12 overlapped regulatory elements in muscle and adipose (23). The biological function of the glycoprotein EPCAM, a transmembrane glycoprotein mediating Ca -independent homotypic cell–cell adhesion, has been explored in cancer studies (55) and it is known to participate in cell signaling (56), but little is known about its role in cattle.

3.8. Final considerations.

Here, we adopted the RRBS approach due to its advantages as base-level resolution and the accessible cost. However, this approach limited the analysis to part of the genome, which might contribute to the small variation in methylation profile found between the groups when compared to another study in Angus (20). Despite this limitation, we were able to provide new molecular insights to beef tenderness, which is a complex trait with a lot of its underlying mechanisms still undeciphered. We identified and suggested candidate regions and CpGs that may shed light on epigenetic mechanisms modulating beef tenderness by screening methylation on the Bos indicus muscle genome and taking advantage of previous information about regulatory elements and QTLs. Further studies are needed to answer complementary questions, among them, whether these methylation patterns are the cause or consequence of gene expression and beef tenderness, whether it was affected by the environment and how much it is inherited across generations.

We identified changes in DNAm level associated with tenderness, and functional analysis suggests that this link may be through modulating signal transduction, more specifically the G-protein-dependent signaling pathways, that have a specific function in muscle homeostasis. We also found strong evidence that differences in methylation level associated with tenderness are modulating GNAS and EBF3 gene expressions, thus making these genes potential candidates associated with the phenotype. The FAANG pilot project data in Bos taurus helped to identify potential functional role of the DMCs and DMR overlapping. Further studies of the DMCs and DMRs overlapping QTLs for SF may help complement the knowledge of the epigenetic mechanisms participating in the makeup of beef tenderness.

5.1. Experimental design and animals

This study aimed to identify epigenetics marks that may contribute to muscle development and resulting in differences in beef tenderness of Bos indicus. For this we investigated the association of DNA methylation profile in Nelore (Bos indicus) skeletal muscle by comparing contrasting animals for tenderness. The study population was composed of Nelore steers raised and kept on feedlots at Embrapa Southeast Livestock in São Carlos-SP, Brazil, as described by Tizioto et al. (57). The steers were sired by 33 registered Nelore bulls belonging to the main family lineages that make up the Nelore breed in Brazil. Steers were raised in feedlots with identical nutritional rations and handling conditions until slaughter at an average age of 25 months.

5.2. Phenotypes and sample selection.

Beef tenderness measurements were described by Tizioto et al. (6). Muscle samples were collected from the Longissimus thoracis skeletal muscle located between the 12th and 13th ribs were collected 24 h after slaughter. Steaks were used to measure the shear force (SF) using a TA XT2i (Stable Micro Systems Ltd., Surrey, United Kingdom) texture analyzer coupled to a Warner-Bratzler blade, at 1.016 mm thickness (58). Samples were kept in a 2°C chamber for aging and measurements of SF were obtained at 7 days after slaughter. Muscle samples were selected from the initial population by ranking the estimated breeding value (EBV) for SF (EBVSF) of each animal. The EBVs were estimated using standard best linear unbiased prediction (BLUP) procedures under an animal model performed using the mixed procedure of SAS (59–61) as reported by Gonçalves et al. (8). Twelve samples were selected and divided into two groups with extreme values and were named according to the level of tenderness: highest values of EBVSF (n = 6), corresponding to the lowest tenderness, were named as “Tough” while the lowest values of EBVSF (n = 6), corresponding to the highest tenderness, were named as “Tender”. Only animals born from different sires were considered for each group to avoid the confounding between phenotype and infinitesimal effect of sire.

5.3. DNA and RNA extraction and quality evaluation.

Longissimus thoracis skeletal muscle samples were collected at slaughter, stored in liquid nitrogen, and kept in a freezer at -80°C until processing. DNA extraction was carried out using the DNeasy® Blood & Tissue kit (Qiagen Inc., Germantown, MD, USA), according to the manufacturer's protocol. DNA concentration was measured using Qubit dsDNA High Sensitivity Assay (Thermo Fisher Scientific Waltham, MS, EUA) and the quality measured by Fragment Analyzer™ and the DNF-487 Standard Sensitivity or the DNF-488 High Sensitivity genomic DNA Analysis Kit (Advanced Analytical).

5.4. Reduced representation bisulfite sequencing (RRBS) and processing.

The RRBS experiments were performed by Diagenode® (Seraing, Belgium). Each sample library was prepared from 100 ng of genomic DNA using the Premium RRBS kit (Diagenode), according to the manufacturer's protocol. Bisulfite conversion efficiency independent of CpG context was assessed by adding methylated and unmethylated spike-in controls (0.1% concentration). Briefly, the protocol consists in the digestion of DNA by MspI enzyme followed by the fragment end repair, and addition of adaptors. Next, samples were quantified by qPCR and the Ct values were used to pool samples by similarity. Bisulfite conversions were performed using the Premium RRBS kit (Diagenode) followed by library enrichment by PCR, based on the manufacturer’s protocol. Adequate fragment size distributions were confirmed by Bioanalyzer High Sensitivity DNA chips (Agilent). Libraries were sequenced on Illumina HiSeq 3000 using single-end 50 bp reads. Sequencing read quality control was performed using FastQC version 0.10.1(https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and adaptors were removed by Trim Galore! version 0.4.1. (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). Bismark version 0.23.0 (62) was used to align the reads to the Bos taurus reference genome UCD1.2/bosTau9 and to identify methylated cytosines. The reference genome UCD1.2/bosTau9 fasta file was downloaded from UCSC Genome Browser in February/2021(http://hgdownload.soe.ucsc.edu/goldenPath/bosTau9/bigZips/bosTau9.fa.gz).

5.5. RNA sequencing and processing.

Eleven animals, out of the 12 used to analyze for DNA methylation, also had RNA sequencing data. These 11 RNA samples were a subset of the 200 reported by Diniz et al. (63). The methods for sample collection, total RNA extraction, library preparation, sequencing and read processing are described by Diniz et al. (63). In summary, muscle samples were collected immediately after slaughter, snap frozen in liquid nitrogen and kept at -80°C until RNA extraction. From 100 mg of frozen tissue, total RNA was isolated using Trizol® standard protocol (Life Technologies, Carlsbad, CA, United States) and the resulting mRNA was evaluated in the Bioanalyzer 2100® (Agilent, Santa Clara, CA, United States). The cDNA libraries were generated using Illumina TruSeq® RNA Sample Preparation Kit v2 (San Diego, CA, United States), followed by purification and validation using Agilent 2100 Bioanalyzer (Santa Clara, CA, United States). Library preparation and sequencing were conducted by ESALQ Genomics Center (Piracicaba, SP Brazil). Paired-end (PE) sequencing was executed on Illumina Hiseq 2500® (San Diego, CA, United States) platform following standard protocols. Read processing protocol consisted in filtration using Seqyclean package version 1.4.13 (64), and quality control using FastQC version 0.11.2 (65) and MultiQC version 1.4 (66).

5.6. Differential DNA methylation analysis between extremes.

The comparison between the dataset from two extreme groups for EBVSF phenotypes was executed using methylKit version 1.8.1 (24) in R version 3.5.1 (67). First, all cytosines were filtered for minimum coverage of 10 reads and maximum percentile of 99.9%. The methylKit function that normalize the coverage was applied using default settings. A Chisq-test was used to define the cytosines or genomic regions that were differentially methylated (DMCs or DMRs, respectively) between the two groups of animals. A correction for overdispersion was applied, as suggested by methylKit user guide. The DMRs window size was 100 bp and both DMCs and DMRs were filtered for at least 25% of difference of methylation between groups and q-value ≤ 0.05.

5.7. Gene annotation and Functional enrichment analysis.

Gene annotation was carried out using the package Genomation version 3.8 (68) through Methylkit version 1.8.1 (24) in R. The NCBI RefSeq annotation bed file used was obtained from UCSC Genome Browser. The webtool Faangmine (http://128.206.116.18:8080/faangmine/) was also used to find genes overlapping the DMCs and DMRs. Promoters were considered as 1,000 bp upstream and downstream from the gene’s TSS, following the default of Methylkit. CpG island and shore regions’ information were obtained from UCSC Genome Browser (69). To retrieve information whether DMCs and DMRs were overlapping regulatory elements, we used the FAANG pilot project (23). These data provided information about the regulatory element status of the region based on enrichment of histone marks (H3K27me3, H3K4me1, H3K27ac, H3K4me3), ATAC-Seq and CTCF data in 15 tissues of two male Bos taurus. The method to classify the genomic region regarding the regulatory element state was described by Pan et al. (70). In summary, the epigenetic marks data described above were used to train a chromatin state prediction model in ChromHMM69 (v.1.20) and 15 states were chosen. We also used these data to find the target genes of the DM regions and the parameters used to download this information from UCSC were Assembly Apr. 2018 (ARS-UCD1.2/bosTau9), group UC Davis FAANG Pilot Project and track Cattle Regulator-Gene Interactions to retrieve the target genes. The UCSC Table Browser tool (69) was accessed in May of 2021. The list of predicted target genes was submitted to functional enrichment analysis using ClueGO v. 2.5.5. software and CluePedia v. 1.5.5 (71). Significantly enriched gene ontology (GO) and biological process (BP) terms were declared with a right-sided hypergeometric test at a Benjamini and Hochberg adjusted p < 0.05 (72). We searched for DMCs and DMRs whose position overlapped QTLs previously associated with tenderness by GWAS (6) in the original population. QTL positions obtained from Tizioto et al. (6) using the UMD 3.1.1 reference genome were converted for the ARS-UCD1.2 version using Lift Genome Annotations UCSC tool (https://genome.ucsc.edu/cgi-bin/hgLiftOver). We also used the webtool Faangmine (http://128.206.116.18:8080/faangmine/) to find tenderness QTLs from other populations in the CattleQTLdb (73), and overlapped the positions with the DMCs and DMRs genomic positions.

5.8. Development of a manually curated list of candidate methylation regions and cytosines affecting tenderness.

We used the initial genome-wide DNAm screen together with information available in the literature to create a manually selected list of regions or CpGs that are potentially acting as regulators of tenderness and whose differences in methylation level may lead to phenotypic changes. First, we selected all the DMCs and DMRs that overlapped QTLs for shear force. We included the target genes and the closest gene of these DMCs and DMR. Then, all target genes that were clustered in any of the enriched gene ontology biological processes were included to the list. Finally, we complemented the list with DMCs/DMRs whose target genes had function related with muscle development or tenderness, based on literature information. The manually curated list of candidate DMCs/DMRs and genes were submitted to statistical association tests between methylation and expression of the target genes

5.9. Correlation between methylation percentage and gene expression.

To investigate the relationship between DMCs or DMRs and its respective target genes, we carried out statistical associations between methylation level and gene expression. The methylation percentage of DMCs and DMRs were obtained from methylKit version 1.8.1 (24) in R using the function percMethylation(). To access the gene expression of target genes, the pre-processed reads obtained from RNA-Seq were used. High-quality reads were filtered to remove extremely low expressed genes among the 11 samples, retaining genes with at least seven non-zero samples. Count normalization was executed using the standard method from DESeq2 version 1.28.1 (74). The correlation between the DMC or DMR percentage of methylation and read counts was analyzed using Pearson correlation in R by cor.test() function. For those genes whose expression were correlated to methylation, we also tested whether the expression was different between the groups of extreme phenotypes. For this we used Wilcoxon test applying the R function wilcox_test() with paired = False.

cAMP, Cyclic adenosine monophosphate

CpG, cytosine followed by guanine in a linear sequence

CpGi, CpG island

DMC, differentially methylated cytosines

DMR, differentially methylated cytosines

DNAm, DNA methylation

EBF3, EBF transcription factor 3

EBVSF, estimated breeding value for shear force

EPCAM, epithelial cell adhesion molecule

GNAS, GNAS complex locus

GO, gene ontology

NMJ, neuromuscular junction

PDE4B, phosphodiesterase 4B

QTL, quantitative trait loci

RRBS, Reduced Representation Bisulfite Sequencing

SF, shear force

Ethical approval

All procedures involving steers were approved by the Institutional Animal Care and Use Committee Guidelines from Brazilian Agricultural Research Corporation – EMBRAPA (process number: Macroprograma 1, 01/2005) and sanctioned by the president Dr. Rui Machado to ensure compliance with international guidelines for animal welfare.

Consent for publication

Not applicable

Competing interests

The authors declare that they have no competing interests

Availability of data and material

The dataset generated and analyzed during the current study regarding the methylation experiments is available in the NCBI’s Gene Expression Omnibus (GEO) under accession numbers GSE190966. The data will remain private until the publication of the manuscript. The data is available to review at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE190966 using the secure token kverccyshduxzkx.

The RNA-Seq dataset is available in the European Nucleotide Archive (ENA) repository (EMBL-EBI), under accession PRJEB13188, PRJEB10898, and PRJEB19421 (https://www.ebi.ac.uk/ena/submit/sra/).

The dataset regarding the histone marks, ATAC-Seq and regulatory elements are available at http://farm.cse.ucdavis.edu/~ckern/Nature_Communications_2020/.

Funding

Research presented in this publication was funded by National Council for Scientific and

Technological Development, CNPq (grant number: 456191/2014-3) and Foundation to Support the

Research of São Paulo State, FAPESP (grant numbers: 2012/23638-8 and 2019/04089-2). MS was granted with CNPq [167672/2017-7]. GM, LR and LC are recipients of a CNPq research productivity fellowship.

Author’s contributions

M.S. contributed to proposal elaboration for financial support, to animal and laboratory experiments, carried out bioinformatics analysis, and wrote the manuscript. S.N. contributed to proposal elaboration for financial support, designed the study, supervised the study, read, discussed, and approved the manuscript. M.R., contributed to proposal elaboration for financial support, to animal and laboratory experiments, carried out bioinformatics analysis and reviewed the manuscript. Z.P. and H.Z. provided the data and performed the analysis to classify regulatory element regions. W.D., J.A., J.B., and P.S.N.O. contributed to animal and laboratory experiments, part of the data analysis, and reviewed the manuscript. A.Z. was responsible for data storage and management, read, discussed, and approved the manuscript. G.M. was responsible for the phenotypic analysis, read, discussed, and approved the manuscript. L.C. participated on experimental design, read, discussed, and approved the manuscript. J.K. coordinated the bioinformatics analysis, read, discussed, and approved the manuscript. L.R. performed experimental design, collected samples, supervised and coordinated the study, read, discussed, and approved the manuscript.

Acknowledgements

We would like to thank the Bioinformatic Multi-user Laboratory, Embrapa Informática Agropecuária (Campinas – Brazil) for the infrastructure for Bioinformatic analysis. We also thank the support of Leandro Carrijo Cintra from Embrapa Informática Agropecuária (Campinas – Brazil).

Verbeke W, Van Wezemael L, de Barcellos MD, Kugler JO, Hocquette JF, Ueland O, et al. European beef consumers’ interest in a beef eating-quality guarantee Insights from a qualitative study in four EU countries. Appetite [Internet]. 2009/12/08. 2010;54(2):289–96. Available from: http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Retrieve&db=PubMed&dopt=Citation&list_uids=19961887
Judge MM, Conroy S, Hegarty PJ, Cromie AR, Fanning R, Kelly D, et al. Eating quality of the longissimus thoracis muscle in beef cattle – Contributing factors to the underlying variability and associations with performance traits. Meat Sci. 2021;
O’Connor SF, Tatum JD, Wulf DM, Green RD, Smith GC. Genetic effects on beef tenderness in Bos indicus composite and Bos taurus cattle. J Anim Sci [Internet]. 1997/07/01. 1997;75(7):1822–30. Available from: http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Retrieve&db=PubMed&dopt=Citation&list_uids=9222838
USDA USD of A. Livestock and Poultry: World Markets and Trade [Internet]. 2021 [cited 2021 Aug 23]. Available from: https://apps.fas.usda.gov/psdonline/circulars/livestock_poultry.pdf
Leal-Gutiérrez JD, Elzo MA, Johnson DD, Hamblen H, Mateescu RG. Genome wide association and gene enrichment analysis reveal membrane anchoring and structural proteins associated with meat quality in beef. BMC Genomics. 2019;
Tizioto PC, Decker JE, Taylor JF, Schnabel RD, Mudadu MA, Silva FL, et al. Genome scan for meat quality traits in Nelore beef cattle. Physiol Genomics [Internet]. 2013;45(21):1012–20. Available from: http://physiolgenomics.physiology.org/cgi/doi/10.1152/physiolgenomics.00066.2013
de Souza MM, Zerlotini A, Rocha MIP, Bruscadin JJ, Diniz WJ da S, Cardoso TF, et al. Allele-specific expression is widespread in Bos indicus muscle and affects meat quality candidate genes. Sci Rep. 2020;
Gonçalves TM, de Almeida Regitano LC, Koltes JE, Cesar ASM, da Silva Andrade SC, Mourão GB, et al. Gene Co-expression Analysis Indicates Potential Pathways and Regulators of Beef Tenderness in Nellore Cattle. Front Genet [Internet]. 2018 Oct 5 [cited 2018 Nov 27];9:441. Available from: https://www.frontiersin.org/article/10.3389/fgene.2018.00441/full
Kappeler BIG, Regitano LCA, Poleti MD, Cesar ASM, Moreira GCM, Gasparin G, et al. MiRNAs differentially expressed in skeletal muscle of animals with divergent estimated breeding values for beef tenderness. BMC Mol Biol. 2019;
Silva VH da, Regitano LC de A, Geistlinger L, Pértille F, Giachetto PF, Brassaloti RA, et al. Genome-Wide Detection of CNVs and Their Association with Meat Tenderness in Nelore Cattle. Davoli R, editor. PLoS One. 2016 Jun;11(6):e0157711.
Tizioto PC, Gasparin G, Souza MM, Mudadu MA, Coutinho LL, Mourão GB, et al. Identification of KCNJ11 as a functional candidate gene for bovine meat tenderness. Physiol Genomics [Internet]. 2013 Dec 15 [cited 2015 Jul 3];45(24):1215–21. Available from: http://www.ncbi.nlm.nih.gov/pubmed/24151244
Bradbury J. Human epigenome project–up and running. PLoS Biol [Internet]. 2003/12/24. 2003;1(3):E82. Available from: http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Retrieve&db=PubMed&dopt=Citation&list_uids=14691553
Lieb JD, Beck S, Bulyk ML, Farnham P, Hattori N, Henikoff S, et al. Applying whole-genome studies of epigenetic regulation to study human disease. Cytogenet Genome Res [Internet]. 2006/05/24. 2006;114(1):1–15. Available from: http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Retrieve&db=PubMed&dopt=Citation&list_uids=16717444
Enright BP, Jeong BS, Yang X, Tian XC. Epigenetic characteristics of bovine donor cells for nuclear transfer: levels of histone acetylation. Biol Reprod [Internet]. 2003/06/13. 2003;69(5):1525–30. Available from: http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Retrieve&db=PubMed&dopt=Citation&list_uids=12801976
Tost J. DNA methylation: an introduction to the biology and the disease-associated changes of a promising biomarker. Mol Biotechnol [Internet]. 2009/10/21. 2010;44(1):71–81. Available from: http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Retrieve&db=PubMed&dopt=Citation&list_uids=19842073
Attar N. The allure of the epigenome. Genome Biol [Internet]. 2012/10/25. 2012;13(10):419. Available from: http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Retrieve&db=PubMed&dopt=Citation&list_uids=23092546
Lokk K, Modhukur V, Rajashekar B, Märtens K, Mägi R, Kolde R, et al. DNA methylome profiling of human tissues identifies global and tissue-specific methylation patterns. Genome Biol [Internet]. 2014;15(4):r54. Available from: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=4053947&tool=pmcentrez&rendertype=abstract
Varley KE, Gertz J, Bowling KM, Parker SL, Reddy TE, Pauli-Behn F, et al. Dynamic DNA methylation across diverse human cell lines and tissues. Genome Res. 2013;
Wang M, Ibeagha-Awemu EM. Impacts of Epigenetic Processes on the Health and Productivity of Livestock. Frontiers in Genetics. 2021.
Zhao C, Ji G, Carrillo JA, Li Y, Tian F, Baldwin RL, et al. The Profiling of DNA Methylation and Its Regulation on Divergent Tenderness in Angus Beef Cattle. Front Genet. 2020;
Huang YZ, Sun JJ, Zhang LZ, Li CJ, Womack JE, Li ZJ, et al. Genome-wide DNA methylation profiles and their relationships with mRNA and the microRNA transcriptome in bovine muscle tissue (Bos taurine). Sci Rep. 2014;3.
Liu L, Amorín R, Moriel P, DiLorenzo N, Lancaster PA, Peñagaricano F. Maternal methionine supplementation during gestation alters alternative splicing and DNA methylation in bovine skeletal muscle. BMC Genomics 2021 221 [Internet]. 2021 Oct 30 [cited 2021 Nov 11];22(1):1–11. Available from: https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-021-08065-4
Kern C, Wang Y, Xu X, Pan Z, Halstead M, Chanthavixay G, et al. Functional annotations of three domestic animal genomes provide vital resources for comparative and agricultural research. Nat Commun. 2021 Dec 1;12(1).
Akalin A, Kormaksson M, Li S, Garrett-Bakelman FE, Figueroa ME, Melnick A, et al. MethylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles. Genome Biol. 2012;13(10).
Mudadu MA, Porto-Neto LR, Mokry FB, Tizioto PC, Oliveira PSN, Tullio RR, et al. Genomic structure and marker-derived gene networks for growth and meat quality traits of Brazilian Nelore beef cattle. BMC Genomics [Internet]. 2016 Dec 15 [cited 2016 Apr 3];17(1):235. Available from: http://www.biomedcentral.com/1471-2164/17/235
McKay S, Betancourt F, Bhattarai S, Buttolph T, White S, Lachance H, et al. 115 Profiling Conservation of DNA Methylation in Cattle. J Anim Sci. 2018;
Oe M, Ojima K, Muroya S. Difference in potential DNA methylation impact on gene expression between fast- and slow-type myofibers. Physiol Genomics. 2021;
Lu X, Yang Y, Zhang Y, Mao Y, Liang R, Zhu L, et al. The relationship between myofiber characteristics and meat quality of Chinese Qinchuan and Luxi cattle. Anim Biosci. 2021;
Guttridge DC. Making muscles grow by G protein-coupled receptor signaling. Science Signaling. 2011.
Kuo IY, Ehrlich BE. Signaling in muscle contraction. Cold Spring Harb Perspect Biol. 2015;
Sah VP, Seasholtz TM, Sagi SA, Brown JH. The role of Rho in G protein-coupled receptor signal transduction. Annual Review of Pharmacology and Toxicology. 2000.
Puetz S, Lubomirov LT, Pfitzer G. Regulation of smooth muscle contraction by small GTPases. Physiology. 2009.
Takano H, Komuro I, Oka T, Shiojima I, Hiroi Y, Mizuno T, et al. The Rho Family G Proteins Play a Critical Role in Muscle Differentiation. Mol Cell Biol. 1998;
Minetti GC, Feige JN, Rosenstiel A, Bombard F, Meier V, Werner A, et al. Gαi2 signaling promotes skeletal muscle hypertrophy, myoblast differentiation, and muscle regeneration. Sci Signal. 2011;
Stefani E, Chiarandini DJ. Ionic channels in skeletal muscle. Annual review of physiology. 1982.
Taverna E, Saba E, Rowe J, Francolini M, Clementi F, Rosa P. Role of Lipid Microdomains in P/Q-type Calcium Channel (Cav2.1) Clustering and Function in Presynaptic Membranes *. J Biol Chem [Internet]. 2004 Feb 13 [cited 2021 Aug 2];279(7):5127–34. Available from: http://www.jbc.org/article/S0021925820748362/fulltext
Leroy J, Richter W, Mika D, Castro LRV, Abi-Gerges A, Xie M, et al. Phosphodiesterase 4B in the cardiac L-type Ca2 + channel complex regulates Ca2 + current and protects against ventricular arrhythmias in mice. J Clin Invest [Internet]. 2011 [cited 2021 Aug 2];121(7):2651. Available from: /pmc/articles/PMC3223827/
Ouali A, Gagaoua M, Boudida Y, Becila S, Boudjellal A, Herrera-Mendez CH, et al. Biomarkers of meat tenderness: Present knowledge and perspectives in regards to our current understanding of the mechanisms involved. Meat Science. 2013.
Sunahara RK, Dessauer CW, Whisnant RE, Kleuss C, Gilman AG. Interaction of G(sα) with the cytosolic domains of mammalian adenylyl cyclase. J Biol Chem. 1997;
Berdeaux R, Stewart R. cAMP signaling in skeletal muscle adaptation: Hypertrophy, metabolism, and regeneration. American Journal of Physiology - Endocrinology and Metabolism. 2012.
Khatib H. Imprinting of Nesp55 gene in cattle. Mamm Genome. 2004;
Chen Z, Hagen DE, Wang J, Elsik CG, Ji T, Siqueira LG, et al. Global assessment of imprinted gene expression in the bovine conceptus by next generation sequencing. Epigenetics. 2016;
Sikora KM, Magee D a, Berkowicz EW, Berry DP, Howard DJ, Mullen MP, et al. DNA sequence polymorphisms within the bovine guanine nucleotide-binding protein Gs subunit alpha (Gsα)-encoding (GNAS) genomic imprinting domain are associated with performance traits. BMC Genet [Internet]. 2011;12(1):4. Available from: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=3025900&tool=pmcentrez&rendertype=abstract
Voigt P, Tee WW, Reinberg D. A double take on bivalent promoters. Genes and Development. 2013.
Uytterhaegen L, Claeys E, Demeyer D, Lippens M, Fiems LO, Boucqué CY, et al. Effects of double-muscling on carcass quality, beef tenderness and myofibrillar protein degradation in Belgian Blue White bulls. Meat Sci. 1994;
Cassar-Malek I, Passelaigue F, Bernard C, Léger J, Hocquette JF. Target genes of myostatin loss-of-function in muscles of late bovine fetuses. BMC Genomics. 2007;
Imumorin IG, Kim EH, Lee YM, De Koning DJ, van Arendonk JA, De Donato M, et al. Genome scan for parent-of-origin QTL effects on bovine growth and carcass traits. Front Genet. 2011;
Weinstein LS, Xie T, Qasem A, Wang J, Chen M. The role of GNAS and other imprinted genes in the development of obesity. International Journal of Obesity. 2010.
Rochtus A, Martin-Trujillo A, Izzi B, Elli F, Garin I, Linglart A, et al. Genome-wide DNA methylation analysis of pseudohypoparathyroidism patients with GNAS imprinting defects. Clin Epigenetics. 2016;
Moorthy BS, Gao Y, Anand GS. Phosphodiesterases catalyze hydrolysis of cAMP-bound to regulatory subunit of protein kinase A and mediate signal termination. Mol Cell Proteomics. 2011;
Bloom TJ. Cyclic nucleotide phosphodiesterase isozymes expressed in mouse skeletal muscle. Can J Physiol Pharmacol. 2002;
Jin S, Kim J, Willert T, Klein-Rodewald T, Garcia-Dominguez M, Mosqueira M, et al. Ebf factors and MyoD cooperate to regulate muscle relaxation via Atp2a1. Nat Commun. 2014;
Wood JD, Enser M, Fisher A V., Nute GR, Sheard PR, Richardson RI, et al. Fat deposition, fatty acid composition and meat quality: A review. Meat Science. 2008.
Scollan N, Hocquette JF, Nuernberg K, Dannenberger D, Richardson I, Moloney A. Innovations in beef production systems that enhance the nutritional and health value of beef lipids and their relationship with meat quality. Meat Science. 2006.
Keller L, Werner S, Pantel K. Biology and clinical relevance of EpCAM. Cell Stress. 2019.
Maetzel D, Denzel S, Mack B, Canis M, Went P, Benk M, et al. Nuclear signalling by tumour-associated antigen EpCAM. Nat Cell Biol. 2009;
Tizioto PC, Taylor JF, Decker JE, Gromboni CF, Mudadu MA, Schnabel RD, et al. Detection of quantitative trait loci for mineral content of Nelore longissimus dorsi muscle. Genet Sel Evol [Internet]. 2015;47(1):1–9. Available from: http://www.gsejournal.org/content/47/1/15
Wheeler TL, Shackelford SD, Johnson LP, Miller MF, Miller RK, Koohmaraie M. A Comparison of Warner-Bratzler Shear Force Assessment Within and among Institutions. J Anim Sci. 1997;
Henderson CR. SIRE EVALUATION AND GENETIC TRENDS Introduction. In: ASAS and ADSA, editor. Proc Animal Breeding and Genetics Symposium in Honor of Dr Jay L Lush [Internet]. Savoy, IL; 1972 [cited 2018 Dec 3]. p. 10–41. Available from: https://academic.oup.com/jas/article-abstract/1973/Symposium/10/4697208
Mrode RA, Thompson R. Linear models for the prediction of animal breeding values: Second Edition. Linear Models For the Prediction of Animal Breeding Values: Second Edition. 2005.
Sanford DA. Genetic Analysis of Complex Traits Using SAS. Crop Sci. 2005;
Krueger F, Andrews SR. Bismark: A flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics. 2011;27(11):1571–2.
Diniz WJS, Mazzoni G, Coutinho LL, Banerjee P, Geistlinger L, Cesar ASM, et al. Detection of co-expressed pathway modules associated with mineral concentration and meat quality in nelore cattle. Front Genet. 2019;
Zhbannikov IY, Hunter SS, Foster JA, Settles ML. Seqyclean: A pipeline for high-throughput sequence data preprocessing. In: ACM-BCB 2017 - Proceedings of the 8th ACM International Conference on Bioinformatics, Computational Biology, and Health Informatics. 2017.
Andrews S, Babraham Bioinformatics. FastQC: A quality control tool for high throughput sequence data [Internet]. Manual. 2010. p. http://www.bioinformatics.babraham.ac.uk/projects/. Available from: http://scholar.google.com/scholar?hl=en&btnG=Search&q=intitle:FastQC+a+quality+control+tool+for+high+throughput+sequence+data.#0
Ewels P, Magnusson M, Lundin S, Käller M. MultiQC: Summarize analysis results for multiple tools and samples in a single report. Bioinformatics. 2016;
Ricci V, Ricci V. R: un ambiente opensource per l’analisi statistica dei dati. Econ e Commer [Internet]. 2004;1:69–82. Available from: http://www.dsa.unipr.it/soliani/allegato.pdf
Akalin A, Franke V, Vlahoviček K, Mason CE, Schübeler D. Genomation: A toolkit to summarize, annotate and visualize genomic intervals. Bioinformatics. 2015;31(7):1127–9.
Karolchik D, Hinricks AS, Furey TS, Roskin KM, Sugnet CW, Haussler D, et al. The UCSC table browser data retrieval tool. Nucleic Acids Res. 2004;
Pan Z, Yao Y, Yin H, Cai Z, Wang Y, Bai L, et al. Pig genome functional annotation enhances the biological interpretation of complex traits and human disease. Nat Commun 2021 121 [Internet]. 2021 Oct 6 [cited 2021 Oct 27];12(1):1–15. Available from: https://www.nature.com/articles/s41467-021-26153-7
Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, et al. ClueGO: A Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics. 2009;25(8):191–1093.
Benjamini Y, Hochberg Y. Controlling the False Discovery Rate - a Practical and Powerful Approach to Multiple Testing. Journal of the Royal Statistical Society Series B-Methodological 1995.pdf. J R Stat Soc Ser B. 1995;(June).
Hu ZL, Park CA, Reecy JM. Building a livestock genetic and genomic information knowledgebase through integrative developments of Animal QTLdb and CorrDB. Nucleic Acids Res. 2019;
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

No competing interests reported.

Download PDF

Editorial decision: Major revision
21 Mar, 2022
Reviews received at journal
19 Mar, 2022
Reviews received at journal
16 Mar, 2022
Reviewers agreed at journal
10 Mar, 2022
Reviewers agreed at journal
10 Mar, 2022
Reviewers invited by journal
10 Mar, 2022
Editor assigned by journal
10 Mar, 2022
Submission checks completed at journal
03 Mar, 2022
First submitted to journal
03 Mar, 2022

You are reading this latest preprint version

DNA Methylation May Affect Beef Tenderness Through Signal Transduction in Bos Indicus

Status:

Version 1

Abstract

Background

Results

Conclusion

Figures

1. Background

2. Results

2.1. Animals and Phenotypic Data.

2.2. RRBS DNA sequencing and Nelore muscle methylation profile.

2.3. Most of the differentially methylated CpG and regions were hypermethylated in animals with higher tenderness.

2.4. Genome distribution of differentially methylated cytosines and regions.

2.5. Enrichment analysis of target genes suggests different methylation levels controlling genes acting on signal transduction.

2.6. DMCs and DMRs overlapped QTLs for shear force.

2.7. Manually curated list of candidate methylation regions and cytosines affecting tenderness.

2.8. Four genes had expression correlated with methylation.

2.8.1. GNAS

2.8.2. PDE4B

2.8.3. EBF3

2.8.4. EPCAM

3. Discussion

3.1. A new layer of knowledge regarding tenderness: methylation.

3.2. DMCs and DMRs are more highly methylated in tender than in tougher beef.

3.3. Signal transduction and postmortem processing may be affected by methylation in muscle.

3.4. DMC23 as a regulatory candidate site affecting tenderness through GNAS expression modulation

3.5. Differentially methylated CpG was negatively correlated with PDE4B expression.

3.6. Changes in methylation level ofEBF3 promoter as a potential candidate affecting beef tenderness.

3.7. DMC12 overlapped QTLs for beef tenderness.

3.8. Final considerations.

4. Conclusion

5. Methods

5.1. Experimental design and animals

5.2. Phenotypes and sample selection.

5.3. DNA and RNA extraction and quality evaluation.

5.4. Reduced representation bisulfite sequencing (RRBS) and processing.

5.5. RNA sequencing and processing.

5.6. Differential DNA methylation analysis between extremes.

5.7. Gene annotation and Functional enrichment analysis.

5.8. Development of a manually curated list of candidate methylation regions and cytosines affecting tenderness.

5.9. Correlation between methylation percentage and gene expression.

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1