Integrative ATAC-seq and RNA-seq Analysis of the Longissimus Muscle of Luchuan and Duroc Pigs

Luchuan pig is a typical obese pig breed in China, and the diameter and area of its longissimus dorsi muscle fibers are significantly smaller than those of Duroc (lean) pig. Skeletal muscle fiber characteristics are related to meat quality of livestock. There is a significant correlation between the quality of different breeds of pork and the characteristics of muscle fiber, which is an important factor affecting the quality of pork. The diameter and area of muscle fibers are related to muscle growth and development. Therefore, we used the assay for transposase-accessible chromatin using sequencing (ATAC-seq) and RNA sequencing (RNA-seq) analysis to investigate the potential mechanism underlying the difference in skeletal muscle growth and development between the two types of pigs. First, transposase-accessible chromatin was analyzed to map the landscape of open chromatin regions and transcription factor binding sites. We identified several transcription factors that potentially affected muscle growth and development, including TFAP4, MAX, NHLH1, FRX5, and TGIF1. We also found that transcription factors with basic helix-loop-helix structures had a preference for binding to genes involved in muscle development. Then, by integrating ATAC-seq and RNA-seq, we found that the Wnt signaling pathway, the mTOR signaling pathway, and other classical pathways regulate skeletal muscle development. In addition, some pathways that might regulate skeletal muscle growth, such as parathyroid hormone synthesis, secretion, and action, synthesis and degradation of ketone bodies, and the thyroid hormone signaling pathway, which were significantly enriched. After further study, we identified a number of candidate genes (ASNS, CARNS1, G0S2, PPP1R14C, and SH3BP5) that might be associated with muscle development. We also found that the differential regulation of chromatin openness at the level of some genes was contrary to the differential regulation at the level of transcription, suggesting that transcription factors and transcriptional repressors may be involved in the regulation of gene expression. Our study provided an in-depth understanding of the mechanism behind the differences in muscle fibers from two species of pig and provided an important foundation for further research on improving the quality of pork.

Introduction

Porcine skeletal muscle is one of the main sources of protein for the Chinese people, and skeletal muscle accounts for 40% of the body weight (1). The development and growth of skeletal muscle determine the quality and yield of meat (2). With the improvement of people's living standards, meat of increasingly high quality is pursued, and the meat from Duroc and other lean pigs can no longer meet the people's needs. Previous researchers have found that meat quality was related to muscle fiber characteristics and intramuscular fat content (3, 4). The character of muscle fiber is an important factor affecting meat quality (5, 6). One reason is that muscle fiber size affects the size of the muscle bundle, resulting in coarser muscle transverse sections (7, 8). Another reason is that the size of muscle fiber is related to tenderness; the smaller the muscle fiber, the lower the shear force, the better the tenderness and thus the better the meat quality (6). The shear force of pork is related to tenderness. Some researchers found that the meat quality of Lantang pig was better than Landrace pig. The shear force of Lantang pig was 61.23N, which significantly lower than that of Landrace pig (78.8N) (9). Luchuan pig is a typical fat breed, which is characterized by good meat quality, high reproductive performance, and a slow growth rate (10). Meanwhile, the diameter and area of the longissimus dorsi muscle fibers are significantly smaller than those of Duroc pigs. They are good models for studying the molecular mechanism underlying differences in meat quality. Therefore, studies on the differences between the longissimus dorsi muscle fibers of Luchuan and Duroc pigs can provide a direction for improving meat quality.

The muscle fiber diameter is related to muscle growth and development. It has been reported that the growth and development of skeletal muscle were regulated by interactions of multiple functional genes. The roles of transcription factors of the myogenic regulatory factor (MRF) family (11) and the myocyte enhancer factor-2 (MEF2) family (12) in the development of skeletal muscle have been intensively studied. Members of the MRF family, including myogenic differentiation antigen (MYOD), myogenic factor 5 (MYF5), myogenin (MYOG), and myogenic factor 6 (MYF6), which are key inducers of skeletal muscle development (13), are characterized by basic helix-loop-helix (bHLH) structures. The MEF2 family also has four members, namely MEF2a, MEF2b, MEF2c, and MEF2d. They all share an N-terminal MADS domain and a binding domain of MEF2, which can specifically bind A/T-enriched conserved sequences (14). There is a direct interaction between the MEF2 gene family and the MRF gene family (15). They can regulate each other's expression through the bHLH and MADS domains, and in most cases, the activation of genes related to muscle development is regulated synergistically. Many classical pathways such as the Notch signaling pathway (16), the Wnt signaling pathway (17), and the mTOR signaling pathway (18) are all key pathways for muscle proliferation and differentiation. In addition, some hormones such as parathyroid hormone (19) and thyroid hormone (20) also regulate skeletal muscle development. Recently, some researchers have studied the skeletal muscle development at different embryonic stages of the same pig breed from the perspective of chromatin openness and RNA sequencing (RNA-seq). However, the molecular mechanisms underlying differences in muscle growth between breeds remain unclear. Particularly, the molecular mechanisms underlying differences in meat quality have not been analyzed from the perspective of chromatin openness characteristics.

In this study, we performed an assay for transposase-accessible chromatin using sequencing (ATAC-seq) (21) for the detection of open chromatin regions of the longissimus dorsi muscle in two pig species in order to identify transcription factors that regulate muscle growth and development. By comparing ATAC-seq results with RNA-seq data, we further identified a number of potential core genes that regulated muscle proliferation and differentiation. In addition, some pathways that regulated this process were identified. These results will provide a useful resource for further exploration of the mechanisms underlying differences in meat quality.

Materials and Methods

Animals and Sample Collection

Three Luchuan and three Duroc boars were used in this study. Animals had free access to food and water and were kept under the same conditions. All pigs were put to death at the age of 180 days, with overnight fasting prior to death. Samples of the longissimus dorsi muscle were collected and frozen in liquid nitrogen for subsequent experiments. All animal experiments were approved by the Institutional Animal Care and Use Committee of Guangxi University (GXU-2015-003).

Histological Analysis

To investigate the size of muscle fibers, the muscle was immobilized in 4% paraformaldehyde solution for 24 h and completely embedded in paraffin. A solidified paraffin block was cut crosswise and stained with hematoxylin and eosin (H&E). Histological examination was conducted under a light microscope. Three tissue sections of each pig were included, and five visual fields were randomly selected for each section under 40× magnification. Images were collected using an upright microscope. At least 50 muscle fibers were randomly selected by Image-Pro Plus 6.0 Image analysis software, and the diameter and area of muscle fibers in the collected images were measured.

ATAC-seq

The original Illumina high-throughput sequencing image data files were converted into original sequences after base group recognition by bcl2fastq software (22). The results were stored in FASTQ file format. Trimmomatic (v0.36) software was used to filter the original sequences to obtain high quality clean reads, and Bowtie2 software was used to compare the sequencing sequences with the reference genome (Sus scrofa 11.1) (23). MACS2 software was used for peak detection to obtain an overview of the open chromatin regions of the whole genome of each sample (24). Different peaks were identified based on the following criteria: |log₂(fold change)| ≥ 1 and P < 0.00001. The genes corresponding to the different peaks of each sample were identified and annotated. Three biological replicates were used. The motif was examined using the Multiple EM for Motif Elicitation (MEME) suite and transcription factors were predicted using the motif database scanning algorithm TOMTOM (25).

RNA-seq

Transcriptome data were assembled using StringTie and compared with the reference genome (Sus scrofa 10.2) using CuffCompare (26). The calculated gene expression was directly used to compare gene expression between different samples. Difference analysis software DEGseq was used to analyze the differences between groups. Significantly differentially expressed genes (DEGs) were identified based on the following criteria: |log₂(fold change)| ≥ 1 and q < 0.001. Three biological replicates were used.

Integration Analysis of ATAC-seq and RNA-seq

The genes associated with the open chromatin regions in Duroc and Luchuan overlapped with the up and down regulated genes in the Duroc and Luchuan transcriptomes, respectively.

Gene Functional Annotation

Gene Ontology (GO) enrichment analysis of DEGs was performed using topGO software. KOBAS software was used for Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis among DEGs.

Analysis of Gene Expression

Three duplicate tissue samples were obtained from each pig. Total RNA was isolated using TRIzol reagent (OMEGA, USA; Genstar, Beijing, China) and then reverse transcribed using a first strand cDNA synthesis kit (Takara, Dalian, China) according to the manufacturer's instructions. The primers used in this study are shown in Table 1. Thermal cycling conditions were as follows: 95°C for 1 min, followed by 40 cycles of 10 s at 95°C, 34 s at 60°C, and 1 min at 60°C. Three biological replicates were used. Primers for selected genes were designed by primer 5.0. Cycle threshold values were collected and normalized to that of TBP.

TABLE 1

Table 1. qPCR primer sequence.

Statistics

All data are expressed as mean ± SEM. Statistical analysis was performed by the unpaired two-tailed Student's t-test. The probability level of P < 0.05 (^*) or P < 0.01 (^**) is used to indicate significance.

Results

Phenotype and ATAC-seq Quality Control of the Duroc and Luchuan Pig Muscle Tissues

It is well-known that the meat quality of Luchuan pigs is higher than that of Duroc pigs, partially because the diameter and area of muscle fiber of Luchuan pigs are smaller than those in Duroc pigs (Figure 1A).

FIGURE 1

Figure 1. Phenotype and ATAC-seq quality control. (A) Muscle tissue sections at 40× and 200× magnification. Quantification of the diameter and area of muscle fibers. The data are expressed as mean ± SD. **P < 0.01. (B) Distribution of insert size. (C) Enrichment of ATAC-seq signals around the TSS. The x-axis represents the normalized gene or peak length, and the y-axis represents the read enrichment. The larger the value, the higher the enrichment. TSS, transcription start site; TES, transcription end site. −2.0 represents 2 kb upstream of the TSS, and 2.0 represents 2 kb downstream of the TES. (D) The Pearson's correlation results shown by a heat map scatterplot. (E) Principal component analysis (PCA).

To unveil the mechanism behind the muscle fiber difference, ATAC-seq was performed to investigate the chromatin accessibility of the genome. A total of 767,777,278 raw reads were obtained. After filtration, 723,124,192 clean reads were obtained (Table 2). All libraries produced fragment lengths with expected distribution (Figure 1B). The highest peak on the left was a nucleosome-free fragment corresponding to the open chromatin region (the nucleosome cleaved fragment <147), and the highest peak on the right was the open chromatin fragment representing the fragment including two nucleosomes (the nucleosome cleaved fragment >147 and <147 × 2), which is consistent with the characteristics of open region of chromatin revealed by other studies. Most of the identified accessible areas were enriched in 2 kb of transcription start site (Figure 1C), indicating that open regions of chromatin were involved in transcriptional regulation. To further determine the quality of our ATAC-seq data, pairwise Pearson correlation between any pair of ATAC-seq samples was calculated according to the read signals from the combined ATAC-seq peaks of all samples, and the results showed that there were similarities between pigs of the same species and differences between the two breeds (Figure 1D). Principal component analysis (PCA) revealed the similarity between biological replicates and the difference between the two pig breeds (Figure 1E). All these results suggested that the quality of the sequencing data was very high.

TABLE 2

Table 2. ATAC-seq data statistics.

Chromatin Accessibility in the Muscle of Duroc and Luchuan Pigs

We identified 44,532 peaks specific to Duroc, 26,773 peaks specific to Luchuan, and 68,776 common peaks. The MA plot in Figure 2A shows the relationship between the open intensity and fold change of all peaks. As can be seen from the distribution map of all peaks on the chromosome, the signal intensity of Duroc peaks was higher than that of Luchuan peaks in total (Figure 2B). Genome-wide functional regions are divided into introns, intergenic regions, promoters, exons, 5'UTRs, and 3'UTRs to annotate all peaks (Figure 2B). We found that the peaks of the promoter regions in Luchuan accounted for 20.62% of the total area, while those of Duroc accounted for 16.90% (Figure 2C). As shown in the heat map, the peaks of Duroc and Luchuan were divided into four clusters. The clusters revealed differences in chromatin openness between the two pig breeds (Figure 2D). Compared with cluster1 and cluster2, the signal intensity was not only higher but also had the biggest difference. Compared with cluster3 and cluster4, although the signal intensity was different, the signal intensity was weaker than that of cluster1 and cluster2, which suggested that the key factors that might regulate muscle development were these strong signal peaks. There was a big difference between the genome-wide Duroc peak and Luchuan peak, but whether there was a chromosomal difference is not clear. Thus, the distribution of Duroc peaks and Luchuan peaks on chromosomes (the length of each different peak was standardized according to the length of each chromosome) was further studied. It was found that the peaks of Duroc were much more enriched on the X chromosome than those of Luchuan. It seems that downregulated chromatin accessibility regions most likely enriched in X chromosome. Then, the genes corresponding to downregulated chromatin accessibility regions on chromosome X were subjected to KEGG pathway enrichment analysis (Figure 2E). Our data revealed that regulation of the actin cytoskeleton, the Hippo signaling pathway, and the Wnt signaling pathway were enriched pathways involved in muscle development, which indicated that the X chromosome was an important chromosome that regulated muscle development, which is consistent with other studies (27, 28).

FIGURE 2

Figure 2. Analyses of the peaks. (A) MA plot. The x-axis represents signal strength in the open region of chromatin, and the y-axis represents fold change. (B) Chromosomal distribution of all peaks. (C) Distribution ratio of peaks in gene functional elements. (D) Top, enrichment of ATAC-seq signals around the 3 kb upper and lower reaches of the peak center. Bottom, peaks are divided into four categories according to the signal intensity. (E) Left, chromosome distribution of Duroc peaks and Luchuan peaks (P < 0.01). Right, KEGG pathway enrichment analysis of the genes corresponding to Duroc peaks on the X chromosome.

Corresponding Genes and Motif Analysis in Different Peaks

In order to identify the open sites related to phenotype, different peaks were filtered with P < 0.00001 |log₂(fold change)| ≥ 1. There were 3,412 different peaks in duroc pigs, 1,639 different peaks in Luchuan pigs (Figure 3A). By annotating each different peak, the different peaks of Duroc correspond to 1,376 genes, and the different peaks of Luchuan correspond to 628 genes, and a total of 1,905 genes were obtained. GO analysis and KEGG pathway enrichment analysis were conducted with these genes. Pathways with P ≤ 0.05 were considered significantly enriched. In the biological process category of GO enrichment, mainly processes related to muscle cell differentiation and muscle development were enriched, among which actin cytoskeleton organization was the most significant (Figure 3B). Moreover, 74 KEGG pathways were significantly enriched (P ≤ 0.05), many of which are involved in muscle development, such as the Wnt signaling pathway, the AMPK signaling pathway, and the Hippo signaling pathway (29–31). We showed 12 classical pathways associated with muscle development in Figure 3C. Duroc muscles are better developed than those of Luchuan. On this foundation, we analyzed the motifs on the different peaks of Duroc using the MEME suite and predicted the transcription factors by using the motif database scanning algorithm TOMTOM. Finally, 373 transcription factors (P < 0.05) were predicted. The transcription factors MYOD, MYOG, MYF5, MYF6, and members of the MEF2 family were key factors, which regulated the expression of many genes related to muscle development. Based on our data, they may bind to a total of 99 sites in different Duroc peaks (Figure 3D). Interestingly, the MEF2 family did not intersect with the other four transcription factors, suggesting it was regulated to muscle development in different ways. Among all the predicted transcription factors, 160 transcription factors shared common binding sites with these key factors. Sorted by the number of common binding sites, the top 10 transcription factors are shown in Figure 3E. Interestingly, nine of these 10 transcription factors have been reported to be associated or potentially associated with muscle development.

FIGURE 3

Figure 3. Different peaks corresponding to target gene enrichment analyses and transcription factor prediction. (A) Scatter plot of differential peaks (P < 0.00001, |log₂(fold change)| ≥ 1). (B) GO enrichment analysis of target genes corresponding to different peaks. (C) KEGG pathway enrichment analysis of target genes corresponding to different peaks. (D) Unified Venn diagrams of predicted binding sites of star factors (MRFs and MEF2s) involved in the regulation of muscle development according to motif alignment. (E) Predicted transcription factors. Left, the number of co-binding sites between the predicted transcription factors and the star factors. Right, the motif type of the transcription factors.

RNA-seq Data From Duroc and Luchuan Pig Muscle Tissues

High-throughput mRNA sequencing was performed to evaluate the transcriptome. PCA showed that there was a correlation within groups and a difference between groups of mRNA abundance among biological repeats (Supplementary Figure 1A). In total, 22,787 genes were identified in Duroc and 22,778 genes were identified in Luchuan by RNA-seq. The abundance of mRNA in two groups followed an approximately lognormal distribution (Supplementary Figures 1B,C). In total, 21,560 genes were shared by Duroc and Luchuan (Supplementary Figure 1D). The average abundance of mRNA in Duroc had a high correlation with that in Luchuan (R² = 0.899; Supplementary Figure 1E).

In order to identify the crucial functional genes related to phenotype, differentially expressed mRNAs (DEmRNAs) were filtered based on the criteria |log₂(fold change)| ≥ 1 and q < 0.001. Overall, 1,680 upregulated and 1,127 downregulated mRNAs were screened (Figure 4A). GO analysis and KEGG pathway enrichment analysis were used to analyze these up and downregulated mRNAs. All DEmRNAs were mapped to terms of the GO database, enabling annotation. These upregulated DEmRNAs were enriched in protein-related processes, etc. (Figure 4B). These downregulated DEmRNAs were enriched in ATP binding, etc. (Figure 4C). Some KEGG pathways related to protein metabolism and energy metabolism were significantly enriched. These upregulated DEmRNAs were enriched in cAMP signaling pathway, etc. (Figure 4D). These downregulated DEmRNAs were enriched in protein digestion and absorption, starch and sucrose metabolism, the mTOR signaling pathway, ABC transporters, and the TGF-beta signaling pathway (Figure 4E) (29–31). This suggested that DEmRNAs may be involved in maintaining intramuscular homeostasis and muscle development.

FIGURE 4

Figure 4. Analyses of RNA-seq. (A) Scatter plot of transcriptome data (q < 0.001, |log₂(fold change)| ≥ 1). (B) GO enrichment analysis of upregulated mRNAs. (C) GO enrichment analysis of downregulated mRNAs. (D) KEGG classification of upregulated mRNAs. (E) KEGG classification of downregulated mRNAs.

Integration of ATAC-seq Results With RNA-seq

To determine whether changes in open chromatin regions in ATAC-seq analysis were associated with changes in gene expression, we integrated ATAC-seq data with RNA-seq data. In total, 234 DEGs were identified (Figure 5A), including 121 upregulated genes and 113 downregulated genes.

FIGURE 5

Figure 5. Analysis results of integrated ATAC-seq and RNA-seq results. (A) Intersection of different mRNAs and different peaks corresponding to the target genes. (B) GO enrichment of upregulated DEGs. (C) GO enrichment of downregulated DEGs. (D) KEGG classification of upregulated DEGs. (E) KEGG classification of downregulated DEGs. (F) Peak diagram showing the relationship between genes in open regions of chromatin and gene expression. Among them, the PVALB gene of Duroc had a different peak and its transcriptional expression level was high; the PDK4 gene of Luchuan had a different peak and its transcriptional expression level was high; the CARNS1 gene of Duroc had different peaks but its transcriptional expression level was low. (G) A radar map of the top 20 candidate genes.

According to the GO annotation of these DEGs, it is interesting to note that the 121 upregulated genes were significantly enriched in pathways related to muscle homeostasis and muscle maintenance, including ion homeostasis and cellular homeostasis (Figure 5B). Meanwhile, the 113 downregulated genes were enriched in the Wnt signaling pathway, ATP binding, the Notch signaling pathway, and so on, which are related to muscle development regulation (Figure 5C). After annotating the DEGs to KEGG pathways, the data showed that some pathways related to energy metabolism and protein metabolism, such as the cAMP signaling pathway, galactose metabolism, starch and sucrose metabolism, alanine, aspartate, and glutamate metabolism, and phenylalanine, tyrosine, and tryptophan biosynthesis were significantly enriched (Figures 5D,E). Interestingly, parathyroid hormone synthesis, secretion, and action was the most significantly enriched pathway. In order to further determine the relationship between chromatin openness and gene expression and to further investigate how transcription factors regulate downstream genes, we analyzed the relationship between gene openness and gene expression. As shown in Table 3, 67 genes were upregulated and 17 genes were downregulated in the region with high chromatin openness, while 96 genes were downregulated and 54 genes were upregulated in the region with low chromatin openness. This may be due to the involvement of transcriptional factors or transcriptional repressors in the regulation of gene expression. Some examples are shown in Figure 5F. The chromatin accessibility and transcription level of PVALB in Duroc were higher than those in Luchuan. On the contrary, the chromatin accessibility and transcription level of PDK4 in Luchuan were higher than those in Duroc. Interestingly, the chromatin accessibility of CARNS1 in Duroc was higher than that in Luchuan, however, the transcription level was lower. According to the P-value, we focused on the top 20 DEGs (Figure 5G). Interestingly, nine of the 20 candidate genes have been reported to be associated with muscle development, and five have been reported to be potentially associated with muscle development.

TABLE 3

Table 3. Chromatin openness of differentially expressed genes.

To validate the accuracy of the RNA-seq data, eight genes randomly selected from the top 20 DEGs (PVALB, LGMN, PPP1R14C, EGF, ASNS, G0S2, CCNYL1, and PDK4) were analyzed by qRT-PCR. These total RNA samples were taken from the longissimus dorsi muscle of three Duroc and three Luchuan pigs. The gene expression patterns were similar in our RNA-seq and qRT-PCR analyses (Figure 6), which confirmed the accuracy of our RNA-seq data in both pig muscles.

FIGURE 6

Figure 6. Verification of RNA-seq data by qRT-PCR. The left y-axis represents the relative expression level as determined by qRT-PCR, and the right y-axis represents the FPKM value as determined by RNA-seq. All data represent the average value of three biological replicates. The error line represents three repeated standard errors, and all data are normalized. *P < 0.05, **P < 0.01.

Discussion

Pig skeletal muscle growth and development is one of the key factors that influence the quality of pork. Genome-wide chromatin accessibility regulate Luchuan pigs and Duroc pigs skeletal muscle growth and development by binding different TFs. The number of open sites and the signal intensity had great difference partly interpret the difference in muscle development between the two pig breeds (32). At the same time, we identified pathways and genes that may regulate the growth and development of skeletal muscle through ATAC-seq and RNA-seq analysis.

ATAC-seq identified three methylated open chromatin regions of H3K4, H3K36, and H3K79 (33). MRF and MEF2 transcription factors were involved in muscle proliferation and differentiation, and they all had a bHLH structure. In general, all bHLH transcription factors have a common motif or domain. Some transcription factors that share common binding sites with MRF and MEF2 families were screened out by the co-localization method. The more common binding sites there were, the more likely they were to have similar functions (34). We focused on the top 10 transcription factors with the most binding sites. Interestingly, these 10 transcription factors all had a bHLH structure, among which four transcription factors, including TCF3, ID4, ASCL2, and TCF4, have already been reported to regulate muscle development. The TCF3 protein forms a heterodimer with muscle regulators such as MYOD, which then binds to DNA and regulates the transcription of target genes related to muscle differentiation (35). ID4 inhibits DNA binding of E47 homotypes and E47/MYOD heterodimers, and is a key transcriptional regulator of muscle-related genes (36). The expression of ASCL2 is regulated by the Notch signal, and inhibits myogenesis by antagonizing the transcriptional activity and blocking the regeneration of injured muscles (37). TCF4 is strongly expressed in connective tissue fibroblasts, and low levels of TCF4 in myoblasts promote slow and rapid myogenesis, thereby promoting the overall maturation of muscle fiber types (38, 39). Besides, five other transcription factors, including TFAP4 (40), MAX (41, 42), NHLH1 (43), FRX5 (44, 45), and TGIF1 (46, 47), have potential roles in muscle development. Whether ATOH1 is involved in muscle development is still unclear.

The integration of ATAC-seq and RNA-seq results showed that 121 genes were upregulated and 113 genes were downregulated in Luchuan vs. Duroc. We focused on the top 20 genes with transcriptome P-values. Nine of these genes, including BCAP31 (48), CD59 (49), EGF (50), GLUL (51), GOT1 (52), MYBPC1 (53), PDK4 (54), PVALB (55), and SHISA2 (56), have been reported to be involved in muscle development. Five genes, namely ASNS (57), CARNS1 (58), G0S2 (59), PPP1R14C (60, 61), and SH3BP5 (62), have been reported to be potentially associated with muscle development. For six other genes, including ABCA5, CCNYL1, LGMN, ENSSSCG00000005481, RDH16, and TDRD7, so far there are no reports that showed a link with muscle development, thus their function in muscle development and function remains to be explored in the future.

Skeletal muscle development is regulated by various signaling pathways. Our KEGG pathway analysis results revealed some classical pathways that participated in the regulation of muscle development (Figures 5D,E). In addition, we found that some hormones, such as parathyroid hormone and thyroid hormone, may also regulate muscle development, because organs involved in calcium and phosphorus metabolism are related to parathyroid hormone, which directly responds to myocyte factors (19). Skeletal muscle is the primary target of thyroid hormone signaling and thyroid hormone signaling is involved in important biological functions, including energy consumption, heat production, development and growth (20). These findings provided new direction for further exploring meat quality.

Furthermore, ATAC-seq and RNA-seq integration analysis showed that not all genes with increased chromatin openness had increased gene expression. Some genes showed the opposite result, which may be caused by the binding of transcription repressors to the open chromatin area (63). Our data indicated that chromatin openness was not directly correlated with gene expression. The specific reasons remain to be further studied.

Conclusion

Overall, in this study we present a novel method for identifying regions of open chromatin and predicting transcription factors (TFAP4, MAX, NHLH1, FRX5, and TGIF1), involved in the regulation of muscle development in different breeds of pigs. With the combination of ATAC-Seq and RNA-Seq, we identified several candidate genes (ASNS, CARNS1, G0S2, PPP1R14C, and SH3BP5) that may regulate muscle development. In addition to this, we found that muscle development was maybe related to some of the hormonal signaling pathways (parathyroid hormone synthesis and action and thyroid hormone signaling pathway) in two breeds of pigs. This study can therefore be used as a reference for future research on the meat quality differences between Luchuan and Duroc pigs.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI GEO; GSE180840.

Ethics Statement

All animal experiments were approved by the Institutional Animal Care and Use Committee of Guangxi University (GXU-2015-003).

Author Contributions

WM, ZM, and LZ conceived the project, designed the protocol, and wrote the manuscript. ZT, LY, SL, and TH analyzed the data. PW, TW, ZS, HZ, and YL performed the experiments. All authors have read and approved the final manuscript.

Funding

This work was supported by grants from the National Key R&D Program of China (2018YFD0500402), Guangxi Science Foundation for Distinguished Young Scholars (2020GXNSFFA297008), Guangxi Science and Technology Base and Talents Project (AD18281085), Guangxi Natural Science Foundation (2019GXNSFDA245029), Guangxi Hundred-Talent Program, State Key Laboratory for Conservation and Utilization of Subtropical Agro-bioresources (SKLCUSA-a202006), and Training Project of High-level Professional and Technical Talents of Guangxi University.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnut.2021.742672/full#supplementary-material

Supplementary Figure 1. The mRNA data characteristics. (A) Principal component analysis (PCA). (B,C) Distribution of the mRNA abundance in Duroc and Luchuan. Genes were classified based on their rounded log₂(FPKM) abundance. Each bar represents the gene count as detected by RNA-seq in each category. (D) The intersection of genes expressed in Duroc and Luchuan. (E) Correlation between Duroc abundance and Luchuan abundance.

References

1. Zhang X, Sun W, He L, Wang L, Qiu K, Yin J. Global DNA methylation pattern involved in the modulation of differentiation potential of adipogenic and myogenic precursors in skeletal muscle of pigs. Stem Cell Res Ther. (2020) 11:536. doi: 10.1186/s13287-020-02053-3

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Fu Y, Shang P, Zhang B, Tian X, Nie R, Zhang R, et al. Function of the porcine TRPC1 gene in myogenesis and muscle growth. Cells Basel. (2021) 10:147. doi: 10.3390/cells10010147

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Zhang Z, Liao Q, Sun Y, Pan T, Liu S, Miao W, et al. Lipidomic and transcriptomic analysis of the longissimus muscle of Luchuan and Duroc pigs. Front Nutr. (2021) 8:667622. doi: 10.3389/fnut.2021.667622

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Joo ST, Kim GD, Hwang YH, Ryu YC. Control of fresh meat quality through manipulation of muscle fiber characteristics. Meat Sci. (2013) 95:828–36. doi: 10.1016/j.meatsci.2013.04.044

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Kim GD, Jeong JY, Jung EY, Yang HS, Lim HT, Joo ST. The influence of fiber size distribution of type IIB on carcass traits and meat quality in pigs. Meat Sci. (2013) 94:267–73. doi: 10.1016/j.meatsci.2013.02.001

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Huo W, Weng K, Gu T, Zhang Y, Zhang Y, Chen G, et al. Effect of muscle fiber characteristics on meat quality in fast- and slow-growing ducks. Poult Sci. (2021) 100:101264. doi: 10.1016/j.psj.2021.101264

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Chen XD, Ma QG, Tang MY, Ji C. Development of breast muscle and meat quality in Arbor Acres broilers, Jingxing 100 crossbred chickens and Beijing fatty chickens. Meat Sci. (2007) 77:220–7. doi: 10.1016/j.meatsci.2007.03.008

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Kokoszynski D, Piwczynski D, Arpasova H, Hrncar C, Saleh M, Wasilewski R. A comparative study of carcass characteristics and meat quality in genetic resources Pekin ducks and commercial crossbreds. Asian Australas J Anim Sci. (2019) 1753–62. doi: 10.5713/ajas.18.0790

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Yu K, Shu G, Yuan F, Zhu X, Gao P, Wang S, et al. Fatty acid and transcriptome profiling of longissimus dorsi muscles between pig breeds differing in meat quality. Int J Biol Sci. (2013) 9:108–18. doi: 10.7150/ijbs.5306

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Ran ML, He J, Tan JY, Yang AQ, Li Z, Chen B. The complete sequence of the mitochondrial genome of Luchuan pig (Sus scrofa). Mitochondrial DNA A DNA Mapp Seq Anal. (2016) 27:1880–1. doi: 10.3109/19401736.2014.947588

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Aase-Remedios ME, Coll-Llado C, Ferrier D. More than one-to-four via 2R: evidence of an independent amphioxus expansion and two-gene ancestral vertebrate state for MyoD-related myogenic regulatory factors (MRFs). Mol Biol Evol. (2020) 37:2966–82. doi: 10.1093/molbev/msaa147

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Taylor MV, Hughes SM. Mef2 and the skeletal muscle differentiation program. Semin Cell Dev Biol. (2017) 72:33–44. doi: 10.1016/j.semcdb.2017.11.020

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Kassar-Duchossoy L, Gayraud-Morel B, Gomes D, Rocancourt D, Buckingham M, Shinin V, et al. Mrf4 determines skeletal muscle identity in Myf5:Myod double-mutant mice. Nature. (2004) 431:466–71. doi: 10.1038/nature02876

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Zia A, Imran M, Rashid S. In silico exploration of conformational dynamics and novel inhibitors for targeting MEF2-associated transcriptional activity. J Chem Inf Model. (2020) 60:1892–909. doi: 10.1021/acs.jcim.0c00008

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Brand-Saberi B. Genetic and epigenetic control of skeletal muscle development. Ann Anat. (2005) 187:199–207. doi: 10.1016/j.aanat.2004.12.018

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Buas MF, Kadesch T. Regulation of skeletal myogenesis by Notch. Exp Cell Res. (2010) 316:3028–33. doi: 10.1016/j.yexcr.2010.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Girardi F, Le Grand F. Wnt signaling in skeletal muscle development and regeneration. Prog Mol Biol Transl Sci. (2018) 153:157–79. doi: 10.1016/bs.pmbts.2017.11.026

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Sakai H, Murakami C, Usuki T, Lu Q, Matsumoto KI, Urano T, et al. Diacylglycerol kinase eta regulates C2C12 myoblast proliferation through the mTOR signaling pathway. Biochimie. (2020) 177:13–24. doi: 10.1016/j.biochi.2020.07.018

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Lombardi G, Ziemann E, Banfi G, Corbetta S. Physical activity-dependent regulation of parathyroid hormone and calcium-phosphorous metabolism. Int J Mol Sci. (2020) 21:5388. doi: 10.3390/ijms21155388

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Salvatore D, Simonides WS, Dentice M, Zavacki AM, Larsen PR. Thyroid hormones and skeletal muscle–new insights and potential implications. Nat Rev Endocrinol. (2014) 10:206–14. doi: 10.1038/nrendo.2013.238

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Shashikant T, Ettensohn CA. Genome-wide analysis of chromatin accessibility using ATAC-seq. Methods Cell Biol. (2019) 151:219–35. doi: 10.1016/bs.mcb.2018.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Buenrostro JD, Wu B, Chang HY, Greenleaf WJ. ATAC-seq: a method for assaying chromatin accessibility genome-wide. Curr Protoc Mol Biol. (2015) 109:21.29.1–9. doi: 10.1002/0471142727.mb2129s109

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. (2012) 9:357–9. doi: 10.1038/nmeth.1923

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. (2008) 9:R137. doi: 10.1186/gb-2008-9-9-r137

PubMed Abstract | CrossRef Full Text

25. Bailey TL, Williams N, Misleh C, Li WW. MEME: discovering and analyzing DNA and protein sequence motifs. Nucleic Acids Res. (2006) 34:W369–73. doi: 10.1093/nar/gkl198

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. (2015) 33:290–5. doi: 10.1038/nbt.3122

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Viggiano E, Ergoli M, Picillo E, Politano L. Determining the role of skewed X-chromosome inactivation in developing muscle symptoms in carriers of Duchenne muscular dystrophy. Hum Genet. (2016) 135:685–98. doi: 10.1007/s00439-016-1666-6

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Moreira-Marconi E, Sa-Caputo DC, Dionello CF, Guedes-Aguiar EO, Sousa-Goncalves CR, Morel DS, et al. Whole-body vibration exercise is well tolerated in patients with Duchenne muscular dystrophy: a systematic review. Afr J Tradit Complement Altern Med. (2017) 14(4 Suppl.):2–10. doi: 10.21010/ajtcam.v14i4S.1

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Wang J, Ren Q, Hua L, Chen J, Zhang J, Bai H, et al. Comprehensive analysis of differentially expressed mRNA, lncRNA and circRNA and their ceRNA Networks in the Longissimus Dorsi muscle of two different pig breeds. Int J Mol Sci. (2019) 20:1107. doi: 10.3390/ijms20051107

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Liu Y, Yang X, Jing X, He X, Wang L, Liu Y, et al. Transcriptomics analysis on excellent meat quality traits of skeletal muscles of the chinese indigenous min pig compared with the large white breed. Int J Mol Sci. (2017) 19:21. doi: 10.3390/ijms19010021

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Cai C, Li M, Zhang Y, Meng S, Yang Y, Gao P, Guo X, et al. Comparative transcriptome analyses of Longissimus thoracis between pig breeds differing in muscle characteristics. Front Genet. (2020) 11:526309. doi: 10.3389/fgene.2020.526309

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Yue J, Hou X, Liu X, Wang L, Gao H, Zhao F, et al. The landscape of chromatin accessibility in skeletal muscle during embryonic development in pigs. J Anim Sci Biotechnol. (2021) 12:56. doi: 10.1186/s40104-021-00577-z

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Lowe EK, Cuomo C, Voronov D, Arnone MI. Using ATAC-seq and RNA-seq to increase resolution in GRN connectivity. Methods Cell Biol. (2019) 151:115–26. doi: 10.1016/bs.mcb.2018.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Jolma A, Yin Y, Nitta KR, Dave K, Popov A, Taipale M, et al. DNA-dependent formation of transcription factor pairs alters their binding specificity. Nature. (2015) 527:384–8. doi: 10.1038/nature15518

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Sun L, Trausch-Azar JS, Ciechanover A, Schwartz AL. E2A protein degradation by the ubiquitin-proteasome system is stage-dependent during muscle differentiation. Oncogene. (2007) 26:441–8. doi: 10.1038/sj.onc.1209793

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Miyata K, Miyata T, Nakabayashi K, Okamura K, Naito M, Kawai T, et al. DNA methylation analysis of human myoblasts during in vitro myogenic differentiation: de novo methylation of promoters of muscle-related genes and its involvement in transcriptional down-regulation. Hum Mol Genet. (2015) 24:410–23. doi: 10.1093/hmg/ddu457

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Wang C, Wang M, Arrington J, Shan T, Yue F, Nie Y, et al. Ascl2 inhibits myogenesis by antagonizing the transcriptional activity of myogenic regulatory factors. Development. (2017) 144:235–47. doi: 10.1242/dev.138099

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Mathew SJ, Hansen JM, Merrell AJ, Murphy MM, Lawson JA, Hutcheson DA, et al. Connective tissue fibroblasts and Tcf4 regulate myogenesis. Development. (2011) 138:371–84. doi: 10.1242/dev.057463

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Murphy MM, Lawson JA, Mathew SJ, Hutcheson DA, Kardon G. Satellite cells, connective tissue fibroblasts and their interactions are crucial for muscle regeneration. Development. (2011) 138:3625–37. doi: 10.1242/dev.064162

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Cramer ML, Xu R, Martin PT. Soluble heparin binding epidermal growth factor-like growth factor is a regulator of GALGT2 expression and GALGT2-dependent muscle and neuromuscular phenotypes. Mol Cell Biol. (2019) 39:e00140-19. doi: 10.1128/MCB.00140-19

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Calura E, Cagnin S, Raffaello A, Laveder P, Lanfranchi G, Romualdi C. Meta-analysis of expression signatures of muscle atrophy: gene interaction networks in early and late stages. BMC Genomics. (2008) 9:630. doi: 10.1186/1471-2164-9-630

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Grandori C, Cowley SM, James LP, Eisenman RN. The Myc/Max/Mad network and the transcriptional control of cell behavior. Annu Rev Cell Dev Biol. (2000) 16:653–99. doi: 10.1146/annurev.cellbio.16.1.653

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Manetopoulos C, Hansson A, Karlsson J, Jonsson JI, Axelson H. The LIM-only protein LMO4 modulates the transcriptional activity of HEN1. Biochem Biophys Res Commun. (2003) 307:891–9. doi: 10.1016/S0006-291X(03)01298-1

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Buttice G, Miller J, Wang L, Smith BD. Interferon-gamma induces major histocompatibility class II transactivator (CIITA), which mediates collagen repression and major histocompatibility class II activation by human aortic smooth muscle cells. Circ Res. (2006) 98:472–9. doi: 10.1161/01.RES.0000204725.46332.97

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Londhe P, Davie JK. Gamma interferon modulates myogenesis through the major histocompatibility complex class II transactivator, CIITA. Mol Cell Biol. (2011) 31:2854–66. doi: 10.1128/MCB.05397-11

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Guca E, Sunol D, Ruiz L, Konkol A, Cordero J, Torner C, et al. TGIF1 homeodomain interacts with Smad MH1 domain and represses TGF-beta signaling. Nucleic Acids Res. (2018) 46:9220–35. doi: 10.1093/nar/gky680

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Abrigo J, Campos F, Simon F, Riedel C, Cabrera D, Vilos C, et al. TGF-beta requires the activation of canonical and non-canonical signalling pathways to induce skeletal muscle atrophy. Biol Chem. (2018) 399:253–64. doi: 10.1515/hsz-2017-0217

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Cacciagli P, Sutera-Sardo J, Borges-Correia A, Roux JC, Dorboz I, Desvignes JP, et al. Mutations in BCAP31 cause a severe X-linked phenotype with deafness, dystonia, and central hypomyelination and disorganize the Golgi apparatus. Am J Hum Genet. (2013) 93:579–86. doi: 10.1016/j.ajhg.2013.07.023

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Kaminski HJ, Kusner LL, Richmonds C, Medof ME, Lin F. Deficiency of decay accelerating factor and CD59 leads to crisis in experimental myasthenia. Exp Neurol. (2006) 202:287–93. doi: 10.1016/j.expneurol.2006.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Nagata Y, Ohashi K, Wada E, Yuasa Y, Shiozuka M, Nonomura Y, et al. Sphingosine-1-phosphate mediates epidermal growth factor-induced muscle satellite cell activation. Exp Cell Res. (2014) 326:112–24. doi: 10.1016/j.yexcr.2014.06.009

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Unal O, Ceylaner S, Akin R. A very rare etiology of hypotonia and seizures: congenital glutamine synthetase deficiency. Neuropediatrics. (2019) 50:51–3. doi: 10.1055/s-0038-1675637

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Wang Y, Zhao ZJ, Kang XR, Bian T, Shen ZM, Jiang Y, et al. lncRNA DLEU2 acts as a miR-181a sponge to regulate SEPP1 and inhibit skeletal muscle differentiation and regeneration. Aging. (2020) 12:24033–56. doi: 10.18632/aging.104095

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Stavusis J, Lace B, Schafer J, Geist J, Inashkina I, Kidere D, et al. Novel mutations in MYBPC1 are associated with myogenic tremor and mild myopathy. Ann Neurol. (2019) 86:129–42. doi: 10.1002/ana.25494

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Pin F, Novinger LJ, Huot JR, Harris RA, Couch ME, O'Connell TM, et al. PDK4 drives metabolic alterations and muscle atrophy in cancer cachexia. FASEB J. (2019) 33:7778–90. doi: 10.1096/fj.201802799R

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Ghahramani SM, Trollet C, Athanasopoulos T, Graham IR, Hu P, Dickson G. Transcriptomic analysis of dystrophin RNAi knockdown reveals a central role for dystrophin in muscle differentiation and contractile apparatus organization. BMC Genomics. (2010) 11:345. doi: 10.1186/1471-2164-11-345

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Liu Z, Wang C, Liu X, Kuang S. Shisa2 regulates the fusion of muscle progenitors. Stem Cell Res. (2018) 31:31–41. doi: 10.1016/j.scr.2018.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Brearley MC, Li C, Daniel Z, Loughna PT, Parr T, Brameld JM. Changes in expression of serine biosynthesis and integrated stress response genes during myogenic differentiation of C2C12 cells. Biochem Biophys Rep. (2019) 20:100694. doi: 10.1016/j.bbrep.2019.100694

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Wang-Eckhardt L, Bastian A, Bruegmann T, Sasse P, Eckhardt M. Carnosine synthase deficiency is compatible with normal skeletal muscle and olfactory function but causes reduced olfactory sensitivity in aging mice. J Biol Chem. (2020) 295:17100–13. doi: 10.1074/jbc.RA120.014188

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Ma T, Lopez-Aguiar AG, Li A, Lu Y, Sekula D, Nattie EE, et al. Mice lacking G0S2 are lean and cold-tolerant. Cancer Biol Ther. (2014) 15:643–50. doi: 10.4161/cbt.28251

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Lang I, Virk G, Zheng DC, Young J, Nguyen MJ, Amiri R, et al. The evolution of duplicated genes of the Cpi-17/Phi-1 (ppp1r14) family of protein phosphatase 1 inhibitors in teleosts. Int J Mol Sci. (2020) 21:5709. doi: 10.3390/ijms21165709

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Wu Y, Erdodi F, Muranyi A, Nullmeyer KD, Lynch RM, Hartshorne DJ. Myosin phosphatase and myosin phosphorylation in differentiating C2C12 cells. J Muscle Res Cell Motil. (2003) 24:499–511. doi: 10.1023/B:JURE.0000009810.36038.53

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Lamers I, Reijnders M, Venselaar H, Kraus A, Jansen S, de Vries B, et al. Recurrent de novo mutations disturbing the GTP/GDP binding pocket of RAB11B cause intellectual disability and a distinctive brain phenotype. Am J Hum Genet. (2017) 101:824–32. doi: 10.1016/j.ajhg.2017.09.015

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Truskina J, Han J, Chrysanthou E, Galvan-Ampudia CS, Laine S, Brunoud G, et al. A network of transcriptional repressors modulates auxin responses. Nature. (2021) 589:116–9. doi: 10.1038/s41586-020-2940-2

PubMed Abstract | CrossRef Full Text | Google Scholar