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

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.

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

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.

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

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.

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.

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

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.

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

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.

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.

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.

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.

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.