Skip to main content

Diverse WGBS profiles of longissimus dorsi muscle in Hainan black goats and hybrid goats



Goat products have played a crucial role in meeting the dietary demands of people since the Neolithic era, giving rise to a multitude of goat breeds globally with varying characteristics and meat qualities. The primary objective of this study is to pinpoint the pivotal genes and their functions responsible for regulating muscle fiber growth in the longissimus dorsi muscle (LDM) through DNA methylation modifications in Hainan black goats and hybrid goats.


Whole-genome bisulfite sequencing (WGBS) was employed to scrutinize the impact of methylation on LDM growth. This was accomplished by comparing methylation differences, gene expression, and their associations with growth-related traits.


In this study, we identified a total of 3,269 genes from differentially methylated regions (DMR), and detected 189 differentially expressed genes (DEGs) through RNA-seq analysis. Hypo DMR genes were primarily enriched in KEGG terms associated with muscle development, such as MAPK and PI3K-Akt signaling pathways. We selected 11 hub genes from the network that intersected the gene sets within DMR and DEGs, and nine genes exhibited significant correlation with one or more of the three LDM growth traits, namely area, height, and weight of loin eye muscle. Particularly, PRKG1 demonstrated a negative correlation with all three traits. The top five most crucial genes played vital roles in muscle fiber growth: FOXO3 safeguarded the myofiber’s immune environment, FOXO6 was involved in myotube development and differentiation, and PRKG1 facilitated vasodilatation to release more glucose. This, in turn, accelerated the transfer of glucose from blood vessels to myofibers, regulated by ADCY5 and AKT2, ultimately ensuring glycogen storage and energy provision in muscle fibers.


This study delved into the diverse methylation modifications affecting critical genes, which collectively contribute to the maintenance of glycogen storage around myofibers, ultimately supporting muscle fiber growth.

Peer Review reports


The domestication of goats (Capra hircus) has a long history dating back to the Neolithic era, providing humans with essential resources such as meat, skin, and milk [1]. Goat meat, known for its nutrient-rich composition characterized by low lipids, low cholesterol, and high-quality proteins, has garnered increasing attention from consumers due to its quality. Over time, goat populations have expanded significantly during the process of domestication, resulting in the rapid emergence of numerous goat breeds worldwide. Some of these breeds have demonstrated rapid growth, while others are renowned for their flavors and nutritional value [2]. Throughout the breeding process, it has become evident that cross bred goats inherit advantageous traits from their parent breeds. One key factor influencing muscle quality is glycogen, a crucial component for myocyte growth and nutritional maintenance. Notably, excessive glycogenolysis during slaughter leads to the production of pale, soft, and exudative meat [3]. It has also been recognized that, methylation can impact muscle development [4], although the specific effect of methylation on myocyte growth mediated by glycogen storage in cross bred goats remains unclear.

Muscle glycogen serves not only as an energy source, but also as a regulator of whole-body energy metabolism [5], For instance, bodybuilders manipulate their carbohydrate intake to increase muscle thickness, thus enhancing cross-sectional area, muscle volume, and overall body weight [6]. Conversely, a decrease in muscle glycogen content, as seen during the post-transportation rest period in alpacas before slaughter, can lead to increased drip loss and purge [7], underscoring the significance of glycogen for maintaining meat quality. Moreover, methylation has been linked to porcine meat pH [8], and lamb meat quality [4]. Whole-genome bisulfite sequencing (WGBS) has proven to be a valuable tool for exploring methylation in various tissues due to its ability to provide comprehensive information about methylcytosine [9]. Methylation levels can impact numerous physical functions, including growth, productive performance, and meat quality [10], as well as muscle development [11]. Furthermore, hypomethylated regions found in promoters and enhancers have the potential to enhance the overall health of ovine tissues, leading to increased production of meat, milk, and wool [12]. The GTPase regulator activation has been identified as valuable in differentially methylated regions (DMRs), as it is closely associated with oxidative and glycolytic muscles [13]. Importantly, differences in DNA methylation levels between breeds can result in various production capacities, such as milk yield and muscle development [14]. Methylation has also been linked to meat tenderness through its influence on alternative splicing of muscle genes [15]. Conversely, DNA methylation might cause insulin resistance by affecting lipid droplets and metabolism genes in muscle [16].

Nubian goats, widely distributed in arid regions of Northern Sudan [17], are characterized by rapid growth, excellent meat quality, and high meat production [18]. These qualities make them superior to many local goat breeds in other countries, leading to their use in hybridization programs to improve growth traits. For example, the F1 hybrid offspring resulting from the crossing of Yunling black goats and Nubian goats, have exhibited specific genes expression related to bone development, muscle cell differentiation, lipid metabolism, and adaptation to hot and humid environments [19]. In contrast, Hainan black goats, a unique local breed exclusive to the tropical island in China, have evolved under the distinct climatic conditions of high temperature and high humidity. They are known for their tender meat, high nutritional value, and disease resistance [20]. However, their growth rate and body size are smaller than those of Nubian goats. Hainan black goats typically reach sexual maturity between 4 and 6 months of age and only give birth to one kid per year [21]. Nubian goats were chosen as ideal crossbreeding partners for Hainan black goats to address their subpar reproductive performance and preserve high-quality meat traits.

Furthermore, hybridization involving geographically different species has yielded significant growth advantages in recent years, with hybrid goats (Nubian goats × Hainan black goats), displaying significantly greater weight compared to Hainan black goats at the same age. While the concept of heterosis is well-established, the genetic mechanism responsible for the expression of desirable traits remains unclear, as do the effects of DNA methylation on skeletal muscle and gene expression patterns. Compared with growth traits of other livestock species, such as cattle [22], pigs [23], and sheep [24], the investigation of DNA methylation in goat has been relatively limited to reproductive performance [25, 26]. Therefore, it is imperative to explore the regulatory mechanisms related to DNA methylation in goat growth traits, particularly in the context of glycogen storage in hybrid goat muscle. The primary aim of this study is to identify the key genes responsible for regulating the growth of the longissimus dorsi muscle (LDM) and to assess the effects of DNA methylation on muscle glycogen storage in hybrid goats.


Identification and distribution of mC levels

WGBS was conducted on LDM tissues from six individuals, and yielded approximate 553.21 Gb of raw data across 1,844,041,561 raw reads, with an average sequencing depth of ~30× for each sample (Table S1). After trimming adaptors and removing low-quality reads, the clean data were mapped to the reference genome (ARS1) to obtain an average mapping ratio of 81.08% (Table S2). The sequencing depth was most often seen to be 24× (Fig. S1), with the highest coverage being 50× (Fig. S2). In general, the methylation density and level of the three mC contexts (CpG, CHG, and CHH) were consistent in LDM between Hainan black goats and hybrid goats (Fig. 1A and B). The CG context occupied a greater relevance on chromosomes 1 compared to other chromosomes (Fig. S3), meanwhile, the CG sites exhibited the largest number and highest methylation proportion across the three mC contexts in Hainan black goats and hybrid goats (Fig. 1C and Fig. S4), as well as methylation levels (Fig. S5).

Fig. 1
figure 1

Overview of mC contexts in LDM of Hainan black goats and hybrid goats. A Density and (B) levels of mC within three contexts (CG, CHG, CHH) across the genome of LDM tissues in two goat species, with Hainan black goats at the center in A and hybrid goats in B, from innermost to outermost, represent gene density, hybrid mCHH, Hainan mCHH, hybrid mCHG, Hainan mCHG, hybrid mCG, Hainan CG. C Comparison of the average log2(mC number) in three mC contexts between Hainan black goats and hybrid goats, with the concentric circles indicating log2(mC number) values of 0, 4, 8, 12, and 16. D Mean mC levels within the three contexts (CG, CHG, and CHH) in various genomic functional regions (promoter, exon, intron, CGI, CGI shore, 3’UTR, 5’UTR, repeat) for LDM tissues in both goat species

The average mC levels in genomic functional regions (promoter, exon, intron, CGI, CGI shore, 3’UTR, 5’UTR, repeat) were consistent across the two species. The mCHG and mCHH displayed a little fluctuation in all functional regions of the LDM in Hainan black goats and hybrid goats, while mCG exhibited a gradual decline in the promoter region and relatively lower levels in the 5’UTR region (Fig. 1D and Fig. S6). Additionally, a higher methylation level of mCG was observed within the gene body (exon, intron) and downstream (3’UTR) regions in both goat species (Fig. 1D and Fig. S6).

The identification and annotation of DMR

The DMR length distributed from 50 bp to 5000 bp, and primarily focused on the 50 bp and 200 bp, the CG length was most often observed at 52 bp, while both CHG and CHH peaked near 53 bp or 200 bp (Fig. S7). The CG content in hybrid LDM tissue was distributed at different regions with diverse methylation levels on a minority of chromosomes (Fig. 2A), while the CHG and CHH presented several significant differences in methylation level over all chromosomes (Fig. S8). The most significant difference (p < 0.05) in hypermethylated (hyper) DMRs between the two goat species was a length of 507 bp of CG content in the CGI region on chromosome 13 (areaStat = 1291.31), followed by a length of 384 bp in the CHH repeat region on chromosome 1 (areaStat = 1288.92) (Fig. 2B and Table S3), suggesting the maximal difference of hyper DMRs was located within regulatory regions. However, the most significant difference in hypomethylated (hypo) DMRs was a 655 bp region of CHH content in intron on chromosome 19 (areaStat = -3134.54), demonstrating the hypo DMRs mainly regulated the gene body region (Fig. 2B and Table S3).

Fig. 2
figure 2

The DMR analysis between Hainan black goats and hybrid goats. A The outermost layer and innermost layers of the circle depict DNA methylation levels for the CG content in hybrid goats and Hainan black goats, respectively. The middle layer represents differences in DNA methylation levels between the two goat species. Each bin in the circos figure was 6,553,600bp. Dark blue bars indicate significantly lower (p < 0.05) DNA methylation levels in hybrid goats compared to Hainan black goats, while red bars signify significantly higher (p < 0.05) DNA methylation levels. B The outer and the inner layers distinguish between hyper and hypo DMR, while the two middle layers show the density of TE and genes. The small red, blue, and purple circles represent CG, CHG, and CHH contexts, respectively. The position of these circles relative to the center indicated the magnitude of the difference in methylation levels (diff.methylation) within DMRs. The size of the circles corresponds to the extent of the difference (areaStat) within DMRs. C Distribution of genes within DMRs across various genomic functional regions. D Comparison of gene functional regions between hyper and hypo DMRs

Moreover, the DMRs represented a total of 3,269 annotated genes, including 1,586 genes in hyper DMRs, while 2566 genes in hypo DMRs, and 883 genes appeared in both DMRs (Fig. S9). The hyper DMRs in the CG context contained more genes than that of hypo DMRs in the promoter and 5’UTR regions (Fig. 2C), with fewer in the exon, intron, and 3’UTR, indicating the key role of hyper DMRs on the upstream regulatory region, while hypo DMRs predominantly functioned in the gene body and downstream (Fig. 2C). Overall, the hyper DMRs accounted for 38.26% of the gene body, promoter, 3’UTR and 5’UTR regions, whereas hypo DMRs occupied the remaining 61.74% (Fig. 2D), suggesting the total number of hypo DMR genes in the gene body and regulatory regions was significantly higher than those of hyper DMRs. However, the areaStat of the most significantly different genes in hyper DMRs was significantly elevated compared to hypo DMRs (Fig. 2D).

The enriched pathways of DMRs related with muscle growth in LDM tissues

The genes of hyper DMRs in LDM tissue were enriched in the Rap1 signaling pathway and calcium signaling pathway (Fig. 3A). Moreover, the biological processes from GO analysis were separated into three categories: bone morphogenesis, synaptic signaling, and extracellular matrix organization (Fig. 3B), All of these impact myofiber formation and differentiation. Importantly, bone morphogenesis is closely tied to skeletal muscle growth as it took up the largest proportion of GO terms. Additionally, genes of hypo DMRs were primarily enriched in the mitogen-activated protein kinase (MAPK) signaling pathway, Phosphatidylinositol-3-kinase (PI3K-Akt) signaling pathway, calcium signaling pathway (Fig. 3C), and the molecular function of actin filament binding and GTPase activity (Fig. 3D), which played an important role in regulating myofiber activity and energy metabolism. Notably, Adenylate Cyclase 5 (ADCY5) and AKT Serine/Threonine Kinase 2 (AKT2) were involved in almost all the top pathways of hyper DMRs (Tables S4 and S5). AKT2 was considered an essential kinase for most vital activities, including cellular proliferation and differentiation, translation, and protein synthesis. Thus, the AKT2 would cause nutritional deficiency, inflammation, and immune dysfunction across tissues.

Fig. 3
figure 3

Functional analysis of the gene sets within DMRs and DEGs. The A KEGG and B GO enrichment of hyper DMR genes. C KEGG and D GO enrichment of hypo DMR genes. E PPI network of genes intersecting between DMRs and DEGs. The 11 hub genes are highlighted by a red outline

Furthermore, the DEGs intersected with the genes of hyper and hypo DMRs, and a gene set of 100 DEGs was obtained within hyper DMRs and 89 DEGs within hypo DMRs, respectively, resulting in 166 unique DMR genes (Table S6). Notably, the PPI network constructed by the 166 genes identified 11 proteins with strong correlation with others (Fig. 3E and Table S7), including Protein Kinase CGMP-Dependent 1 (PRKG1), AKT2, Forkhead Box O3 (FOXO3), Forkhead Box O6 (FOXO6), G Protein Subunit Alpha O1 (GNAO1), ADCY5, Enhancer of Zeste 2 Polycomb Repressive Complex 2 Subunit (EZH2), Zinc Finger And BTB Domain Containing 16 (ZBTB16), Nuclear Receptor Subfamily 4 Group A Member 1 (NR4A1), Insulin Receptor Substrate 2 (IRS2) and RUNX Family Transcription Factor 1 (RUNX1). The 11 proteins regulated diverse functions, including nutrient metabolism, growth hormone synthesis, inflammation, and cellular immune activity. For instance, ADCY5, AKT2, and IRS2 are important proteins in muscle growth hormone synthesis, lipolysis, and glycogen storage, while GNAO1 is involved in cholinergic synapse regulation, as well as FOXO3 and FOXO6, which regulates the balance of apoptosis and autophagy.

The selection of DMR genes associated with muscle growth

Furthermore, the gene sets intersecting with hyper DMRs and DEGs was enriched in nine pathways, focused on the regulation of lipolysis in adipocytes, cholinergic synapse, growth hormone synthesis, secretion and action (Fig. 4A), indicating the predominant function of energy metabolism within the LDM tissue. When comparing the pathways and the PPI network of hyper DMRs and DEGs, the critical genes that appeared across almost all the pathways also appeared in the PPI hub region, suggesting the consistency between the PPI and enriched pathways. The 11 hub genes were clearly separated into two clusters representing Hainan black goats and hybrid goats (Fig. 4B), with the methylation sites being dispersed across the gene body and regulatory regions, which contained all three contexts of CG, CHG and CHH (Fig. 4C, Table S8). Although it was clear that hypermethylation occurring in gene promoter regions would inhibit gene expression [27], the majority of the expression of 100 hyper DMR genes from the intersection DMRs and DEGs were not linked to their methylation level, nor related to gene promoter methylation (Fig. S10). However, 11 genes were identified with their expression significantly correlated with gene methylation levels in hypo DMRs (Fig. 4D, Table S9), and six genes (AKT2, FOXO3, FOXO6, PRKG1, ZBTB16, NR4A1) were also hub genes, functioning in myofiber growth, inflammation, and endocytosis.

Fig. 4
figure 4

Functional analysis of gene sets intersecting between DMRs and DEGs. A KEGG enrichment of the intersection gene sets. B Expression patterns of the 11 hub genes from Hainan black goats and hybrid goats. C A depiction of the connections among genes, functional regions, and C contents of hub genes. D Correlations between differential gene expression and the degree of differential methylation levels within hypo DMRs

Correlations between gene expression and growth traits in the LEA

Three growth traits, including the area, height, and weight of the loin eye muscle area (LEA), were observed in both hybrid goats and Hainan black goats (Table S10). Correlations were established between genes expression and growth traits. PRKG1 was the only gene exhibiting strongly negative correlations with all three growth traits (Fig. 5A), and NR4A1 displayed the strongest negative correlation with two growth traits, namely height and weight (R = -0.91, P = 0.01; R = -0.89, P = 0.01) (Fig. 5B). Other genes, including FOXO3, FOXO6, and GNAO1 were positively correlated with two growth traits of LEA (Fig. 5C-E). Specifically, FOXO3 and FOXO6 were both positively correlated with area and weight. Additionally, ADCY5 was positively linked to weight (Fig. 5F), and participated in insulin-glucose signaling for glucose transferase [28]. AKT2 gene accounted for the largest AKT isoform in human skeletal muscle [29] and was positively correlated with the LEA (Fig. 5G). Transcription factor ZBTB16 and the target gene of transcription factor EZH2 were linked to weight and LEA, respectively (Fig. 5H and I). Aside from two transcription factor-associated genes, the remaining genes in the PPI node scoring over 5 were ADCY5, AKT2, PRKG1 and FOXO3, and the homologous genes of FOXO3 and FOXO6 frequently cooperated in glucose metabolism. Therefore, the five genes were considered to be the most critical genes for muscle fiber growth.

Fig. 5
figure 5

Correlations of hub genes with growth traits of LEA. APRKG1 was negatively correlated with the weight, height, and LEA. BNR4A1 was negatively correlated with the height and weight. C-EFOXO3, FOXO6 and GNAO1 were positively correlated with two of the three growth traits. F-IADCY5, AKT2, ZBTB16 were positively correlated with one trait, and EZH2 was negatively correlated with one trait. All hollow circles in A-E indicate Hainan black goats, and solid circles indicate hybrid goats

The cooperation of critical genes in LDM

The muscle fibers were surrounded by numerous blood vessels and macrophages, constructing a microenvironment that ensured normal development and movement of myofibers. According to the enrichment pathways of DMR and DEGs, the protection of vascular endothelial cells may be critical for vasodilatation and the growth of LDM. The ADCY family, including ADCY1, ADCY5, ADCY8, and ADCY9, all took part in insulin secretion and growth hormone synthesis, secretion and action pathway, and vascular smooth muscle contraction upon DMR enrichment (Tables S4 and S5), especially ADCY5 not only participated in these functions, also attributed to the regulation of glucose transport. AKT2 was a hub gene linked with diverse signaling pathways, including phagocytosis and cellar growth (Tables S4 and S5). PRKG1 and FOXO3 were markers for the balance of autophagy of vascular endothelial cells [30, 31], while PRKG1 also influenced the smooth vessel activation. Moreover, FOXO6 maintained myotube proliferation and differentiation, and protected against atrophy [32]. Therefore, the healthy microenvironment for muscle growth was largely hinged on vascular activation and nutrient transport (Fig. 6).

Fig. 6
figure 6

Identification and cooperation of hub genes associated with muscle growth. The left side depicts muscle fibers densely interspersed with blood vessels. The right side illustrates a microenvironment encompassing muscle fibers, blood vessels, and phagocytes. The central pink circle, filled with small myofibrils, represents a muscle fiber, surrounded by red and blue blood vessels, and several genes are interconnected around the muscle fiber, jointly contributing to glycogen storage. Red arrows indicate genes with increased expression, while blue arrows indicate genes with decreased expression

The five most important genes, including ADCY5, AKT2, FOXO3, FOXO6, and PRKG1, were identified by RT-qPCR, and top four genes were significantly upregulated (p < 0.05) in hybrid goats compared to Hainan goats, while PRKG1 was downregulated (p < 0.05) (Fig. S11). Therefore, the increasing expression of ADCY5, AKT2, FOXO3 and FOXO6 may be responsible for promoting muscle fiber growth. In addition, the decreasing expression of PRKG1 suppressed vasoconstriction to release glucose more easily from the blood.


DNA methylation is an epigenetic regulation approach that is important for gene expression, tissue development [33], inflammation [34], autophagy and apoptosis [35]. The purpose of this study was to investigate the effects of DNA methylation on the LDM growth-related genes. Although the genome-wide methylation patterns of LDM between Hainan black goats and hybrid goats were generally similar (Fig. 1A and B), the gene number and functions were massively different during the DMR. The number of methylation sites in the gene body within hypo DMR was significantly higher than that within hyper DMR (Fig. 2C). In contrast, genes in hyper DMR also acted as hypermethylated genes with retrocopies available for glycan and lipid synthesis [36, 37]. The transcription factor EZH2, triggered DNA methylation in most cell activations [38], which were significantly downregulated expression in hybrid goats (p < 0.05), consisting of a negative correlation with LEA (Fig. 5I).

The DMR gene functions were focused on nutrient metabolism, growth hormone synthesis, vasodilatation, and endocytosis. Skeletal muscle is highly vascularized, and the muscle stem cell proliferation was accompanied by remarkable vascular and neuronal regeneration and differentiation [38, 39]. Vasodilatation could promote the skeletal muscle cells to uptake glucose [40]. Moreover, the inhibition of overexpressed PRKG1 could limit vascular sclerosis to reverse aortopathy in Marfan syndrome mice [31], suggesting that PRKG1 might be important in blood vessel activity. In this study, the expression of PRKG1 was negatively correlated with the height and weight of LEA (Fig. 5A), which might be due to its decreasing expression benefitting vasodilatation to facilitate nutrient release, and acting as a switch to turn on glucose metabolism.

Similarly to previous studies, PI3K-Akt and MAPK signaling pathways activation could stimulate LDM development and differentiation in sheep and goats [41, 42], and these pathways were enriched in goat hypo DMR in this study. Notably, FOXO3, FOXO6 and AKT2 genes were related to PI3K-AKT signaling pathway (Table S5), which was important for skeletal muscle protein synthesis [43]. FOXO gene family is evolutionarily conserved, and FOXO3 has been linked to longevity due to its vascular protection and regeneration [44]. Moreover, FOXO3 acted as a cell surveillance mechanism to balance the apoptosis and autophagy [30], and their imbalance caused dystrophic features of skeletal muscle [45]. Furthermore, FOXO3 eliminated senescent cells and harmful substances through hypo methylation [46]. Similarly, the homologous gene FOXO6 maintained myotubes and protected against skeletal muscle atrophy [32]. In this study, both FOXO3 and FOXO6 were positively correlated with the area and weight of loin eye muscle, and both were located in hypo DMR. Therefore, it was clear that these two genes promoted myofiber development of skeletal muscle. In contrast, the methylation caused by maternal undernutrition treatment had no significant effects on the methylation of PI3K and AKT genes in goat semitendinosus (ST) and vastus lateralis (VL) muscles, nor did influence most of the gene expression levels of myogenic factors (MYFs) [47].

Besides the important genes in the top pathways, several genes were correlated to numerous functions and associated with the growth traits of muscle. For example, NR4A1 was a glucose indicator as its expression was related to the insulin sensitivity of skeletal muscle [48], and its reduced methylation level influenced gene expression in adult sheep LDM [43]. Moreover, IRS2 exhibited diverse functional features of insulin-dependent regulation of glucose and lipid metabolism in skeletal muscle and liver. Therefore, reduced expression of NR4A1 and IRS2 could inhibit sugar usage, whereas promote glycogen storage in LDM. Contrary to sheep containing opposite methylation levels and gene expression [43], these genes were downregulated (p < 0.05) in hybrid goats LDM compared to Hainan goats, indicating the different effects on gene expression caused by methylation in goats and sheep. Similarly, ZBTB16 expression was accompanied by insulin-stimulated glucose transitioning into skeletal muscle glycogen, it acted as a surveillant for glucose quantity in muscles [49], therefore, the upregulated (p < 0.05) ZBTB16 was response for muscular energy provision. EZH2 was critical for the maintenance of smooth muscle physical form and function, and the inhibition of EZH2 expression through promoter methylation could improve smooth muscle cell development [50]. However, its function in goat muscle was not reported. In this study, the downregulated (p < 0.05) EZH2 gene may enhance the transport function of vascular smooth muscle in hybrid goats LDM. The missense variants of GNAO1 present with paroxysmal dystonic manifestations, thus, this gene plays a role in muscular tension [51], and its increasing expression might enhance the LDM tension of hybrid goats. RUNX1 contributed to myofiber development and differentiation, and induced muscle regeneration by promoting satellite cell proliferation [52], indicating the increased expression of RUNX1 improved hybrid goats muscular growth.

Importantly, glucose transport and glycogen synthesis preferred to select ADCY5-mediated insulin-glucose pattern in LDM [53]. Glycogen storage was essential for myofiber movement [54], and the deletion of ADCY5 was related to abnormal muscle activity during active movement as neurotransmitter disorders caused by glucose transporter deficiency [55]. In this study, upregulated (p < 0.05) ADCY5 expression was positively correlated with the weight of LEA (Fig. 5F), which might be due to its linkage to glucose transport. The AKT2 gene accounts for the largest percentage of AKT isoform in skeletal muscle to regulate myoblast differentiation and control the morphology of the muscle fibers [29], which may explain its positive correlation with the LEA. Further, other genes related to myofiber proliferation and differentiation also positively correlated with LEA and contributed to goat LDM growth.


This study investigated the DNA methylation differences between the LDM of Hainan black goats and hybrid goats, acquired the critical genes, and analyzed their functions on myofiber growth. Notably, nine genes, including PRKG1, AKT2, FOXO3, FOXO6, GNAO1, ADCY5, EZH2, ZBTB16, and NR4A1, were significantly correlated (p < 0.05) with three growth traits of LDM partially or totally, and in particular, PRKG1 was negatively related to all growth traits. The five key genes identified in DMR connected nutrient transport, vasoconstriction, and inflammation. Specifically, FOXO3 and FOXO6 monitored the balance of apoptosis and autophagy to maintain the immune microenvironment; PRKG1 regulated the vasoconstriction to release nutrients, and its decreased expression promoted the glucose transfer and glycogen metabolism mediated by ADCY5 and AKT2 genes. Therefore, this study highlights the importance of methylation on the glycogen synthesis and myofiber growth-related genes in goat LDM, as well as provides the basic molecular mechanism for LDM development.


Tissue sampling

Given that the weight of hybrid goats (Nubian goats × Hainan black goats) was significantly higher than that of Hainan black goats (p < 0.05), we investigated the molecular mechanism ofzLDM growth, and examined the correlation between critical genes and growth traits of weight and LEA. The 12 individuals of 7-months-old male goats were purchased from Hainan Chengmai Fuyang Animal Husbandry Company (Hainan, China), and theses goats were maintained in a similar environment of free food and water for an adaptive week to reduce stress before slaughtering. There were 6 individuals for each group of Hainan black goats and hybrid goats, after the transitional period, all the individuals were experienced 24h of feed withdrawal, and three individuals for each group were randomly selected and transferred to the slaughterhouse of the HAAS YongFa goat experimental facility. The goats were stunned by electric shock, and exsanguinated, flushed and split at unconsciousness. All methods are reported in accordance with ARRIVE guidelines ( for the reporting of animal experiments. The LDM were collected from each goat after slaughtering, and stored in liquid nitrogen for further analysis.

Whole genome bisulfite sequencing and quality control

The DNA was extracted from LDM tissues by DNA extraction kit (Qiagen, 69504), and the concentration was detected by Nanodrop and Qubit 4.0, and then the high-quality DNA was sent to Novogene Company and subjected to WGBS, which was treated by bisulfite to separate methylated and unmethylated base C firstly, and then performed library construction of PE150 [56]. The inserted size of library was verified using an Agilent DNA-1000 Kit on an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA), and each library was sequencing on the Illumina NovaSeq 6000 platform. The raw sequenced data was evaluated by FastQC, and trimmed by fastp [57].

Methylation stat analysis

The clean data was mapped to goat genome reference (ARS1, ASM170441v1) by Bismark [58], which transformed both the sequenced reads and genomic data from C to T or from G to A, separately, and four transforms were compared in pairs and removed duplications, and the mapping stats, containing unique mapping rate, sequencing depth and coverage, and the coverage of C, were evaluated to determine the high-quality data.

The CpG, CHG and CHH are three sequence contexts of methylated C (mC), where H is indicated as A, T, or C. The unique high quality bisulfite reads are classified into 6 types, which are methylated/unmethylated C in CpG, CHG, and CHH, separately. The accuracy of methylation detection is evaluated by binomial distribution [59], and the qualified methylation is required by sequencing depth more than 5 and q-value≤0.01. The methylation stats were summarized based on the genomic functional regions (promoter, exon, intron, CpG island (CGI)), where promoter was defined as the region of 2Kb upstream the gene transcription start site (TSS), and CGI and CGI shore were predicted by the cpgIslandExt (src/utils/cpgIslandExt/). DMR resulted from hybrid goats compared with Hainan black goats were identified by DSS package [60], taking account of main three factors, that are spatial correlation, sequencing depth of the sites, and the variance among biological replicates. The circos figures displayed methylation level, density and DMR were draw by TBtools [61].

Expression of methylated genes

The RNA-seq data for each group were replicated for 3 individuals. Raw reads of RNA-seq were filtered by Trimmomatic (Version 0.36.6) [62] to removed adapters and low-quality reads, and clean reads were mapped to the goat genome reference (ARS1, ASM170441v1) using HISAT2 with default parameters (Version 2.2.1) [63], and read counts were calculated by HTSeq v. 0.9.1 [64]. Differentially expressed genes (DEGs) were identified by DEseq2 [65], requiring |log2(foldchange) |≥1, and adjusted by p ≤ 0.05. The primers sequences for RT-qPCR were designed and listed (Table S11).

Enrichment analysis of methylation

The structure and gene annotation of DMR were further enriched for functions. The pathway analysis and functional classification were conducted by R package clusterProfiler [66] based on KEGG [67] and GO database, and the results were statistically adjusted by p ≤ 0.05. To simplify enriched GO terms, GOSemSim was used to calculate the similarity of GO terms and remove those redundancy terms, and only keep representative GO terms [68]. A protein-protein interactions (PPIs) network was constructed using the intersecting genes from DMRs and DEGs by STRING ( Subsequently, the visualization of enrichment results was implemented by R software.

Statistical analysis

Both the correlation coefficient (R) either between methylation level and the gene expression, or between the gene expression and LDM traits, were calculated by R packages and detected by statistical method of pearson test, and the results were statistically adjusted by p≤0.05.

The statistical analysis for RT-qPCR was using Student's t test method within two-sided for 3 replicated individuals, and adjusted by p ≤ 0.05.

Availability of data and materials

The raw data that support the findings of this study are available from the accession number PRJNA924997 in the NCBI BioProject database (



longissimus dorsi muscle


whole-genome bisulfite sequencing


differentially methylated regions


differentially expressed genes


mitogen-activated protein kinase




adenylate cyclase 5

AKT2 :

AKT serine/threonine Kinase 2


protein kinase CGMP-dependent 1


forkhead box O3


forkhead box O6


G protein subunit alpha O1

EZH2 :

enhancer Of zeste 2 polycomb repressive complex 2 subunit

ZBTB16 :

zinc finger and BTB domain containing 16

NR4A1 :

nuclear receptor subfamily 4 group A member 1

IRS2 :

insulin receptor substrate 2


RUNX family transcription factor 1


loin eye muscle area


methylated C


CpG island


transcription start site


protein-protein interactions


  1. Daly KG, Mattiangeli V, Hare AJ, Davoudi H, Fathi H, Doost SB, et al. Herded and hunted goat genomes from the dawn of domestication in the Zagros Mountains. PNAS. 2021;118(25):e2100901118.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Migdal W, Kawecka A, Sikora J, Migdal L. Meat quality of the native Carpathian goat breed in comparison with the Saanen breed. Animals (Basel). 2021;11(8):2220.

    Article  PubMed  PubMed Central  Google Scholar 

  3. Gonzalez-Rivas PA, Chauhan SS, Ha M, Fegan N, Dunshea FR, Warner RD. Effects of heat stress on animal physiology, metabolism, and meat quality: A review. Meat Sci. 2020;162:108025.

    Article  CAS  PubMed  Google Scholar 

  4. Luo R, Dai X, Zhang L, Li G, Zheng Z. Genome-Wide DNA Methylation Patterns of Muscle and Tail-Fat in DairyMeade Sheep and Mongolian Sheep. Animals (Basel). 2022;12(11):1399.

    Article  PubMed  PubMed Central  Google Scholar 

  5. Iwayama K, Tanabe Y, Tanji F, Ohnishi T, Takahashi H. Diurnal variations in muscle and liver glycogen differ depending on the timing of exercise. J Physiol Sci. 2021;71(1):35.

    Article  CAS  PubMed  Google Scholar 

  6. de Moraes WMAM, de Almeida FN, Dos Santos LEA, Cavalcante KDG, Santos HO, Navalta JW, et al. Carbohydrate Loading Practice in Bodybuilders: Effects on Muscle Thickness, Photo Silhouette Scores, Mood States and Gastrointestinal Symptoms. J Sports Sci Med. 2019;18(4):772–9 (PMID: 31827362).

    PubMed  PubMed Central  Google Scholar 

  7. Biffin TE, Hopkins DL, Bush RD, Hall E, Smith MA. The effects of season and post-transport rest on alpaca (Vicunga pacos) meat quality. Meat Sci. 2020;159:107935.

    Article  CAS  PubMed  Google Scholar 

  8. Park H, Seo KS, Lee M, Seo S. Identification of meat quality-related differentially methylated regions in the DNA of the longissimus dorsi muscle in pig. Anim Biotechnol. 2020;31(3):189–94.

    Article  CAS  PubMed  Google Scholar 

  9. Kurdyukov S, Bullock M. DNA methylation analysis: choosing the right method. Biology (Basel). 2016;5(1):3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Wang K, Wu P, Wang S, Ji X, Chen D, Jiang A, et al. Genome-wide DNA methylation analysis in Chinese Chenghua and Yorkshire pigs. BMC Genom Data. 2021;22(1):21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Barazandeh A, Mohammadabadi M, Ghaderi-Zefrehei M, Rafeie F, Imumorin IG. Whole genome comparative analysis of CpG islands in camelid and other mammalian genomes. Mammalian Biology. 2019;98:73–9.

    Article  Google Scholar 

  12. Davenport KM, Massa AT, Bhattarai S, McKay SD, Mousel MR, Herndon MK, et al. Characterizing genetic regulatory elements in ovine tissues. Front Genet. 2021;12:628849.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Shen L, Du J, Xia Y, Tan, Tan Z, Fu Y, Yang Q, et al. Genome-wide landscape of DNA methylomes and their relationship with mRNA and miRNA transcriptomes in oxidative and glycolytic skeletal muscles. Sci Rep. 2016;6:32186.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Denoyelle L, de Villemereuil P, Boyer F, Khelifi M, Gaffet C, Alberto F, et al. Genetic variations and differential DNA methylation to face contrasted climates in small ruminants: an analysis on traditionally-managed sheep and goats. Front Genet. 2021;12:745284.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Yuan Z, Ge L, Zhang W, Lv X, Wang S, Cao X, et al. Preliminary results about lamb meat tenderness based on the study of novel isoforms and alternative splicing regulation pathways using Iso-seq, RNA-seq and CTCF ChIP-seq data. Foods. 2022;11(8):1068.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Guo X, Puttabyatappa M, Domino SE, Padmanabhan V. Developmental programming: Prenatal testosterone-induced changes in epigenetic modulators and gene expression in metabolic tissues of female sheep. Mol Cell Endocrinol. 2020;514:110913.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Rahmatalla SA, Arends D, Reissmann M, Said Ahmed A, Wimmers K, Reyer H, et al. Whole genome population genetics analysis of Sudanese goats identifies regions harboring genes associated with major traits. BMC Genet. 2017;18(1):92.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Osman OA, Elkhair NM, Abdoun KA. Effects of dietary supplementation with different concentration of molasses on growth performance, blood metabolites and rumen fermentation indices of Nubian goats. BMC Vet Res. 2020;16(1):411.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Cao Y, Xu H, Li R, Gao S, Chen N, Luo J, et al. Genetic basis of phenotypic differences between Chinese Yunling black goats and Nubian goats revealed by allele-specific expression in their F1 hybrids. Front Genet. 2019;10:145.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Jiang S, Huo D, You Z, Peng Q, Ma C, Chang H, et al. The distal intestinal microbiome of hybrids of Hainan black goats and Saanen goats. PLoS One. 2020;15(1):e0228496.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Hua R, Zhou L, Zhang H, Yang H, Peng W, Wu K. Studying the variations in differently expressed serum proteins of Hainan black goat during the breeding cycle using isobaric tags for relative and absolute quantitation (iTRAQ) technology. J Reprod Dev. 2019;65(5):413–21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. 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;4:6546.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Ponsuksili S, Trakooljul N, Basavaraj S, Hadlich F, Murani E, Wimmers K. Epigenome-wide skeletal muscle DNA methylation profiles at the background of distinct metabolic types and ryanodine receptor variation in pigs. BMC Genomics. 2019;20(1):492.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Namous H, Peñagaricano F, Del Corvo M, Capra E, Thomas DL, Stella A, et al. Integrative analysis of methylomic and transcriptomic data in fetal sheep muscle tissues in response to maternal diet during pregnancy. BMC Genomics. 2018;19(1):123.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Deng M, Zhang G, Cai Y, Liu Z, Zhang Y, Meng F, Wang F, Wan Y. DNA methylation dynamics during zygotic genome activation in goat. Theriogenology. 2020;156:144–54.

    Article  CAS  PubMed  Google Scholar 

  26. Kang B, Wang J, Zhang H, Shen W, El-Mahdy Othman O, Zhao Y, Min L. Genome-wide profile in DNA methylation in goat ovaries of two different litter size populations. J Anim Physiol Anim Nutr (Berl). 2022;106(2):239–49.

    Article  CAS  PubMed  Google Scholar 

  27. Hughes AL, Kelley JR, Klose RJ. Understanding the interplay between CpG island-associated gene promoters and H3K4 methylation. Biochim Biophys Acta Gene Regul Mech. 2020;1863(8):194567.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. de Gusmao CM, Silveira-Moriyama L. Paroxysmal movement disorders - practical update on diagnosis and management. Expert Rev Neurother. 2019;19(9):807–22.

    Article  CAS  PubMed  Google Scholar 

  29. Matheny RW Jr, Geddis AV, Abdalla MN, Leandry LA, Ford M, McClung HL, et al. AKT2 is the predominant AKT isoform expressed in human skeletal muscle. Physiol Rep. 2018;6(6):e13652.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Fitzwalter BE, Thorburn A. FOXO3 links autophagy to apoptosis. Autophagy. 2018;14(8):1467–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Cardona A, Pagani L, Antao T, Lawson DJ, Eichstaedt CA, Yngvadottir B, et al. Genome-wide analysis of cold adaptation in indigenous Siberian populations. PLoS One. 2014;9(5):e98076.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Zhang L, Zhang Y, Zhou M, Wang S, Li T, Hu Z, et al. Role and mechanism underlying FOXO6 in skeletal muscle in vitro and in vivo. Int J Mol Med. 2021;48(1):143.

    Article  CAS  PubMed  Google Scholar 

  33. Moore LD, Le T, Fan G. DNA methylation and its basic function. Neuropsychopharmacology. 2013;38(1):23–38.

    Article  CAS  PubMed  Google Scholar 

  34. Ciechomska M, Roszkowski L, Maslinski W. DNA methylation as a future therapeutic and diagnostic target in rheumatoid arthritis. Cells. 2019;8(9):953.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Jeon M, Park J, Yang E, Baek HJ, Kim H. Regulation of autophagy by protein methylation and acetylation in cancer. J Cell Physiol. 2022;237(1):13–28.

    Article  CAS  PubMed  Google Scholar 

  36. Fang C, Zou C, Fu Y, Li J, Li Y, Ma Y, et al. DNA methylation changes and evolution of RNA-based duplication in Sus scrofa: based on a two-step strategy. Epigenomics. 2018;10(2):199–218.

    Article  CAS  PubMed  Google Scholar 

  37. Emran AA, Chatterjee A, Rodger EJ, Tiffen JC, Gallagher SJ, Eccles MR, et al. Targeting DNA methylation and EZH2 activity to overcome melanoma resistance to immunotherapy. Trends Immunol. 2019;40(4):328–44.

    Article  CAS  PubMed  Google Scholar 

  38. Almada AE, Wagers AJ. Molecular circuitry of stem cell fate in skeletal muscle regeneration, ageing and disease. Nat Rev Mol Cell Biol. 2016;17(5):267–79.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Ottone C, Krusche B, Whitby A, Clements M, Quadrato G, Pitulescu ME, et al. Direct cell-cell contact with the vascular niche maintains quiescent neural stem cells. Nat Cell Biol. 2014;16(11):1045–56.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Walsh LK, Ghiarone T, Olver TD, Medina-Hernandez A, Edwards JC, Thorne PK, et al. Increased endothelial shear stress improves insulin-stimulated vasodilatation in skeletal muscle. J Physiol. 2019;597(1):57–69.

    Article  CAS  PubMed  Google Scholar 

  41. Fan Y, Liang Y, Deng K, Zhang Z, Zhang G, Zhang Y, et al. Analysis of DNA methylation profiles during sheep skeletal muscle development using whole-genome bisulfite sequencing. BMC Genomics. 2020;21(1):327.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Zou J, Shen Y, Zou J, Yu J, Jiang Y, Huang Y, et al. Transcriptome-wide study revealed that n6-methyladenosine participates in regulation meat production in goats. Foods. 2023;12(6):1159.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Zhang ZK, Li J, Liu J, Guo B, Leung A, Zhang G, et al. Icaritin requires Phosphatidylinositol 3 kinase (PI3K)/Akt signaling to counteract skeletal muscle atrophy following mechanical unloading. Sci Rep. 2016;6:20300.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Yan P, Li Q, Wang L, Lu P, Suzuki K, Liu Z, et al. FOXO3-engineered human ESC-derived vascular cells promote vascular protection and regeneration. Cell Stem Cell. 2019;24(3):447-461.e8.

    Article  CAS  PubMed  Google Scholar 

  45. Yazid MD, Hung-Chih C. Perturbation of PI3K/Akt signaling affected autophagy modulation in dystrophin-deficient myoblasts. Cell Commun Signal. 2021;19(1):105.

    Article  PubMed  PubMed Central  Google Scholar 

  46. Choi S, Jeong HJ, Kim H, Choi D, Cho SC, Seong JK, et al. Skeletal muscle-specific Prmt1 deletion causes muscle atrophy via deregulation of the PRMT6-FOXO3 axis. Autophagy. 2019;15(6):1069–81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Zhou X, Yan Q, Liu L, Chen G, Tang S, He Z, et al. Maternal undernutrition alters the skeletal muscle development and methylation of myogenic factors in goat offspring. Anim Biosci. 2022;35(6):847–57.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Kasch J, Kanzleiter I, Saussenthaler S, Schürmann A, Keijer J, van Schothorst E, et al. Insulin sensitivity linked skeletal muscle NR4A1 DNA methylation is programmed by the maternal diet and modulated by voluntary exercise in mice. J Nutr Biochem. 2018;57:86–92.

    Article  CAS  PubMed  Google Scholar 

  49. Krupkova M, Liska F, Kazdova L, Šedová L, Kábelová A, Křenová D, et al. Single-gene congenic strain reveals the effect of ZBTB16 on dexamethasone-induced insulin resistance. Front Endocrinol (Lausanne). 2018;9:185.

    Article  PubMed  PubMed Central  Google Scholar 

  50. Maleszewska M, Gjaltema RA, Krenning G, Harmsen MC. Enhancer of zeste homolog-2 (EZH2) methyltransferase regulates transgelin/smooth muscle-22α expression in endothelial cells in response to interleukin-1β and transforming growth factor-β2. Cell Signal. 2015;27(8):1589–96.

    Article  CAS  PubMed  Google Scholar 

  51. Dzinovic I, Škorvánek M, Necpál J, Boesch S, Švantnerová J, Wagner M, et al. Dystonia as a prominent presenting feature in developmental and epileptic encephalopathies: A case series. Parkinsonism Relat Disord. 2021;90:73–8.

    Article  PubMed  Google Scholar 

  52. Bao M, Liu S, Yu XY, Wu C, Chen Q, Ding H, et al. Runx1 promotes satellite cell proliferation during ischemia - Induced muscle regeneration. Biochem Biophys Res Commun. 2018;503(4):2993–7.

    Article  CAS  PubMed  Google Scholar 

  53. Huang LO, Rauch A, Mazzaferro E, Preuss M, Carobbio S, Bayrak CS, et al. Genome-wide discovery of genetic loci that uncouple excess adiposity from its comorbidities. Nat Metab. 2021;3(2):228–43.

    Article  CAS  PubMed  Google Scholar 

  54. Migocka-Patrzałek M, Elias M. Muscle glycogen phosphorylase and its functional partners in health and disease. Cells. 2021;10(4):883.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Pérez-Dueñas B, Gorman K, Marcé-Grau A, Ortigoza-Escobar JD, Macaya A, Danti FR, et al. The genetic landscape of complex childhood-onset hyperkinetic movement disorders. Mov Disord. 2022;37(11):2197–209.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Miura F, Enomoto Y, Dairiki R, Ito T. Amplification-free whole-genome bisulfite sequencing by post-bisulfite adaptor tagging. Nucleic Acids Res. 2012;40(17): e136.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):i884–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics. 2011;27(11):1571–2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Habibi E, Brinkman Arie B, Arand J, Kroeze LI, Kerstens HH, Matarese F, et al. Whole-genome bisulfite sequencing of two distinct interconvertible DNA methylomes of mouse embryonic stem cells. Cell Stem Cell. 2013;13(3):360–9.

    Article  CAS  PubMed  Google Scholar 

  60. Wu H, Xu T, Feng H, Chen L, Li B, Yao B, et al. Detection of differentially methylated regions from whole-genome bisulfite sequencing data without replicates. Nucleic Acids Res. 2015;43(21):e141.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Chen C, Chen H, Zhang Y, Thomas HR, Frank MH, He Y, et al. TBtools: an integrative toolkit developed for interactive analyses of big biological data. Mol Plant. 2020;13(8):1194–202.

    Article  CAS  PubMed  Google Scholar 

  62. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Anders S, Pyl PT, Huber W. HTSeq–a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31(2):166–9.

    Article  CAS  PubMed  Google Scholar 

  65. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Kanehisa M, Goto S. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 2000;28:27–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Yu G, Li F, Qin Y, Bo X, Wu Y, Wang S. GOSemSim: an R package for measuring semantic similarity among GO terms and gene products. Bioinformatics. 2010;26(7):976–8.

    Article  CAS  PubMed  Google Scholar 

Download references


Not applicable.


This research was funded by the Hainan Provincial Natural Science Foundation of China (320CXTD445); and the National Key Research and Development Program of China (2021YFD1200802); and the Special Fund for the Development of Local Science and Technology Guided by Central Government (ZY2022HN09). The funding bodies played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Author information

Authors and Affiliations



YR and ZC conceived and designed the study. YR and XC performed the bioinformatics. YR performed statistical analyses and prepared figures. FW, XZ, HL, YL, LH and XH contributed to the material preparation. YR and ZC wrote original draft. YR, ZC, RS, LW and YZ reviewed the manuscript. All authors have read and approved the final version of the manuscript.

Corresponding author

Correspondence to Zhe Chao.

Ethics declarations

Ethics approval and consent to participate

The animal used in this study was Capra hircus, that were Hainan black goats and its F1 hybrid offspring crossed with Nubian goats. Capra hircus is an important economic goat in China, which can be bought from the Chinese market. All the purchased animals were reared on standard water and food condition. All procedures were conducted according to the Regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China, revised in March 2017), and approved by the Institutional Animal Care and Use Committee at the Experimental Animal Center of Hainan Academy of Agricultural Science (HNXMSY-20210503).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1. 

The quality control of sequencing data. Table S2. The mapping rate of clean data of LDM tissues. Table S3. The distribution and the conditions of the top DMR between two goat species. Table S4. The enriched top pathways of hyper DMR. Table S5. The enriched top pathways of hypo DMR. Table S6. The intersection set of DMR genes and DEGs. Table S7. The PPI node score of 11 hub genes. Table S8. The methylation sites of 11 hub genes. Table S9. The expressions and methylation levels of 11 genes with correlations in hypo DMRs. Table S10. The growth traits of hybrid goats and Hainan black goats in LEA. Table S11. The primer sequences for RT-qPCR. Fig. S1. The sequencing depth distribution of all bases for Hainan black goats and hybrid goats. Fig. S2. The sequencing depth distribution of base accumulative in Hainan black goats and hybrid goats. The x-axis meant the sequencing depth; the y-axis meant the percentage of accumulative fraction of bases. Fig. S3. The average CG context number of chromosomes in Hainan black goats and hybrid goats. Blue color meant Hainan black goats, and orange color meant Hybrid goats. Fig. S4. The classification of mean mC proportions in LDM. Fig. S5. The methylation level distribution of mC contexts. Fig. S6. The average mC levels in different functional regions. Fig. S7. The length density distribution of three mC contexts in LDM. Fig. S8. The DMR analysis of LDM between Hainan black goats and hybrid goats. Fig. S9. The number of annotated genes in DMRs. Fig. S10. The correlations between methylation levels and the expression of DEGs in hyper DMR. Fig. S11. The RT-qPCR results of 5 most important genes.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Ren, Y., Chen, X., Zheng, X. et al. Diverse WGBS profiles of longissimus dorsi muscle in Hainan black goats and hybrid goats. BMC Genom Data 24, 77 (2023).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: