Figures
Abstract
DNA methylation is an important epigenetic mechanism involved in expression of genes in many biological processes including muscle growth and development. Its effects on economically important traits are evinced from reported significant differences in meat quality traits between Japanese black (Wagyu) and Chinese Red Steppes cattle, thus presenting a unique model for analyzing the effects of DNA methylation on these traits. In the present study, we performed whole genome DNA methylation analysis in the two breeds by whole genome bisulfite sequencing (WGBS). Overall, 23150 differentially methylated regions (DMRs) were identified which were located in 8596 genes enriched in 9922 GO terms, of which 1046 GO terms were significantly enriched (p<0.05) including lipid translocation (GO: 0034204) and lipid transport (GO: 0015914). KEGG analysis showed that the DMR related genes were distributed among 276 pathways. Correlation analysis found that 331 DMRs were negatively correlated with the expression levels of differentially expressed genes (DEGs) with 21 DMRs located in promoter regions. Our results identified novel candidate DMRs and DEGs correlated with meat quality traits, which will be valuable for future genomic and epigenomic studies of muscle development and for marker assisted selection of meat quality traits.
Citation: Fang X, Zhao Z, Yu H, Li G, Jiang P, Yang Y, et al. (2017) Comparative genome-wide methylation analysis of longissimus dorsi muscles between Japanese black (Wagyu) and Chinese Red Steppes cattle. PLoS ONE 12(8): e0182492. https://doi.org/10.1371/journal.pone.0182492
Editor: Sanjay B. Jadhao, International Nutrition Inc, UNITED STATES
Received: December 21, 2016; Accepted: July 7, 2017; Published: August 3, 2017
Copyright: © 2017 Fang et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: The authors state that all data necessary for confirming the conclusions presented in the article are represented fully within the article. Sequence data has been submitted to the GEO (Gene Expression Omnibus) and the accession number is GSE98818.
Funding: This work was supported by the National Natural Science Foundation of China (no. 31372278), the National Major Special Project on New Varieties Cultivation for Transgenis Organisms (2016ZX08009003-006), the National High Technology Research and Development Program (863 Program, no.2013AA102505), the Jilin Scientific and Technological Development Program (No. 20170519014JH), the Jilin province industrial technology research and development program (2016C032) and the Graduate Innovation Fund of Jilin University (Grant number 2015073).
Competing interests: The authors have declared that no competing interests exist.
Introduction
Beef, which is rich in protein, iron, zinc, B vitamins and essential polyunsaturated fatty acids, is a source of high-quality nutrition for humans [1]. However, several diseases, such as heart disease, diabetes and obesity, have been linked with the consumption of fat from beef. Thus, beef breeding program should try to find a balance between eating quality and nutritional quality, which demands a better understanding of the genetic and epigenetic mechanisms controlling meat quality traits.
Advancement in DNA sequencing technology has led to rapid accumulation of data on individual sequence variations in livestock, including cattle [2, 3], which, in combination with the development of high-density bovine genotyping based genome-wide association studies [4, 5], makes it possible for selection of production traits and increased rate of genetic progress in livestock. However, quantitative traits such as meat quality traits are generally determined by many genes, thus the benefit of a particular genetic variant in marker-assisted selection is limited to the proportion of the variant contributing to the heritable phenotypic variance. More importantly, many of the variants detected using GWAS have been in regulatory regions which result in changes in gene expression rather than protein sequences, indicating that the regulatory mechanisms are possibly epigenetics in nature. It is imperative to understand how epigenetics affects phenotypes through the regulation of gene expression. Furthermore, sensory quality of meat is affected by many environmental factors (such as nutrition, stress, exposure to pollution) that likely exert their influence through epigenetic modifications [6].
Epigenetics is the study of heritable changes in gene expression that are not caused by DNA sequence changes [7]. Epigenetic mechanisms include DNA methylation, histone modifications, and non-coding RNAs. DNA methylation plays a central role in regulating many different biological processes, such as growth, development, genomic imprinting and X-chromosome inactivation in females [8]. There is a clear link between abnormal DNA methylation and human diseases [9]. DNA methylation is a critical intermediate molecular phenotype which is considered as the linkage between genotypes and other macro-level phenotypes and maybe attributed to the missing heritability [10]. Current data suggest that genetic variation at specific loci is correlated with the quantitative traits of DNA methylation [11], and genetic variants at CpG sites (meSNPs) could cause changes in their methylation status. It is suggested in human studies that meSNPs attributed a large portion of observed signal from methylation-associated loci (meQTLs) and might be crucial in explaining the results of association between genetic variants and epigenetic changes.
In livestock, genome-wide methylation studies have been performed in chicken, pig, cattle, and sheep, especially in the high-value muscle tissues [6, 12]. The objective of the present study is to perform a comparative genome-wide methylation analysis of longissimus dorsi muscles between the Japanese Black (Wagyu) and Chinese Red Steppes cattle, especially in DMRs and DEGs. Since significant differences exist in meat quality traits between the two breeds, further correlation analysis will allow us to identify DNA methylation variants correlated with meat quality traits, which will be valuable for future genomic and epigenomic studies of muscle development and for marker assisted selection of meat quality traits.
Materials and methods
Ethics statement
Animal care and experiments were performed according to the guideline established by the Regulation for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China, 2004) and approved by the Animal Welfare and Research Ethics Committee at Jilin University (Approval ID: 20140310).
Animals and samples
Three longissimus muscle (LM) tissue samples of Japanese black cattle (28 months old) and three LM samples of Chinese Red Steppes cattle (28 months old) were provided by the National Research Center for Animal Transgenic Bio-technology, Inner Mongolia University (Hohhot, China) and the Branch of Animal Science, Jilin Academy of Agricultural Sciences (Gongzhuling, China), respectively. The three Japanese black cattle were randomly chosen from 20 males born from cows that were artificially inseminated with semen stocks from the same bull. The three Chinese Red Steppes cattle were similarly chosen from 20 males with a common father. The two farms housing the two groups of cattle were located at similar altitudes with similar natural weather conditions. The cattle in both groups were raised under similar normal conditions on a diet of corn and hay with access to feed and water ad libitum. The tissue samples were transported in dry ice and stored in the laboratory in liquid nitrogen.
Extraction of DNA and bisulfite conversion
DNA was extracted from longissimus muscle tissues with a DNA extraction kit (Tiangen, Beijing, China) according to the manufacturer’s protocol and the concentration of the genomic DNA was determined with a NanoDrop ND-2000 UV spectrophotometer (Thermo Fisher Scientific Inc., USA). The degradation and contamination of genomic DNA was monitored with 1% agarose gel electrophoresis.
Library preparation and quantification
Three DNA libraries for each breed were constructed and grouped by breed as WC_LC (Wagyu) and RC_LC (Chinese Red Steppes) for subsequent analysis. 5.2 μg genomic DNA spiked with 26 ng lambda DNA were fragmented to 200–300 bp by sonication with a focused-ultra sonicator (Covaris, Woburn, MA, USA), followed by end repair and adenylation. Cytosine-methylated barcodes were ligated to sonicated DNA according to the manufacturer’s protocol. The DNA bisulfite conversion was performed using the EZ DNA Methylation Gold kit (Zymo Research, California, USA) according to the manufacturer’s instruction. After DNA bisulfite conversion, single-strand DNA fragments were PCR amplified using the KAPA HiFi HotStart Uracil + ReadyMix (2X) (Kapa Biosystems, Wilmington, MA, USA). Library size was quantified by quantitative PCR and a Qubit® 2.0 Flurometer (Life Technologies, CA, USA), and the insert size was checked on an Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).
Whole genome bisulfite sequencing and identification of DMRs
The libraries from 6 samples were sequenced on an Illumina Hiseq2500 platform and 125 nt paired-end reads were generated by Novogene (Novogene, Beijing, China). Read sequences were produced by the Illumina pipeline with FastQ format and were pre-processed by in-house perl scripts. At first, all parts of the 3’ adapter oligonucleotide sequences were filtered out. Secondly, if the percentage of unknown bases (Ns) was greater than 10% the read was removed. Thirdly, reads with low quality (more than 50% of bases with PHRED score< = 5) were trimmed. At the same time, Q20, Q30 and GC content were calculated. All of the subsequent analyses were based on the remaining clean reads that passed the filters.
After filtering the low quality reads, the clean reads data were aligned to the Ensembl bovine reference genome (Bos taurus UMD_3.1.1) and the bisulfite mapping of methylation sites were performed by Bismark (version 0.12.5) [13] with score_min L, 0, -0.2 and then indexed using Bowtie2 [14]. At first, the reference genome and clean reads were transformed into bisulfite-converted version (C-to-T and G-to-A converted). Secondly, sequence read which produced a unique best alignment from the two alignment processes (original top and bottom strand) was then compared to the normal genomic sequence and the methylation status of all cytosine positions in the reads were identified. The reads that aligned to the same regions of genome were regarded as duplicated ones. The sequencing depth and coverage were estimated using deduplicated reads. The data of methylation extractor were transformed into bigWig format for visualization using IGV browser. The sodium bisulfite non-conversion rate of bovine genome was calculated as the percentage of cytosine sequenced at cytosine reference positions in the lambda genome.
To identify the differentially methylated regions (DMRs) between Japanese black cattle (Wagyu) and Chinese Red Steppes, we referenced the modeled of Hao. Y. et al [1] to estimate methylation level [15]. DMRs were identified using a sliding-window approach of swDMR software (http://122.228.158.106/swDMR/) [16] and the window is set to 1000 bp and step length at 100 bp. Fisher test is implemented to detect the DMRs.
GO and KEGG enrichment analysis of genes related to DMRs.
Gene related to DMRs were implemented by the GOseq R package [17], in which gene length bias was corrected. GO terms with corrected p values of less than 0.05 were considered significantly enriched. KOBAS software [18] was used to test the statistical enrichment of DMR-related genes in KEGG pathways [19, 20]. Pathway with a corrected p value <0.05 was considered as significantly enriched.
Sample preparation for RNA‑seq analysis
Total RNA was isolated from longissimus dorsi muscles using Trizol Reagent (Invitrogen, USA) according to the manufacturer's instructions. Total RNA was treated with DNase I (NEB, Beijing, China), extracted with phenol-chloroform, and precipitated with ethanol. The quality and quantity of total RNA were determined using an Agilent 2100 Bioanalyzer (Agilent technologies, Palo Alto, CA). The mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. 1.5 μg of mRNA from each sample was used to construct six cDNA libraries for sequencing. The mRNA was treated with a Thermomixer (Eppendorf AG, Hamburg, Germany) to generate fragments with an average size of 200 bp (200 ± 25 bp) for the paired-end libraries. The fragmented mRNA was then used as templates for synthesizing the first-strand cDNA. The double-stranded cDNAs were purified and ligated to adaptors for Illumina paired-end sequencing. Library concentration was quantified by qPCR and a Qubit® 2.0 Flurometer (Life Technologies, CA, USA), and the insert size was checked on an Agilent Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA). The cDNA libraries were sequenced using the Illumina HiSeq2000 platform by the Beijing Genomics Institute (BGI) and 100 nt paired-end reads were generated.
Quantification of gene expression levels and identification of DEGs
According to our data, the EB-Seq algorithm was applied to filter differentially expressed genes for the WC_LC and RC_LC groups and the URL of the GTF file used to annotate the EB-Seq is ftp://ftp.ensembl.org/pub/release-79/gtf/bos_taurus/Bos_taurus.UMD3.1.79.gtf.gz. After significance analysis and FDR (false discovery rate) analysis[21], we selected the DEGs according to the Log2FC >0.585 or <-0.585, FDR <0.05.
Association analysis
For association analysis, a set of differentially expressed genes with differential methylation was obtained from the intersection between the set of differentially methylated genes and the set of differentially expressed genes. Negative correlations were identified by correlation analysis between the methylation level of DMRs and the expression level of the corresponding genes (r with negative value). Due to the small size of the sample, the present data were not further screened with r2 and p values. In addition, GO analysis was performed on negatively correlated genes in Biological Process, Molecular Function, Cellular Component, their sub-categories and KEGG pathway.
Results
Bisulfite sequencing and DNA methylation profiling
The BS conversion rates of genomic DNA ranged from 99.61% to 99.72%, and a total of 668.37 Gbp of raw sequence data was obtained from 6 samples. After filtering out low quality data, 630.39 Gbp clean sequences were mapped to Bos taurus (UMD_3.1.1) with Q30 of clean full-length reads ranging from 86.15 to 90.84%. High quality methylation maps of the two bovine breeds were obtained and the unique mapping rates ranged from 69.57% to 81.82%. The details of sequencing data quality were shown in Table 1. The results of WGBS analyses identified differentially methylated regions (DMRs) covering almost the entire genome with sufficient depth and high resolution. The average distribution coverage of the genome in 6 samples was shown in S1 Fig.
CpGs, CHHs and CHGs (H = C, T and A) were methylated at different levels. 48.27–64.61% of CpGs in the whole genome were methylated while only 0.07–0.10% and 0.8–0.13% cytosines in the CHG and CHH contexts were methylated, respectively (Table 2). The distribution details of methylated cytosines in the contexts of CpGs, CHHs and CHGs in whole genome were shown in S1 Table.
To better understand DNA methylation levels at various functional genomic elements, such as promoters, exons, introns, etc., the 2 kb region upstream of the TSS point (defined as the promoter region) of each gene was equally divided into 20 bins, and the corresponding bin level on C locus of all genes were averaged. For mC sites in each context, the average methylation levels and the genome-wide landscape of distribution of various functional genomic elements of 6 samples were built by ggplot2 R package[22] as shown in Fig 1. Similar tendencies of methylation levels in all samples were observed in different genetic elements in the three mC contexts, without any significant differences among them. Overall, the DNA methylation levels in introns were the highest, followed by exons and 3’ UTR, with 5’ UTR showing the lowest level. In promoter regions, the proximal regions had the lowest DNA methylation levels, followed by the intermediate regions, with the distal regions having the highest levels. As for mC in CHH context, the 5’ proximal regions of intron had the highest methylation levels in all samples which were different from mCs in other contexts.
Abscissa represents different genomic elements, the ordinate represents the level of methylation and different colors represents different sequence contexts (CpG, CHG, CHH). WC_LM1, WC_LM2, WC_LM3: longissimus dorsi muscle samples of Japanese black cattle (Wagyu); RC_LM1, RC_LM2, RC_LM3: longissimus dorsi muscle samples of Chinese Red Steppes.
The difference and trend in genomic methylation levels between samples can be analyzed by regional cluster classification analysis to understand the biological significance [23]. Using hierarchical clustering method, we analyzed the whole genomic methylation in Japanese black cattle (Wagyu) and Chinese Red Steppes cattle, as shown in the heat maps of Fig 2, which revealed a clear separation between the two breeds. Our results also showed that the biological variation among the three samples within the same breed was extremely low. In addition, the cluster analysis results for the genomic methylation in different functional regions were summarized in S2 Fig.
Heat map displays highly methylated loci in red and sparsely methylated loci in blue.
Characterization of DMRs
In total, we compared the methylation levels of longissimus dorsi muscle from the two breeds and identified 23,150 DMRs in gene bodies or promoter regions from 8,596 genes that were significantly different between Wagyu and Chinese Red Steppes (corrected p< 0.05) which have different meat quality traits. Among the DMRs, 11,943 corresponded to 6,072 genes with higher methylation levels in Wagyu than that in Chinese Red Steppes, whereas 11,207 corresponded to 5,034 genes with lower methylation levels in Wagyu (Fig 3). The top 20 DMRs were listed in Table 3 by ascending order of corrected p value. Interestingly, 2,510 genes possessed multiple DMRs with both hyper-methylated (higher methylation level in Wagyu than in Chinese Red Steppes) regions and hypo-methylated (lower methylation level in Wagyu than in Chinese Red Steppes) regions, such as acyl-CoA synthetase long-chain family member 6 (ACSL6), epidermal growth factor receptor (EGFR), 1-acylglycerol-3-phosphate O-acyltransferase 4 (AGPAT4), upstream binding protein 1(UBP1).
(A) Numbers of DMRs; (B) Numbers of genes related to DMRs. DMRs: differentially methylated regions.
The length of DMRs ranged between 3 nt to 3,780 nt (S3 Fig). Most of the DMRs were located at introns (73.66%), followed by exons (14.36%) and regulatory regions (11.98%), such as promoters, 5’UTRs and 3’UTRs, as shown in Table 4.
Gene ontology and KEGG enrichment analyses of genes related to DMRs.
The identified DMRs were functionally classified by GO analysis. Functional enrichments on 23,150 DMRs (corrected p< 0.05) by gene ontology analysis found that 8,596 genes related to DMRs were enriched in 9,922 GO terms, among which 1,046 GO terms were significantly enriched (corrected p< 0.05). The most significant functional terms were binding (GO: 0005488) of hyper-methylated genes and protein binding (GO: 0005515) of hypo-methylated genes. In terms of biological process, the most significant functional terms were localization (GO: 0051179) of hyper-methylated genes and anatomical structure morphogenesis (GO: 0009653) of hypo-methylated genes. In cellular component, the most significant functional terms were plasma membrane parts (GO: 0044459) of hyper-methylated genes and cell projection (GO: 0042995) of hypo-methylated genes. The top 30 GO terms were listed in Fig 4 by ascending order of corrected p value.
The ordinate represents the enriched GO terms; the abscissa represents the number of genes; “*” represents the GO term significantly enriched (corrected p< 0.05) (only top 30 GO terms are listed by ascending order of corrected p values).
KEGG enrichment analysis of genes related to DMRs was also conducted for both hyper-methylated and hypo-methylated genes. Overall, 276 pathways were identified, such as calcium signaling pathway (bta04020), focal adhesion (bta04510), cAMP signaling pathway (bta04024), and cGMP-PKG signaling pathway (bta04022), and the top 20 in ascending order of corrected p value were listed in Fig 5. Results showed that calcium signaling pathway (bta04020, corrected p = 0.002), glutamatergic synapse (bta04724, corrected p = 0.047) and cAMP signaling pathway (bta04024, corrected p = 0.047) were significantly enriched in hyper-methylated genes (corrected p< 0.05), with calcium signaling pathway showed the highest enrichment (S4 Fig). However, no pathway was significantly enriched in hypo-methylated genes (corrected p> 0.05).
The ordinate represents the enriched pathways, and the abscissa represents the rich factor of corresponding pathways; the size of the spots represents the number of genes related to DMRs enriched in each pathway, while the color of the spot represents the corrected p value of each pathway. The rich factors indicate the ratio of the number of DMGs mapped to a certain pathway to the total number of genes mapped to this pathway. Greater rich factor means greater enrichment. DMRs: differentially methylated regions; DMGs: differentially methylated genes.
Association analysis between DMRs and DEGs
We next analyzed the correlation between the methylation levels of DMRs and the expression levels of DEGs. The results indicated that 147 DMR related genes had negative correlation with the expression levels of DEGs (Table 5), of which 19 genes had DMRs in promoter regions, 50 genes had DMRs in exons, 129 genes had DMRs in introns and only 7 genes and 5 genes had DMRs in 3’UTRs and 5’UTRs, respectively (Fig 6). For genes with negative correlation, 56 genes were enriched in 969 GO terms, of which 259 GO terms were significant (corrected p< 0.05) including positive regulation of fat cell differentiation (GO: 0045600, corrected p = 0.005), glycerol metabolic process (GO: 0006071, corrected p = 0.010), negative regulation of low-density lipoprotein particle receptor catabolic process (GO: 0032804, corrected p = 0.016), triglyceride metabolic process (GO: 0006641, corrected p = 0.046), negative regulation of lipoprotein lipase activity (GO: 0051005, corrected p = 0.031) and glycerophospholipid metabolic process (GO: 0006650, corrected p = 0.039). In addition, the 56 genes were enriched in 118 pathways, of which 11 pathways were significant (corrected p< 0.05) such as ECM-receptor interaction (bta04512, corrected p = 0.005), propanoate metabolism (bta00640, corrected p = 0.029), cGMP-PKG signaling pathway (bta04022, corrected p = 0.032), carbon metabolism (bta01200, corrected p = 0.044), fatty acid degradation (bat00071, corrected p = 0.047) (Fig 7). Interestingly, 13 genes with negative correlation were enriched in metabolic pathway term (bta01100). Among negatively correlated genes, 6 genes with DMRs in promoter region were enriched in 10 pathways in which 4 pathways were significantly enriched (corrected p< 0.05) (S5 Fig).
DMRs: differentially methylated regions; DEGs: differentially expressed genes.
Abscissa represents the value of log2 (p value) and ordinate represents the pathway name; Red bars indicate corrected p< 0.05. DMRs: differentially methylated regions; DEGs: differentially expressed genes.
Discussion
DMR is an important epigenetic marker which may be involved in regulation of gene expression by changing chromatin structure or transcription efficiency under different conditions [24]. Although DNA methylation analysis has been performed on cattle [12, 25, 26], this study is the first to systematically compare the genome-wide methylation profiles in longissimus dorsi muscle between Japanese black cattle and Chinese Red Steppes, between which significant difference exist in overall meat quality traits.
DNA methylation profiles in longissimus dorsi muscle
The genome-wide methylation patterns in functional genomic regions were very similar between the two breeds. However, there were differences among three mC contexts which might be related to the sequence differences in different genetic elements. The majority of DMRs were small fragments (50–1000 nt > 90% of DMRs), which suggested that methylation changes within a small regions could play an important role on the regulation of gene expression. There were more DMRs with higher methylation level and less DMRs with lower methylation in Wagyu in comparison to Chinese Red Steppes, and DMRs were mainly concentrated in the intron regions (> 70%) with only a small proportion distributed in the 3’ UTRs and 5’ UTRs. In addition, the distribution of certain DMRs in either hyper-methylated or hypo-methylated status with the same gene indicated that different methylation status in different regions might have different regulatory functions on gene expression.
Correlation analysis between DMRs and signal pathways
Analyses of pathways and genes related to DMRs revealed that the hyper-methylated genes were significantly concentrated in three pathways: calcium signaling pathway (bta04020, corrected p = 0.002), glutamatergic synapse (bta04724, corrected p = 0.047) and cAMP signaling pathway (bta04024, corrected p = 0.047) (S6 Fig). The KEGG enrichment analysis indicated that calcium signaling pathway was the most significantly correlated pathway, which is the key pathway exerting allosteric regulation on many proteins by the calcium ions signaling effects, such as activation of ion channels or as a secondary messenger.
It is reported that Caveolin 1 (CAV1) gene has a higher expression level in adipocytes [27], and CAV1 gene knockout leads to imbalance in lipid metabolism [28]. Phospholipase C (PLC) cleaves phospholipids into free fatty acids (FFA), lysophospholipids and diacylglycerol, and the relative activities of PLC are highly correlated with the decrease of phospholipids and the increase of free fatty acids [29]. These studies indicated that CAV1 and PLC genes are involved in lipid and fatty acid metabolism. A mutation in the ryanodine receptor (RyR) gene is found to be correlated with malignant hyperthermia in lean, heavily muscled swine breeds [30]. In this study, CAV1, PLC and RYR genes were identified as differentially methylated genes between the two breeds, which suggested that certain differentially methylated genes in calcium signaling pathway could affect the metabolism of the skeletal muscles. Peroxisome proliferator-activated receptor gamma (PPARγ), a key molecule found in cAMP signal pathway involving the G protein coupled receptors [31], is involved in adipocyte differentiation and function. The present study identified three types of PPARs: alpha, gamma, and delta [32]. PPARγ is expressed in liver, kidney, heart, muscle, adipose tissue, and others [33]. As a dietary lipid sensor, the activation of PPARγ results in lipid metabolism in muscle [34]. Interestingly, all the three pathways could be regulated by calcium ions concentration, which indicated that this might be the core pathway in the regulation of muscle metabolism [35].
Negative correlation between the methylation levels of DMRs and differential expressions of DEGs
There is a complex relationship between gene methylation and gene expression level. Current research data suggest that methylation in the promoter region blocks transcription initiation, but the function on gene expression of DNA methylation within gene body has not been clearly understood [24].
As seen in Table 5, the DMRs in the promoter and gene body regions were negatively correlated with gene expression level. Enoyl-CoA, hydratase/3-hydroxyacyl CoA dehydrogenase (EHHADH) and 4-N-trimethylaminobutyraldehyde dehydrogenase (ALDH9A1) genes were enriched in the pathway of fatty acid degradation (p = 0.047, p< 0.05), during which fatty acids are broken down into acetyl-CoA to enter into the citric acid cycle for energy production in animals [36]. It seemed that DMRs related to genes in the fatty acid degradation pathway might affect fatty acid composition in meat through regulation of gene expression levels. EHHADH gene was identified by GWAS as one of the 20 promising novel genes associated with milk fatty acid traits in Chinese Holstein [37]. Furthermore, EHHADH and angiopoietin-like protein 4 (ANGPTL4) genes both belong to the PPAR signaling pathway (bta03320), which is a key biological pathway in determining meat quality traits in mammals by involving in lipid metabolism. ANGPTL4 is mainly expressed in liver and adipose tissues [38] and plays important roles in the regulation of lipid metabolism. Previous studies have suggested that a SNP in ANGPTL4 is associated with intramuscular fat [39] and the expression level of ANGPTL4 affects meat quality traits of pigs [40]. Recent studies have shown that acyl-CoA synthetase family member 3 (ACSF3), in which DMRs in promoter region were found in the present study to be negatively correlated with the gene expression level, plays a critical role in mammal fatty acid synthesis in mitochondrion [41]. Thus EHHADH, ALDH9A1, ANGPTL4 and ACSF3 play important roles in fatty acid and lipid metabolism in muscle and regulation of these genes through DNA methylation might be part of the mechanisms in determining the difference in meat quality traits between the two breeds.
The classical theory on animal breeding believes that the appearance of a quantitative trait is the result of multiple genetic loci, in particular, a single polymorphic locus with multiple, differentially expressed alleles can give rise to the quantitative trait variation within a natural population [42]. Many of the differentially methylated genes identified in the present study between Wagyu and Chinese Red Steppes are also found to be differentially methylated in genes related to lipid metabolism and fatty acid metabolism, such as fatty acid synthase (FASN) which is shown to be related to fatty acid composition in a genome-wide association study of Japanese black cattle [43], and carnitine palmitoyltransferase 1A (CPT1A), fatty acid desaturase 1 (FADS1) and acyl-CoA synthetase long-chain family member 1 (ACSL1) genes which are well-known genes in the pathways of fatty acid metabolism and fatty acid degradation. In addition, other genes known to affect meat quality traits of cattle are also related to the lipid metabolism and fatty acid metabolism (insulin like growth factor 2 (IGF2), leptin (OB), lipase hormone sensitive (HSL), diacylglycerol O-acyltransferase 1 (DGAT1) and fatty acid binding protein 3 (FABP3)) [44–49]. We believed that the methylation of these genes might partially contribute to the significant variance in meat quality traits between Japanese black cattle and Chinese Red Steppes. However, the regulatory epigenetic mechanisms of these genes and genetic regions on bovine meat quality traits require further study.
Conclusion
In the present study, whole genome, high resolution DNA methylation profiles of longissimus dorsi muscles were established for Japanese black and Chinese Red Steppes cattle. Combined with transcriptomic analysis, our study investigated the genome-wide regulation of gene expression by DNA methylation at transcription level, and identified a number of novel and important genes associated with DMRs and pathways that might affect muscle development and meat quality traits in cattle, which needed to be further experimentally validated in the future. Our results provided valuable data to further our understanding of the genetic and epigenetic mechanisms that control economic traits in cattle, which could be used in marker-assisted selection programs to improve beef production.
Supporting information
S1 Fig. Distribution of genomic sequencing depth.
(A) The distribution of sequence coverage at the genome of all samples, abscissa represents the depth of coverage and ordinate represents its frequency. (B) The accumulated distribution of sequence coverage at genome of all samples, abscissa represents the depth of coverage and ordinate represents the ratio of coverage not less than the depth of the total number of base sites; different lines with different colors represent different samples.
https://doi.org/10.1371/journal.pone.0182492.s001
(TIF)
S2 Fig. Clusters of different functional regions by methylation levels.
https://doi.org/10.1371/journal.pone.0182492.s002
(TIF)
S4 Fig. Pathways enriched for genes related to DMRs with higher methylation levels in Japanese black cattle (Wagyu) than that in Chinese Red Steppes.
The ordinate represents the enriched pathways, and the abscissa represents the rich factor of corresponding pathways; the size of the spots represented the number of genes related to DMRs enriched in each pathway, while the color of the spot represents the corrected p value of each pathway. The rich factors indicate the ratio of the number of DMGs mapped to a certain pathway to the total number of genes mapped to this pathway. Greater rich factor means greater enrichment. DMRs: differentially methylated regions; DMGs: differentially methylated genes.
https://doi.org/10.1371/journal.pone.0182492.s004
(TIF)
S5 Fig. Pathway enrichment of genes with negative correlation between methylation levels of DMRs in promoter and expression levels of mRNA.
Abscissa represents the value of log2 (p value) and ordinate represents the pathway name; red bars indicate corrected p< 0.05.
https://doi.org/10.1371/journal.pone.0182492.s005
(TIF)
S6 Fig. KEGG pathways with genes related to significantly enriched DMRs with higher methylation in Wagyu.
Genes with red marker are related to DMRs with higher methylation levels in Japanese black cattle (Wagyu) than that in Chinese Red Steppes. DMRs: differentially methylated regions.
https://doi.org/10.1371/journal.pone.0182492.s006
(TIF)
S1 Table. Comparison of DNA methylation patterns between the two groups.
https://doi.org/10.1371/journal.pone.0182492.s007
(DOCX)
Acknowledgments
We are grateful to the Branch of Animal Science, Jilin Academy of Agricultural Sciences and National Research Center for Animal Transgenic Bio-technology, Inner Mongolia University for animal sampling. The authors thank Novel Bioinformatics Co., Ltd for technical help and KEGG for providing pathway figures in S6 Fig.
References
- 1. McNeill S, Van Elswyk ME. Red meat in global nutrition. Meat science. 2012;92(3):166–73. Epub 2012/06/05. pmid:22658075.
- 2. Ramos AM, Crooijmans RP, Affara NA, Amaral AJ, Archibald AL, Beever JE, et al. Design of a high density SNP genotyping assay in the pig using SNPs identified and characterized by next generation sequencing technology. PLoS One. 2009;4(8):e6524. Epub 2009/08/06. pmid:19654876; PubMed Central PMCID: PMCPMC2716536.
- 3. Stothard P, Choi JW, Basu U, Sumner-Thomson JM, Meng Y, Liao X, et al. Whole genome resequencing of black Angus and Holstein cattle for SNP and CNV discovery. BMC genomics. 2011;12:559. Epub 2011/11/17. pmid:22085807; PubMed Central PMCID: PMCPMC3229636.
- 4. Ma L, Runesha HB, Dvorkin D, Garbe JR, Da Y. Parallel and serial computing tools for testing single-locus and epistatic SNP effects of quantitative traits in genome-wide association studies. BMC bioinformatics. 2008;9:315. Epub 2008/07/23. pmid:18644146; PubMed Central PMCID: PMCPMC2503991.
- 5. Matukumalli LK, Lawley CT, Schnabel RD, Taylor JF, Allan MF, Heaton MP, et al. Development and characterization of a high density SNP genotyping assay for cattle. PLoS One. 2009;4(4):e5350. Epub 2009/04/25. pmid:19390634; PubMed Central PMCID: PMCPMC2669730.
- 6. Couldrey C, Brauning R, Bracegirdle J, Maclean P, Henderson HV, McEwan JC. Genome-wide DNA methylation patterns and transcription analysis in sheep muscle. PLoS One. 2014;9(7):e101853. Epub 2014/07/11. pmid:25010796; PubMed Central PMCID: PMCPMC4092064.
- 7. Egger G, Liang G, Aparicio A, Jones PA. Epigenetics in human disease and prospects for epigenetic therapy. Nature. 2004;429(6990):457–63. Epub 2004/05/28. pmid:15164071.
- 8. Smith ZD, Meissner A. DNA methylation: roles in mammalian development. Nature reviews Genetics. 2013;14(3):204–20. Epub 2013/02/13. pmid:23400093.
- 9. Robertson KD. DNA methylation and human disease. Nature reviews Genetics. 2005;6(8):597–610. Epub 2005/09/02. pmid:16136652.
- 10. Maher B. Personal genomes: The case of the missing heritability. Nature. 2008;456(7218):18–21. Epub 2008/11/07. pmid:18987709.
- 11. Zhi D, Aslibekyan S, Irvin MR, Claas SA, Borecki IB, Ordovas JM, et al. SNPs located at CpG sites modulate genome-epigenome interaction. Epigenetics. 2013;8(8):802–6. Epub 2013/07/03. pmid:23811543; PubMed Central PMCID: PMCPMC3883783.
- 12. 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). Scientific reports. 2014;4:6546. Epub 2014/10/14. pmid:25306978; PubMed Central PMCID: PMCPMC4194443.
- 13. Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics (Oxford, England). 2011;27(11):1571–2. Epub 2011/04/16. pmid:21493656; PubMed Central PMCID: PMCPMC3102221.
- 14. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nature Methods. 2012;9(4):357. pmid:22388286
- 15. Hao Y, Cui Y, Gu X. Genome-wide DNA methylation profiles changes associated with constant heat stress in pigs as measured by bisulfite sequencing. Scientific Reports. 2016;6.
- 16. Wang Z, Li X, Jiang Y, Shao Q, Liu Q, Chen BY, et al. swDMR: A Sliding Window Approach to Identify Differentially Methylated Regions Based on Whole Genome Bisulfite Sequencing. PloS one. 2015;10(7).
- 17. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11(2):R14. Epub 2010/02/06. pmid:20132535; PubMed Central PMCID: PMCPmc2872874.
- 18. Mao X, Cai T, Olyarchuk JG, Wei L. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics. 2005;21(19):3787–93. Epub 2005/04/09. pmid:15817693.
- 19. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic acids research. 2000;28(1):27–30. Epub 1999/12/11. pmid:10592173; PubMed Central PMCID: PMCPmc102409.
- 20. Kanehisa M, Sato Y, Morishima K. BlastKOALA and GhostKOALA: KEGG Tools for Functional Characterization of Genome and Metagenome Sequences. Journal of molecular biology. 2016;428(4):726–31. Epub 2015/11/21. pmid:26585406.
- 21. Benjamini Y, Drai D, Elmer G, Kafkafi N, Golani I. Controlling the false discovery rate in behavior genetics research. Behavioural brain research. 2001;125(1–2):279–84. Epub 2001/10/30. pmid:11682119.
- 22.
Wickham H. ggplot2: Elegant Graphics for Data Analysis: Springer Publishing Company, Incorporated; 2009. 180–5 p.
- 23. Smallwood SA, Lee HJ, Angermueller C, Krueger F, Saadeh H, Peat J, et al. Single-cell genome-wide bisulfite sequencing for assessing epigenetic heterogeneity. Nat Methods. 2014;11(8):817–20. Epub 2014/07/22. pmid:25042786; PubMed Central PMCID: PMCPmc4117646.
- 24. Jones PA. Functions of DNA methylation: islands, start sites, gene bodies and beyond. Nature reviews Genetics. 2012;13(7):484–92. Epub 2012/05/30. pmid:22641018.
- 25. Salilew-Wondim D, Fournier E, Hoelker M, Saeed-Zidane M, Tholen E, Looft C, et al. Genome-Wide DNA Methylation Patterns of Bovine Blastocysts Developed In Vivo from Embryos Completed Different Stages of Development In Vitro. PloS one. 2015;10(11):e0140467. Epub 2015/11/05. pmid:26536655; PubMed Central PMCID: PMCPmc4633222.
- 26. Mendonca Ados S, Guimaraes AL, da Silva NM, Caetano AR, Dode MA, Franco MM. Characterization of the IGF2 Imprinted Gene Methylation Status in Bovine Oocytes during Folliculogenesis. PloS one. 2015;10(10):e0142072. Epub 2015/10/31. pmid:26517264; PubMed Central PMCID: PMCPmc4627647.
- 27. Li WP, Liu P, Pilcher BK, Anderson RG. Cell-specific targeting of caveolin-1 to caveolae, secretory vesicles, cytoplasm or mitochondria. Journal of Cell Science. 2001;114(Pt 7):1397–408. pmid:11257005
- 28. Razani B, Combs TP, Wang XB, Frank PG, Park DS, Russell RG, et al. Caveolin-1-deficient mice are lean, resistant to diet-induced obesity, and show hypertriglyceridemia with adipocyte abnormalities. The Journal of biological chemistry. 2002;277(10):8635–47. Epub 2001/12/12. pmid:11739396.
- 29. Wang D, Zhang M, Bian H, Xu W, Xu X, Zhu Y, et al. Changes of phospholipase A(2) and C activities during dry-cured duck processing and their relationship with intramuscular phospholipid degradation. Food chemistry. 2014;145:997–1001. Epub 2013/10/17. pmid:24128575.
- 30. Fujii J, Otsu K, Zorzato F, Leon SD, Khanna VK, Weiler JE, et al. Identification of a mutation in porcine ryanodine receptor associated with malignant hyperthermia. Science. 1991;253(5018):448–51. pmid:1862346
- 31. Gilman AG. G Proteins: Transducers of Receptor-Generated Signals. Annual Review of Biochemistry. 1987;56(56):615–49.
- 32.
Berger J, Moller DE. The mechanisms of action of PPARs: Cold Spring Harbor Laboratory. 1249–59 p.
- 33. Tyagi S, Gupta P, Saini AS, Kaushal C, Sharma S. The peroxisome proliferator-activated receptor: A family of nuclear receptors role in various diseases. Journal of advanced pharmaceutical technology & research. 2011;2(4):236–40. Epub 2012/01/17. pmid:22247890; PubMed Central PMCID: PMCPMC3255347.
- 34. Ebrahimi M, Rajion MA, Goh YM. Effects of oils rich in linoleic and alpha-linolenic acids on fatty acid profile and gene expression in goat meat. Nutrients. 2014;6(9):3913–28. Epub 2014/09/26. pmid:25255382; PubMed Central PMCID: PMCPmc4179195.
- 35. Maclennan DH, Ms P. Malignant Hyperthermia. Science. 1992;256(5058):789–94. pmid:1589759
- 36.
Máximo CO. Fatty acid degradation. Fatty Acids. 2011.
- 37. Li C, Sun D, Zhang S, Wang S, Wu X, Zhang Q, et al. Genome wide association study identifies 20 novel promising genes associated with milk fatty acid traits in Chinese Holstein. PLoS One. 2014;9(5):e96186. Epub 2014/05/27. pmid:24858810; PubMed Central PMCID: PMCPMC4032272.
- 38. Mamedova LK, Robbins K, Johnson BJ, Bradford BJ. Tissue expression of angiopoietin-like protein 4 in cattle. Journal of animal science. 2010;88(1):124–30. Epub 2009/09/29. pmid:19783696.
- 39. Ren ZQ, Wu WJ, Liu WH, Zheng R, Li JL, Zuo B, et al. Differential expression and effect of the porcine ANGPTL4 gene on intramuscular fat. Genetics and molecular research: GMR. 2014;13(2):2949–58. Epub 2014/05/02. pmid:24782129.
- 40. Wang J, Zhu X, Wang Z, Yao J, Zhao B, Liu G. Non-esterified fatty acids promote expression and secretion of angiopoietin-like protein 4 in calf hepatocytes cultured in vitro. Molecular and cellular biochemistry. 2015;401(1–2):141–6. Epub 2014/12/17. pmid:25501302.
- 41. Witkowski A, Thweatt J, Smith S. Mammalian ACSF3 protein is a malonyl-CoA synthetase that supplies the chain extender units for mitochondrial fatty acid synthesis. The Journal of biological chemistry. 2011;286(39):33729–36. Epub 2011/08/19. pmid:21846720; PubMed Central PMCID: PMCPmc3190830.
- 42.
Koning DJd. Identification of (non-) Mendelian factors affecting pork production. [S.l.: s.n.]; 2001.
- 43. Sasago N, Abe T, Sakuma H, Kojima T, Uemoto Y. Genome-wide association study for carcass traits, fatty acid composition, chemical composition, sugar, and the effects of related candidate genes in Japanese Black cattle. Animal science journal = Nihon chikusan Gakkaiho. 2016. Epub 2016/04/27. pmid:27112906.
- 44. Fang XB, Zhang LP, Yu XZ, Li JY, Lu CY, Zhao ZH, et al. Association of HSL gene E1-c.276C>T and E8-c.51C>T mutation with economical traits of Chinese Simmental cattle. Molecular biology reports. 2014;41(1):105–12. Epub 2013/11/12. pmid:24213829.
- 45. Goszczynski DE, Mazzucco JP, Ripoli MV, Villarreal EL, Rogberg-Munoz A, Mezzadra CA, et al. Characterization of the bovine gene LIPE and possible influence on fatty acid composition of meat. Meta gene. 2014;2:746–60. Epub 2015/01/22. pmid:25606458; PubMed Central PMCID: PMCPMC4287880.
- 46. Sherman EL, Nkrumah JD, Murdoch BM, Li C, Wang Z, Fu A, et al. Polymorphisms and haplotypes in the bovine neuropeptide Y, growth hormone receptor, ghrelin, insulin-like growth factor 2, and uncoupling proteins 2 and 3 genes and their associations with measures of growth, performance, feed efficiency, and carcass merit in beef cattle. Journal of animal science. 2008;86(1):1–16. Epub 2007/09/06. pmid:17785604.
- 47. Han RH, Zan LS, Yang DP, Hao RC. [SNPs detection of IGF2 gene and its relationship with carcass and meat quality traits in Qinchuan cattle]. Yi chuan = Hereditas / Zhongguo yi chuan xue hui bian ji. 2008;30(12):1579–84. Epub 2008/12/17. pmid:19073573.
- 48. Aviles C, Polvillo O, Pena F, Juarez M, Martinez AL, Molina A. Associations between DGAT1, FABP4, LEP, RORC, and SCD1 gene polymorphisms and fat deposition in Spanish commercial beef. Journal of animal science. 2013;91(10):4571–7. Epub 2013/08/24. pmid:23965384.
- 49. de Oliveira JA, da Cunha CM, Crispim Bdo A, Seno Lde O, Fernandes AR, Nogueira Gde P, et al. Association of the leptin gene with carcass characteristics in Nellore cattle. Animal biotechnology. 2013;24(3):229–42. Epub 2013/06/20. pmid:23777351.