Analyzing super-enhancer temporal dynamics reveals potential critical enhancers and their gene regulatory networks underlying skeletal muscle development

  1. Yuwen Liu1,2,6
  1. 1Shenzhen Branch, Guangdong Laboratory for Lingnan Modern Agriculture, Key Laboratory of Livestock and Poultry Multi-Omics of MARA, Agricultural Genomics Institute at Shenzhen, Chinese Academy of Agricultural Sciences, Shenzhen 518124, China;
  2. 2Innovation Group of Pig Genome Design and Breeding, Research Centre for Animal Genome, Agricultural Genomics Institute at Shenzhen, Chinese Academy of Agricultural Sciences, Shenzhen 518124, China;
  3. 3Key Laboratory of Agricultural Animal Genetics, Breeding and Reproduction of Ministry of Education and Key Laboratory of Swine Genetics and Breeding of Ministry of Agriculture, College of Animal Science and Technology, Huazhong Agricultural University, Wuhan 430070, China;
  4. 4Guangdong Provincial Key Laboratory of Animal Molecular Design and Precise Breeding, School of Life Science and Engineering, Foshan University, Foshan 528225, China;
  5. 5Guangxi Key Laboratory of Animal Breeding, Disease Control and Prevention, College of Animal Science and Technology, Guangxi University, Nanning 530004, China;
  6. 6Kunpeng Institute of Modern Agriculture at Foshan, Agricultural Genomics Institute at Shenzhen, Chinese Academy of Agricultural Sciences, Foshan 528226, China

  1. ↵7 These authors contributed equally to this work.

  • Corresponding authors: tangzhonglin{at}caas.cn, liuyuwen{at}caas.cn
    •  
    Next Section

    Abstract

    Super-enhancers (SEs) govern the expression of genes defining cell identity. However, the dynamic landscape of SEs and their critical constituent enhancers involved in skeletal muscle development remains unclear. In this study, using pig as a model, we employed cleavage under targets and tagmentation (CUT&Tag) to profile the enhancer-associated histone modification marker H3K27ac in skeletal muscle across two prenatal and three postnatal stages, and investigated how SEs influence skeletal muscle development. We identify three SE families with distinct temporal dynamics: continuous (Con, 397), transient (TS, 434), and de novo (DN, 756). These SE families are associated with different temporal gene expression trajectories, biological functions, and DNA methylation levels. Notably, several lines of evidence suggest a potential prominent role of Con SEs in regulating porcine muscle development and meat traits. To pinpoint key cis-regulatory units in Con SEs, we developed an integrative approach that leverages information from eRNA annotation, genome-wide association study (GWAS) signals, and high-throughput capture self-transcribing active regulatory region sequencing (STARR-seq) experiments. Within Con SEs, we identify 20 candidate critical enhancers with meat and carcass-associated DNA variations that affect enhancer activity, and infer their upstream transcription factors and downstream target genes. As a proof of concept, we experimentally validate the role of one such enhancer and its potential target gene during myogenesis. Our findings reveal the dynamic regulatory features of SEs in skeletal muscle development and provide a general integrative framework for identifying critical enhancers underlying the formation of complex traits.

    Skeletal muscle, which constitutes nearly half of the body weight, is the most abundant tissue and plays essential roles in contraction, movement, energy storage, energy delivery, and homeostasis (Mukund and Subramaniam 2020). For the study of skeletal muscle development, the pig (Sus scrofa) offers significant advantages, due to its anatomical, structural, immunological, genomic, and physiological similarities to humans (Lunney et al. 2021). Additionally, muscle samples from pigs can be easily obtained at various prenatal and postnatal stages from genetically similar individuals.

    Porcine skeletal muscle (PSM) development is divided into prenatal myofiber formation and postnatal myofiber hypertrophy. Our previous study identified four important developmental phases: embryonic (33–65 days after conception, abbreviated as E33–E65), fetal (E70–E105), neonatal (0–60 days after birth, abbreviated as D0–D60), and adult (D80–D180) (Yang et al. 2021a). By E60 of the embryonic phase and E90 of the fetal phase, primary and secondary myofibers are fully formed, respectively (Beermann and Cassens 1977; Muráni et al. 2007). During the neonatal period, the transition between slow/fast myofiber types significantly affects pork quality (Zhang et al. 2022). The adult period, or growing–fattening stage (Van Erp-Van der Kooij et al. 2003; Scheffler et al. 2014), involves rapid myofiber hypertrophy with muscle satellite cells (MuSCs) proliferating and differentiating in response to high-nutrient diets and growth hormones (Brameld 1997; Bi and Kuang 2012). Understanding prenatal and postnatal PSM development is essential for elucidating the transcriptional regulatory networks that govern skeletal muscle growth and development.

    During skeletal muscle development, myogenic differentiation is driven by transcription factors (TFs) from the myogenic regulatory factor (MRF) family (MYOD1, MYOG, MYF5, and MYF6) and myocyte enhancer factor 2 (MEF2) family (MEF2A, MEF2B, MEF2C, and MEF2D) (Hernández-Hernández et al. 2017). These TFs bind cis-regulatory elements (CREs), such as enhancers, to spatiotemporally regulate gene expression (Yan et al. 2019; Zhao et al. 2021). One approach to identifying critical enhancers underlying PMS development is by locating super-enhancers (SEs), which are characterized by dense clusters of H3K27ac modifications and play pivotal roles in regulating genes that define cell identity (Hnisz et al. 2013; Whyte et al. 2013). Another approach involves integrating population genetics data with functional genomics data. Populations with enhancer DNA variations act as natural genetic screens, shedding light on associations between enhancers and their significant biological functions. Along this line of thinking, several studies have successfully identified critical enhancers that contribute to the genetic basis of complex traits (He et al. 2015; Aneas et al. 2021; Liu et al. 2022). However, the search for these enhancers remains in its early stages compared to the identification of critical genes.

    Fully exploiting enhancers to elucidate molecular mechanisms of complex traits relies on accurately and comprehensively annotating DNA regions with enhancer activity. Pioneers in this field include the Encyclopedia of DNA Elements (ENCODE), the modENCODE, and the Roadmap Epigenome Project, which characterize the enhancer landscape in humans and model organisms (Shen et al. 2012; The ENCODE Project Consortium 2012; Yue et al. 2014; Roadmap Epigenomics Consortium et al. 2015). The Functional Annotation of Animal Genomes (FAANG) consortium initiated similar enhancer profiling in livestock in 2019 (Giuffra et al. 2019), followed by other studies that collectively built annotations of potential enhancers in various pig tissues (Kern et al. 2021; Pan et al. 2021; Yuan et al. 2021; Yue et al. 2021; Zhao et al. 2021; Li et al. 2022; Zhu et al. 2022). However, in skeletal muscle, prior studies lack well-characterized temporal enhancer dynamics from prenatal to postnatal stages, and pinpointing critical trait-affecting enhancers remains largely unaddressed. In the present study, we aim to fill this gap by identifying potential critical enhancers and their gene regulatory networks during PSM development.

    Our approach involves analyzing developmental trajectories of PSM SEs and incorporating SE profiling with multiomics data, and high-throughput screening for DNA variants affecting enhancer activity (Fig. 1). We first performed cleavage under targets and tagmentation (CUT&Tag) (Kaya-Okur et al. 2019) to profile H3K27ac and identify SEs in PSM samples across myofiber formation and hypertrophy stages. By developing an innovative approach integrating SE temporal patterns, eRNA annotation, genome-wide association study (GWAS) signals, and self-transcribing active regulatory region sequencing (STARR-seq) experiments, Our study not only reveals the critical enhancer-mediated mechanisms underlying PSM development but also provides a general framework for investigating how enhancers contribute to the genetic basis of complex traits.

    Figure 1.
    View larger version:
    Figure 1.

    Workflow for revealing potential critical enhancers and trait-affecting gene regulatory networks during PSM development. (1) CUT&Tag libraries from frozen PSM samples across five developmental stages are produced and sequenced. (PSM) porcine skeletal muscle. (2) H3K27ac peaks are identified at each developmental stage and filtered for intragroup consistency. (IDR) irreproducibility discovery rate. (3) Identify SEs and typical enhancers (TEs) based on distal H3K27ac peaks. (ROSE) ranked order of super-enhancers, (TSSs) transcription start sites. (4) Analysis of enriched TF motifs on SEs. (5) Using hierarchical clustering to reveal temporal patterns of SEs. (6) Construct gene regulatory networks mediated by potential critical enhancers based on TF motif scanning and coexpression relationships among TFs, eRNAs, and target genes. (7) Refinement of trait-affecting regulatory network by identifying candidate critical enhancers within SEs, characterized by eRNA transcription and GWAS variants modulating enhancer activity. (8) Validate the role of candidate critical enhancers in myogenic differentiation using CRISPR–Cas9-mediated enhancer knockout in C2C12 cells.

    Previous SectionNext Section

    Results

    Dynamic H3K27ac landscape across prenatal and postnatal stages of PSM

    To investigate the enhancer landscape during myofiber formation and hypertrophy phases in skeletal muscle, we generated high-quality CUT&Tag H3K27ac (Supplemental Table S1) profiling in domestic pigs (Landrace) across five well-characterized developmental stages (E60, E90, D0, D60, and D180) (Fig. 2A). By combining two biological replicates via the irreproducible discovery rate (IDR) method (Li et al. 2011), we identified 15,390, 31,138, 30,817, 29,991, and 42,981 H3K27ac peaks at E60, E90, D0, D60, and D180, respectively (Supplemental Fig. S1A; Supplemental Table S2). These peaks overlapped with the majority of skeletal muscle-specific enhancers identified by Zhao et al. (2021), while exhibiting low overlap with their enhancers specific to other tissues (Supplemental Fig. S1B), reinforcing the robustness of our findings. Pearson's correlation analysis (Supplemental Fig. S1C) and principal component analysis (Supplemental Fig. S1D) for H3K27ac signals revealed distinct stage-dependent segregation, with differentiation of H3K27ac signals before and after birth (Supplemental Fig. S1E). This indicates significant variation in H3K27ac modification during PSM development. Further analysis of H3K27ac dynamics, conducted using Short Time-series Expression Miner (STEM) software (Ernst and Bar-Joseph 2006), revealed 14 distinct H3K27ac signal profiles, suggesting a complex functional landscape. Increasing and decreasing H3K27ac signals (Supplemental Fig. S1F) are associated with genes exhibiting corresponding changes in expression levels, respectively (Supplemental Fig. S1G,H; Supplemental Methods). Notably, elevated H3K27ac signals were associated with genes involved in muscle development, whereas reduced signals were associated with genes related to basal metabolism (Supplemental Fig. S1I; Supplemental Methods). Together, these findings indicate that H3K27ac signals and associated gene expression undergo significant changes throughout PSM development, potentially serving as key drivers of this biological process.

    Figure 2.
    View larger version:
    Figure 2.

    Identification and characterization of SEs. (A) Illustration of developmental PSM samples and three types of high-throughput sequencing data. CUT&Tag data were generated in this study, while whole genome bisulfite sequencing (WGBS) and strand-specific RNA-seq data were previously published by our group. (B) Classification of stitched enhancers as TEs and SEs based on H3K27ac signal (taking E90 stage as an example). Several muscle-specific genes proximal to SEs are highlighted. (C) Number of SEs identified at each developmental stage. (D) Comparison of gene expression associated with TEs and SEs at each stage. Statistical significance between groups was assessed using a two-sided Student's t-test. (*) P < 0.05, (**) P < 0.01. Unless explicitly stated, a two-sided Student's t-test is the default method for differential comparisons. (E) Top 30 TF motifs with the highest standard deviation of enrichment significance P in SEs across five developmental stages. Only TFs with significant motif enrichment (P < 1 × 10−10 at least at one stage) and sufficient mRNA abundance (transcripts per million [TPM] ≥ 3 at least at one stage) were retained. (F) Gene Ontology (GO) enrichment analysis of SEs at each stage.

    SEs exhibit temporal regulatory features associated with PSM development

    H3K27ac signals are canonical markers for active enhancers and promoters (Creyghton et al. 2010). We defined H3K27ac peaks within 1000 bp upstream and 500 bp downstream from gene transcription start sites (TSSs) as potential active promoters, while the remaining peaks were classified as potential active enhancer. Across the five developmental stages, enhancers exhibited a broader quantity spectrum and greater activity dynamics than promoters (Supplemental Fig. S2A,B). Enhancers were enriched with muscle-specific TF motifs and displayed developmental stage specificity (Supplemental Fig. S2C), whereas promoters mainly contained general TF motifs and lacked dynamic enrichment (Supplemental Fig. S2D). Additionally, the average number of enhancers per promoter increased from 1.35 to 3.60 during PSM development (Supplemental Fig. S2E), indicating growing redundancy and regulatory complexity among enhancers.

    To investigate critical enhancers, we used ROSE software for SE discovery (Whyte et al. 2013), analyzing IDR-consistent H3K27ac peaks. We observed several key MRF genes (MYOD1, MYOG, and MYF5) with adjacent SEs (Fig. 2B). The number of SEs increased gradually from 321 at E60 to 929 at D180 (Fig. 2C). SE-adjacent genes exhibited significantly higher expression levels than TE-adjacent genes, aligning with previous studies that SEs have a stronger cis transcription-stimulating effect than TEs (Fig. 2D). Although both SEs (Fig. 2E) and TEs (Supplemental Fig. S3A) displayed similar stage-specific enrichment patterns for muscle-specific TF motifs, they differed in GO term enrichment. SEs were increasingly enriched in categories such as muscle development, muscle differentiation, myofibril assembly, and muscle contraction (Fig. 2F), while TEs maintained a stable pattern (Supplemental Fig. S3B). This suggests that SEs play a more significant role in dynamically regulating critical genes during PSM development, a phase characterized by substantial shifts in transcriptional programs and cell states.

    To further investigate the functional significance of SEs compared to TEs, we analyzed their enrichment for various genetic features associated with regulatory activities (see Supplemental Methods). We curated GWAS hits associated with pig meat and carcass traits from the Animal Quantitative Trait Loci (QTL) Database (Animal QTLdb) (Hu et al. 2019) and cis-expression quantitative trait loci (cis-eQTLs) identified in skeletal muscle from the PigGTEx portal (Teng et al. 2024). Both TEs and SEs showed significant enrichment (Permutation test, P < 0.05) for these GWAS hits and cis-eQTLs, but SEs consistently displayed higher enrichment levels (Supplemental Fig. S3C,D), emphasizing their central role in gene regulation and trait formation within PSM. Furthermore, SE sequences showed higher enrichment for conserved DNA elements across 103 mammalian species from the Ensembl database compared to TEs (Supplemental Fig. S3E). Previous studies have highlighted the importance of conserved noncoding DNA elements in vertebrate development (Woolfe et al. 2005). These findings align with the established role of SEs in governing master regulators that define cell identity, prompting us to focus our analysis on critical enhancers within SEs.

    Continuous SEs may play a prominent role during PSM development

    A recent study has demonstrated that SEs exhibit differential temporal patterns associated with distinct biological functions in several cell differentiation models (Kai et al. 2021). Applying their method to standardize SE components based on tertiles of constituent units, we performed hierarchical clustering on 1587 merged SEs, revealing three distinct temporal patterns: 397 continuous (Con) SEs, 756 de novo (DN) SEs, and 434 transient (TS) SEs (Fig. 3A; Supplemental Table S3). The H3K27ac intensity of the “mother” unit, the constituent enhancer unit with the strongest H3K27ac signal in each SE (Bojcsuk et al. 2017), displayed distinct patterns: Con SEs maintained a consistently stable signal, DN SEs gradually increased, and TS SEs progressively declined (Fig. 3A,B). These SE categories were associated with different biological functions (Fig. 3C) and gene expression trajectories (Fig. 3D). Specifically, Con SEs were enriched in processes related to muscle cell differentiation and development, neighboring key MRF genes such as MYOD1, MYOG (Fig. 3B), and MYF5. DN SEs were associated with skeletal muscle contraction and metabolism, while TS SEs mainly contributed to embryonic morphogenesis and vascular development. Additionally, enhancer units in Con SEs exhibited greater fluctuations in H3K27ac activity than those in DN and TS SEs (Supplemental Fig. S4A,B), aligning with the association of dynamic H3K27ac peaks with muscle development and muscle cell differentiation (Supplemental Fig. S4C). Collectively, our findings suggest a crucial type of SEs characterized by stable mother units and highly dynamic constituent units.

    Figure 3.
    View larger version:
    Figure 3.

    Classification and characteristics of temporal SEs. (A) Hierarchical clustering of SE signal vectors across five developmental stages, showing three distinct clusters. Each row represents an SE, characterized by its four representative enhancer units: strongest (maximum value), upper (second tertile), lower (first tertile), and weakest (minimum value), generated by a linear interpolation approach (see Methods for details). (B) Genome browser snapshots illustrating representative examples of Con, DN, and TS SEs. Tracks display H3K27ac CUT&Tag signals. (C) GO enrichment analysis for Con, DN, and TS SEs. (D) Average expression trends (relative to E60 stage) of genes associated with Con, DN, and TS SEs. Gene expression fold change (FC) values are log2 transformed. (E–G) Fold enrichment of Con, DN, and TS SEs for GWAS hits associated with pig meat and carcass traits (E), PSM cis-eQTLs (F), and mammalian conserved DNA elements (G). (H) Average DNA methylation levels of Con, DN, and TS SEs across five PSM stages. (I) Changes in the length of DNA nonmethylated regions within three types of SEs relative to the E60 stage. CpG sites are classified as methylated or nonmethylated, with contiguous nonmethylated sites merged into nonmethylated regions. See Supplemental Methods for details.

    Compared to DN and TS SEs, Con SEs exhibited significantly higher enrichment for GWAS hits associated with pig meat and carcass traits (Fig. 3E), PSM cis-eQTLs (Fig. 3F), and conserved DNA elements (Fig. 3G), underscoring their prominent role in PSM development. Simulation analysis confirmed that these enrichment patterns were not confounded by the size distribution of different SE classes (Supplemental Fig. S4D). Notably, Con SEs had the lowest average DNA methylation levels across all five developmental stages (Fig. 3H), supporting the transition of enhancers from basal to active states (Sharifi-Zarchi et al. 2017) and facilitating changes in cell identity (Guo et al. 2017). Given that average methylation levels might mask local differences within SEs, we further assessed DNA methylation dynamics based on the length of nonmethylated regions within SEs (see Supplemental Methods). Our results indicated a positive correlation between the length of nonmethylated DNA region (Fig. 3I), H3K27ac intensity (Fig. 3A), and expression levels of adjacent genes of SEs (Fig. 3D), indirectly validating our SE categorization. Taken together, we present multiple lines of evidence robustly supporting the prominent regulatory role of Con SEs in PSM development.

    Potential trait-affecting regulatory network mediated by candidate critical enhancers

    To pinpoint critical enhancers influencing muscle-relevant traits within SEs, we constructed a transcriptional regulatory network using an innovative analytical framework (Fig. 4A). We first identified 11,807 bidirectionally transcribed eRNAs (see Supplemental Methods), which exhibited positive correlations with H3K27ac intensity and negative correlations with DNA methylation (Supplemental Fig. S5A–C). Scanning SE eRNAs for the presence of the top 30 dynamic TF motifs (Fig. 2E) and combining TF-target gene coexpression with eRNA-target gene coexpression across 27 developmental stages (Supplemental Fig. S5D,E; Yang et al. 2019), we inferred the upstream TFs and downstream target genes of the enhancers corresponding to these SE eRNAs. We then prioritized the network by retaining only SE eRNAs with GWAS variants associated with pig meat and carcass traits.

    Figure 4.
    View larger version:
    Figure 4.

    Screening strategy for candidate critical enhancers. (A) Candidate critical enhancers were defined as SE units meeting three conditions: transcribing eRNAs, overlapping with GWAS hits related to pig meat and carcass traits (trait classes from Animal QTLdb), and harboring eRNA-located GWAS single-nucleotide polymorphisms (SNPs) that modulated enhancer activity in STARR-seq experiments. (B) Correlation of allelic ratios for SNPs between two STARR-seq biological replicates. (C,D) Identification of enhancer-modulating regulatory SNPs by the change in allelic ratios between the output library, including replicate1 (C) and replicate2 (D), and the input library. Statistical significance (false discovery rate [FDR] < 0.05) was evaluated using Fisher's exact test. (E) DNA-binding motif-breaking scores were assessed for regulatory SNPs and inactive SNPs on binding motifs of TFs expressed (TPM ≥ 3 at least at one stage) and nonexpressed (TPM < 3 at all stages) in PSM. (F) Number of enhancers retained during candidate critical enhancer screening process. (TAD) topologically associating domain.

    To exclude GWAS signals confounded by linkage disequilibrium (LD), we refined the network by filtering for GWAS SNPs that modulate enhancer activity using capture STARR-seq experiments. We designed a capture array targeting 300–400 bp DNA regions encompassing 610 SNPs within the retained SE eRNAs from eight eastern and eight western pigs, successfully capturing polymorphisms in 554 SNPs. The captured DNA fragments were cloned into the hSTARR-seq_ORI vector to construct input libraries, which were then transfected into mouse C2C12 myoblast cells. Fragment abundance was quantified in the output RNA-seq libraries, demonstrating robust screening performance with a high correlation in allelic ratios across biological replicates (Pearson's r = 0.98) (Fig. 4B). At FDR < 0.01 and |log2FC| > 0.58, we identified 55 regulatory SNPs affecting enhancer activity in both replicates (Fig. 4C,D; Supplemental Table S4). These SNPs showed significantly higher DNA-binding motif-breaking scores of TFs expressed in PSM compared to SNPs with no allelic enhancer activity (Fig. 4E; Supplemental Table S5). This difference was absent when considering motifs of TFs not expressed in PSM, highlighting the biological relevance (Fig. 4E). We also observed a fivefold enrichment of PSM cis-eQTLs among these regulatory SNPs (0.31 = 17/55) (Hypergeometric distribution test, P < 6.26 × 10−9), emphasizing their cis-regulatory potential in gene expression. We defined the enhancers containing these regulatory SNPs as candidate critical enhancers.

    We further refined the network by filtering cis-regulatory interactions based on PSM topologically associating domains (TADs), retaining only interactions where the candidate critical enhancer and its potential target gene were within the same TAD (Zhao et al. 2021). The final network consisted of 30 candidate critical enhancers (20, 7, and 3 from Con, DN, and TS SEs, respectively) (Fig. 4F), 23 upstream TFs, and 85 downstream target genes (Fig. 5; Supplemental Table S6). Notably, nearly all regulatory SNPs (16/17) overlapping with PSM cis-eQTLs (Teng et al. 2024) resided in Con candidate critical enhancers (Supplemental Table S7). Among them, four eGenes associated with seven cis-eQTLs have been identified in our network as the target genes of these cis-eQTLs. The effects of these cis-eQTLs on eGenes and their regulatory impact on enhancers are consistent. For example, the cis-eQTL Chr12:34231267_C/T and eGene CUEDC1 showed a negative association (β = −0.19), where the T allele corresponded to reduced CUEDC1 expression. In our network, this cis-eQTL was located in an enhancer (Chr12:34230720–34231974) predicted to regulate CUEDC1, with the T allele consistently downregulating enhancer activity. These findings suggest the reliability of our network, as the observed enhancer-modulating activity aligns with eQTL associations. Overall, our approach identified candidate critical enhancers within SEs, particularly Con SEs, which likely play a pivotal role in mediating transcriptional regulatory networks related to PSM development and meat traits.

    Figure 5.
    View larger version:
    Figure 5.

    Potential trait-affecting regulatory network mediated by candidate critical enhancers. Lavender nodes represent TFs, pink nodes represent candidate critical enhancers, and light blue nodes represent target genes. The size of each node is proportional to its degree, signifying that larger nodes have more connections. Regulatory interactions flow from TFs to enhancers to target genes. The initial network was constructed by integrating TF motif scanning, TF-target gene coexpression, and eRNA-target gene coexpression. The final network was refined through a series of filtering steps: SE identification, eRNA identification, GWAS hit association, regulatory SNP identification, and TAD filtering.

    A candidate critical enhancer is required for myogenic differentiation

    Due to the strong association of Con SEs with muscle cell differentiation and pork traits, we hypothesized that deleting Con candidate critical enhancers using CRISPR–Cas9 would impact myogenic differentiation. To test this, we evaluated the conservation of these enhancers at three levels compared to the mouse genome (Fig. 6A) and conducted functional validation in C2C12 cells, a robust and well-established myogenesis model. Four enhancers met the following conservation criteria (Fig. 6B): having orthologous sequences in the mouse genome (sequence conservation), orthologous sequences overlapping H3K27ac signals from either of the six myogenic differentiation stages of C2C12 cells (Peng et al. 2017) (epigenomics conservation) (Fig. 6C; Supplemental Fig. S6A), and sharing the same adjacent target genes (synteny conservation).

    Figure 6.
    View larger version:
    Figure 6.

    Conservation analysis and functional assay of candidate critical enhancers within Con SEs. (A) Schematic diagram illustrating the mapping of candidate critical enhancers to functional orthologous sequences in the mouse genome. (B) Flowchart illustrating the conservative filtering steps identifying functional orthologous sequences of candidate critical enhancers within Con SEs. The Arabic numerals indicate the number of candidate critical enhancers retained in each filtering step. (C) Genome browser snapshot displaying the mouse orthologous sequence of the candidate critical enhancer Chr12:34230720–34231974, extended 500 bp upstream and downstream. The tracks represent H3K27ac ChIP-seq signals across six stages during C2C12 myogenesis. The 24 h before induction of differentiation are denoted by “pre24hr,” and so forth. The gRNAs targeting the sequence and the predicted target gene were annotated. (D) qPCR analysis of the relative expression of four myogenic markers (Myod1, Myog, Mef2c, and Myhc) and one predicted target gene (Akap1) in wild-type (WT) and knockout (KO) C2C12 cells at four myogenic differentiation time points. KO cells were generated by deleting the orthologous sequence using CRISPR–Cas9. All qPCR data are presented as means ± SD of triplicate experiments. (E) Immunofluorescence staining of WT and KO C2C12 cells after 6 days of differentiation induction. Differentiated myotubes were labeled by immunostaining for MYOD1 (green) and MYHC (red), with nuclei stained by DAPI (blue). (F) qPCR analysis of the relative expression of myogenic markers in C2C12 cells with Akap1 knockdown on day 7 of myogenic differentiation. (G) Immunofluorescence staining for MYOD1 and MYHC was used to assess the effect of Akap1 knockdown on C2C12 differentiation.

    Among these conserved enhancers, the one at Chr12:34230720–34231974 exhibited the highest collinear conservation, sharing the most adjacent target genes with its mouse ortholog. This enhancer is of particular interest as it contains rs81434214, a GWAS SNP associated with the obesity index (OI), a selection index in pig breeding that reflects genetic predispositions toward obesity or leanness (Kogelman et al. 2014). Lean-type pig breeds, such as western commercial Landrace and Yorkshire, exhibit higher myogenic potential compared to obese-type breeds like Chinese domestic Tongcheng and Meishan (Zhao et al. 2011, 2015; Qiu et al. 2020). This enhancer may thus play a crucial role in muscle mass and quality differences between these pig breeds.

    To investigate the function of this enhancer, we performed a luciferase reporter assay, confirming that the A allele of its regulatory SNP (Chr12:34231458_T/A) reduced enhancer activity (Supplemental Fig. S6B,C), consistent with STARR-seq findings. CRISPR–Cas9 knockout of the orthologous sequence of this enhancer in C2C12 cells (Fig. 6C; Supplemental Fig. S7) led to reduced expression of myogenic marker (Myod1, Myog, Myhc) and predicted target gene Akap1 during myogenic differentiation (Fig. 6D), along with impaired differentiation (Fig. 6E). Consistently, siRNA knockdown of Akap1 (Supplemental Fig. S6D) also hindered C2C12 differentiation (Fig. 6F,G). These results serve as a proof of concept for our integrative strategy for identifying candidate trait-affecting enhancers and mapping their associated regulatory network.

    Previous SectionNext Section

    Discussion

    In this study, we identified dynamic SE trajectories closely linked to PSM development, introduced an innovative approach to pinpoint key enhancers driving the formation of complex traits, and applied this method to investigate the genetic basis of skeletal muscle development. To our knowledge, this study also represents the first comprehensive skeletal muscle SE landscape covering multidevelopmental stages before and after birth in any species, providing a valuable resource for skeletal muscle biology research.

    Leveraging H3K27ac data from an extended developmental period, we partitioned SEs based on their temporal patterns and identified Con SEs as crucial regulators for PSM development. In mammals, prenatal myofiber formation determines postnatal myofiber number and plays a decisive role in skeletal muscle development (Ashmore et al. 1973; Picard et al. 2002). Consequently, the noncoding cis-regulatory regions of prenatal PSM were prioritized for profiling and investigation (Yue et al. 2021; Zhao et al. 2021). However, our results indicate that Con SEs, which persist in both prenatal and postnatal stages, are more closely associated with PSM development. These Con SEs are characterized by stable mother units and highly dynamic constituent units throughput PSM development. The mother units serve as the foundational precursors and cornerstones of SE architecture (Bojcsuk et al. 2017), potentially providing a primed platform for the dynamic constituent units to fine-tune gene expression during PSM development. Moreover, the low methylation levels of Con SEs ensure high accessibility of these constituent units to TFs (Moore et al. 2013), which may facilitate a rapid response to external stimuli, such as nutrition or hormones, ensuring timely and adaptive skeletal muscle development.

    Our integrative approach leveraged Con SEs as regulatory hubs to identify critical enhancers affecting PSM development. A key challenge in enhancer profiling is pinpointing enhancers that truly influence phenotypes, as H3K27ac signals alone often mark both essential and nonessential enhancers. To overcome this, we employed an integrative approach that combined SE annotations with multiomics data. This strategy mirrors fine-mapping causal variants in GWAS studies and serves as an initial step to prioritize functional enhancers within SEs. SEs, known for regulating master transcriptional regulators during development, provided a fertile ground for identifying key constituent enhancers. Building on this foundation, we introduced Con SE annotation, which demonstrated a strong association with skeletal muscle GO terms and GWAS signals related to pork traits. To enhance functional annotation, we incorporated allele-specific STARR-seq signal analysis, aimed at filtering potential critical enhancers by evaluating the impacts of GWAS SNPs on enhancer activity. Together, this integrative approach offers an opportunity to identify critical enhancers supported by extensive evidence from epigenomic profiling, high-throughput enhancer screening, and population genetics.

    The effectiveness of our integrative approach was highlighted by experimental validations demonstrating the regulatory impact of a candidate critical enhancer (Chr12:34230720–34231974) and its predicted target gene (AKAP1) in muscle differentiation. AKAP1 regulates mitochondrial respiration via AMPK-dependent phosphorylation (Hoffman et al. 2015), and its deficiency results in mitochondrial dysfunction (Soundararajan et al. 2022). Mitochondria function as cellular energy factories, crucial for metabolic regulation and ATP production, supporting skeletal muscle growth and development (Peterson et al. 2012; Abrigo et al. 2019). Additionally, mitochondrial function correlates with muscle fiber type and antioxidant capacity, influencing meat tenderness and flavor (Zhang et al. 2022). Hence, we hypothesized that AKAP1 might be a causal gene affecting pork yield and quality. Through the OI GWAS SNP within this enhancer, our study advances understanding of the genetic basis of pork production traits.

    The recent publication of PigGTEx research provided an unprecedented amount of cis-eQTL data across various pig tissues (Teng et al. 2024), enabling a comprehensive interpretation of cis-regulatory effects underlying pig GWAS data. Compared to humans, skeletal muscle enhancers in pigs show less pronounced enrichment for skeletal muscle-specific cis-eQTLs (Supplemental Fig. S8A,B). This discrepancy may result from more precise identification of tissue-specific eQTLs in humans and broader LD regions in pigs. Given these challenges, high-throughput functional validation of enhancer-modulating SNPs, as demonstrated in this study, could effectively fine-map potentially causal variants within eQTL loci. Our findings indicate that 31% of identified enhancer-modulating SNPs correlate with PigGTEx skeletal muscle cis-eQTLs, primarily located in Con candidate critical enhancers, which are deemed essential for PSM development and pork traits. Notably, four eGenes with these cis-eQTLs were identified as potential target genes of the corresponding cis-eQTLs in our regulatory network. For these cis-eQTLs, alleles associated with enhanced enhancer activity consistently showed increased gene expression. This consistency suggests that our approach could facilitate the fine-mapping of functional cis-eQTL variants and potentially act as a proxy in situations where limited sample sizes hinder robust eQTL detection.

    Our framework has several limitations and areas for improvements. We initially classified distal H3K27ac peaks as potential enhancers and proximal ones as promoters based on Ensembl TSS annotations (Ling et al. 2024; Wang et al. 2024), which may be incomplete for pigs. To address this, we used a lenient initial screen to retain potential enhancers, followed by STARR-seq for functional validation. Future improvements could involve better TSS annotations and integrating more histone marks. We also refined our regulatory network using published TADs, which assign distal regulatory elements to target genes (Choi et al. 2021; Zuin et al. 2022; Wang et al. 2024). For greater precision, capture Hi-C could improve enhancer–gene mapping, though it is resource-intensive and has limitations (Roberts et al. 2024). Hi-C offers a static view of chromatin, so applying it to PSM across developmental stages would provide deeper insights. Additionally, our approach focuses on enhancer variants and may miss GWAS signals from coding mutations or other cis-regulatory regions like promoters, UTRs, silencers, insulators, or splicing sites. Nonetheless, with a few adjustments, integrating such GWAS signals with high-throughput allele-aware reporter assays targeting other types of cis-regulatory variants could significantly complement the enhancer-centric approach outlined in this study.

    Previous SectionNext Section

    Methods

    Preparation and sequencing of CUT&Tag libraries

    Skeletal muscle (longissimus dorsi) from Landrace pigs was sourced from our previous report (Yang et al. 2019, 2021a,b). The nuclear extraction protocol for frozen tissue was adapted from established methods with minor modifications (Nègre et al. 2006; Corces et al. 2017). Isolated nuclei were used to prepare the H3K27ac-modified library with a commercial CUT&Tag Kit (Vazyme TD901). Each experimental system required ∼50,000 nuclei, 3 μg of primary antibody (Abcam ab4729), and 1 μg of secondary antibody (Abcam ab6702). CUT&Tag libraries were sequenced on the NovaSeq 6000 platform (Illumina), generating 150 bp paired-end reads.

    CUT&Tag data analysis

    Raw reads from CUT&Tag libraries were filtered using Trim Galore! (v0.6.7) (https://github.com/FelixKrueger/TrimGalore), then aligned to the pig (S. scrofa) genome (Ensembl release 112) using Bowtie 2 (v2.4.2) (Langmead and Salzberg 2012). High-quality aligned and deduplicated reads were processed with SAMtools (v0.1.19) (Li et al. 2009) for H3K27ac peak calling via MACS2 (v2.2.7.1) (Zhang et al. 2008). Peaks with high intragroup repeatability were identified using IDR (v2.0.3) (Li et al. 2011). Initially, peaks were called from the merged BAM file with MACS2, using a loose threshold with a cutoff P of 0.01. Next, peaks with high coincidence with IDR peaks of each stage were retained using the BEDTools intersect command (v2.30.0, -F 0.5 -u) (Quinlan and Hall 2010), resulting in final merged peaks. Reads within these peaks were counted with the multiBamSummary function of BEDTools, and depth was normalized by calculating TPM. Finally, edgeR (v3.40.0) (Robinson et al. 2010) of the R software (R Core Team 2023) was used to correct library size and obtain the final peak signal intensity matrix. For motif enrichment analysis, genomic features within 100 bp of peak summits were analyzed using the HOMER findMotifsGenome.pl program (v4.11) (Heinz et al. 2010).

    Classification of the temporal patterns of SEs

    SEs consist of various constituent units, each with different quantities and positions. To cluster SEs, a previously established method was employed, creating four artificial constituent units for each SE (Kai et al. 2021). These artificial units were based on the zeroth, 33.3th, 66.6th, and 100th percentiles of the H3K27ac intensity range of each stage's true constituent units. H3K27ac intensities of all artificial units across all five developmental stages were combined to generate a 20-column signal matrix for SEs. The pheatmap (v1.0.12) R package (https://CRAN.R-project.org/package=pheatmap) was used for hierarchical clustering and visualization of this matrix.

    Construction of a TFs/eRNAs/targets regulatory network

    Seqmonk (v1.47.2) was used to detect bidirectional transcription of eRNAs within ±500 bp of enhancer regions and quantify eRNA expression levels (see Supplemental Methods), following a previous study (Carullo et al. 2020). An initial TFs/eRNAs/targets network was constructed based on the top 30 dynamic TF motifs (Fig. 2E) in SE eRNAs using the HOMER scanMotifGenomeWide.pl tool, coexpression of TFs and genes (Spearman's correlation Rs ≥ 0.3, FDR < 0.05), and coexpression of eRNAs and genes (distance ≤1 Mb, Spearman's correlation Rs ≥ 0.3, FDR < 0.05). The corr.test function in the psych R package (v2.2.5) was used to calculate correlation values and statistical significance. The final trait-affecting network was then refined by retaining only SE eRNAs that overlapped (within ±20 kb) with GWAS hits associated with pig meat and carcass traits (Hu et al. 2019) and harbored regulatory SNPs that modulate enhancer activity.

    Identification of regulatory SNPs by STARR-seq

    The enhancer-modulating capacity of SNPs in the initial network was evaluated by allele-specific STARR-seq strategy referring to our previous work (Liu et al. 2017a,b; Wang et al. 2024). Briefly, genomic regions containing candidate PigVar SNPs (Zhou et al. 2017) were isolated from a pooled genomic DNA sample from eight eastern and eight western pigs using the Enzyme Plus Library Prep Kit (iGeneTech Biotech C10011) and TargetSeq One Hyb and Wash Kit v2.0 (iGeneTech Biotech C10331), with fragment lengths of 100–200 bp. Captured fragments were ligated into the hSTARR-seq_ORI vector (Addgene 99296). Recombinant vectors were transformed into Escherichia coli DH10B cells (Life Technologies) using a Gene Pulser Xcell Electroporation System (Bio-Rad Laboratories) to construct an input plasmid library, which was then electroporated into C2C12 cells using an Amaxa 4D-Nucleofector (Lonza). Cells were harvested after 24 h. Total RNA was extracted, and mRNA was enriched and reverse transcribed. cDNA was amplified by PCR to generate an output library. High-throughput sequencing on the input and output libraries was performed using the PE150 strategy on the Illumina platform. Allelic ratio fold change (output/input) was calculated as the effect size of each SNP, and two-tailed Fisher's exact tests determined the significance of allele ratio differences.

    Generation of enhancer-deleted C2C12 cells by CRISPR–Cas9

    pgRNAs targeting the 5′- and 3′-flanking regions (500 bp upstream and 500 bp downstream) of the targeted regulatory locus were designed using the CRISPOR web tool (Concordet and Haeussler 2018) and synthesized by Tsingke Biotech Co., Ltd. for CRISPR–Cas9-mediated deletion. The pgRNA sequences are listed in Supplemental Figure S7. The pgRNAs were assembled into the lentiviral vector pLV[2gRNA]-Puro-U6>Scramble[gRNA#1]-U6>Scramble[gRNA#2] (VectorBuilder VB220224-1316qur) by Tsingke Biotech Co., Ltd using a seamless cloning approach. Lentivirus was produced by transfecting the lentiviral vector and packaging plasmids into HEK293T cells using jetPRIME reagent (Polyplus). C2C12 myoblasts with stable Cas9 expression were infected with the lentivirus. Positive cells after puromycin selection were seeded into 96-well plates by dilution plating and grown into single-cell colonies. Cell clones with the deletion of the target region were verified by PCR and sequencing.

    Investigation of the enhancer function in myogenesis

    The C2C12 cells were switched to a differentiation medium consisting of 2% horse serum (Gibco) and 1% penicillin–streptomycin in DMEM, followed by a 5-day myogenic induction. Immunofluorescence was conducted using anti-MYHC antibody (DHSB MF20), anti-MYOD1 antibody (Proteintech 18943-1-AP), anti-rabbit FITC antibody (Servicebio GB22301), anti-mouse CY3 antibody (Servicebio GB21303), and DAPI for nuclei labeling. The stained cells were observed under an inverted fluorescence microscope. Total RNA was extracted from the collected cells using M5 Universal Plus RNA Mini Kit (Mei5 Biotechnology MF036) and reverse transcribed using HiScript III RT SuperMix for qPCR (Vazyme R323-01). qPCR was performed using ChamQ Universal SYBR qPCR Master Mix (Vazyme Q341-02) in a QuantStudio3 system (Applied Biosystems). qPCR primers are listed in Supplemental Table S8.

    Investigation of the role of Akap1 in myogenesis

    The siRNAs targeting Akap1 (5′-TGCACAGAGACCATCTCAA-3′) and NC were custom synthesized by Ribo Biological Co., Ltd. C2C12 cells were transfected with the siRNAs using the jetPRIME reagent (Polyplus) when they reached 80% confluency. After 8 h, the differentiation medium was added to induce myogenic differentiation. An additional siRNA transfection was performed after 3 days, continuing the differentiation process for another 3 days. The degree of differentiation was evaluated using immunofluorescence and qPCR analyses.

    Previous SectionNext Section

    Data access

    The raw CUT&Tag and STARR-seq data generated in this work have been submitted to the Genome Sequence Archive (GSA; https://ngdc.cncb.ac.cn/gsa/) under accession number CRA011158.

    Previous SectionNext Section

    Competing interest statement

    The authors declare no competing interests.

    Previous SectionNext Section

    Acknowledgments

    This work was supported by the China National Key R&D Program during the 14th Five-year Plan Period China (2021YFF1200500 to Y.L.), the National Natural Science Foundation of China (32070595 and 32372858 to Y.L.), and the Shenzhen Science and Technology Innovation Commission (JCYJ20190814163803664 to Y.L.).

    Author contributions: Y.L. and S.Z. designed this study. S.Z. worked on the preparation and sequencing of CUT&Tag libraries. S.Z., C.W., Q.L., and Y.Zhao analyzed the data. S.Q., C.C., Y.B., Y.Zhang, S.Z., and L.X. performed the functional validation experiments. Y.L. supervised the study. S.Z. wrote the manuscript. Y.L. guided the writing and revised the manuscript. Z.T. and K.L. provided samples and revised the manuscript. All authors provided feedback during the drafting process and approved the final version.

    Previous SectionNext Section

    Footnotes

    • Received July 29, 2023.
    • Accepted October 15, 2024.

    This article is distributed exclusively by Cold Spring Harbor Laboratory Press for the first six months after the full-issue publication date (see https://genome.cshlp.org/site/misc/terms.xhtml). After six months, it is available under a Creative Commons License (Attribution-NonCommercial 4.0 International), as described at http://creativecommons.org/licenses/by-nc/4.0/.

    Previous Section
     

    References

    1. ↵
      Abrigo J, Simon F, Cabrera D, Vilos C, Cabello-Verrugio C. 2019. Mitochondrial dysfunction in skeletal muscle pathologies. Curr Protein Pept Sci 20: 536–546. doi:10.2174/1389203720666190402100902
      CrossRefMedlineGoogle Scholar
    2. ↵
      Aneas I, Decker DC, Howard CL, Sobreira DR, Sakabe NJ, Blaine KM, Stein MM, Hrusch CL, Montefiori LE, Tena J, et al. 2021. Asthma-associated genetic variants induce IL33 differential expression through an enhancer-blocking regulatory region. Nat Commun 12: 6115. doi:10.1038/s41467-021-26347-z
      CrossRefMedlineGoogle Scholar
    3. ↵
      Ashmore CR, Addis PB, Doerr L. 1973. Development of muscle fibers in the fetal pig. J Anim Sci 36: 1088–1093. doi:10.2527/jas1973.3661088x
      CrossRefMedlineGoogle Scholar
    4. ↵
      Beermann DH, Cassens RG. 1977. Indirect fluorescence of primary and secondary myofibers in developing porcine muscle. J Histochem Cytochem 25: 439–442. doi:10.1177/25.6.881553
      CrossRefMedlineGoogle Scholar
    5. ↵
      Bi P, Kuang S. 2012. Meat Science and Muscle Biology Symposium: stem cell niche and postnatal muscle growth. J Anim Sci 90: 924–935. doi:10.2527/jas.2011-4594
      CrossRefMedlineGoogle Scholar
    6. ↵
      Bojcsuk D, Nagy G, Balint BL. 2017. Inducible super-enhancers are organized based on canonical signal-specific transcription factor binding elements. Nucleic Acids Res 45: 3693–3706. doi:10.1093/nar/gkw1283
      CrossRefMedlineGoogle Scholar
    7. ↵
      Brameld JM. 1997. Molecular mechanisms involved in the nutritional and hormonal regulation of growth in pigs. Proc Nutr Soc 56: 607–619. doi:10.1079/pns19970061
      CrossRefMedlineGoogle Scholar
    8. ↵
      Carullo NVN, Phillips RA III., Simon RC, Soto SAR, Hinds JE, Salisbury AJ, Revanna JS, Bunner KD, Ianov L, Sultan FA, et al. 2020. Enhancer RNAs predict enhancer–gene regulatory links and are critical for enhancer function in neuronal systems. Nucleic Acids Res 48: 9550–9570. doi:10.1093/nar/gkaa671
      CrossRefMedlineGoogle Scholar
    9. ↵
      Choi J, Lysakovskaia K, Stik G, Demel C, Söding J, Tian TV, Graf T, Cramer P. 2021. Evidence for additive and synergistic action of mammalian enhancers during cell fate determination. eLife 10: e65381. doi:10.7554/eLife.65381
      CrossRefMedlineGoogle Scholar
    10. ↵
      Concordet J-P, Haeussler M. 2018. CRISPOR: intuitive guide selection for CRISPR/Cas9 genome editing experiments and screens. Nucleic Acids Res 46: W242–W245. doi:10.1093/nar/gky354
      CrossRefMedlineGoogle Scholar
    11. ↵
      Corces MR, Trevino AE, Hamilton EG, Greenside PG, Sinnott-Armstrong NA, Vesuna S, Satpathy AT, Rubin AJ, Montine KS, Wu B, et al. 2017. An improved ATAC-seq protocol reduces background and enables interrogation of frozen tissues. Nat Methods 14: 959–962. doi:10.1038/nmeth.4396
      CrossRefMedlineGoogle Scholar
    12. ↵
      Creyghton MP, Cheng AW, Welstead GG, Kooistra T, Carey BW, Steine EJ, Hanna J, Lodato MA, Frampton GM, Sharp PA, et al. 2010. Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proc Natl Acad Sci 107: 21931–21936. doi:10.1073/pnas.1016071107
      Abstract/FREE Full Text
    13. ↵
      The ENCODE Project Consortium. 2012. An integrated encyclopedia of DNA elements in the human genome. Nature 489: 57–74. doi:10.1038/nature11247
      CrossRefMedlineGoogle Scholar
    14. ↵
      Ernst J, Bar-Joseph Z. 2006. STEM: a tool for the analysis of short time series gene expression data. BMC Bioinformatics 7: 191. doi:10.1186/1471-2105-7-191
      CrossRefMedlineGoogle Scholar
    15. ↵
      Giuffra E, Tuggle CK, the FAANG Consortium. 2019. Functional Annotation of Animal Genomes (FAANG): current achievements and roadmap. Annu Rev Anim Biosci 7: 65–88. doi:10.1146/annurev-animal-020518-114913
      CrossRefMedlineGoogle Scholar
    16. ↵
      Guo H, Hu B, Yan L, Yong J, Wu Y, Gao Y, Guo F, Hou YU, Fan X, Dong JI, et al. 2017. DNA methylation and chromatin accessibility profiling of mouse and human fetal germ cells. Cell Res 27: 165–183. doi:10.1038/cr.2016.128
      CrossRefMedlineGoogle Scholar
    17. ↵
      He H, Li W, Liyanarachchi S, Srinivas M, Wang Y, Akagi K, Wang Y, Wu D, Wang Q, Jin V, et al. 2015. Multiple functional variants in long-range enhancer elements contribute to the risk of SNP rs965513 in thyroid cancer. Proc Natl Acad Sci 112: 6128–6133. doi:10.1073/pnas.1506255112
      Abstract/FREE Full Text
    18. ↵
      Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, Cheng JX, Murre C, Singh H, Glass CK. 2010. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell 38: 576–589. doi:10.1016/j.molcel.2010.05.004
      CrossRefMedlineGoogle Scholar
    19. ↵
      Hernández-Hernández JM, García-González EG, Brun CE, Rudnicki MA. 2017. The myogenic regulatory factors, determinants of muscle development, cell identity and regeneration. Semin Cell Dev Biol 72: 10–18. doi:10.1016/j.semcdb.2017.11.010
      CrossRefMedlineGoogle Scholar
    20. ↵
      Hnisz D, Abraham BJ, Lee TI, Lau A, Saint-André V, Sigova AA, Hoke HA, Young RA. 2013. Super-enhancers in the control of cell identity and disease. Cell 155: 934–947. doi:10.1016/j.cell.2013.09.053
      CrossRefMedlineGoogle Scholar
    21. ↵
      Hoffman NJ, Parker BL, Chaudhuri R, Fisher-Wellman KH, Kleinert M, Humphrey SJ, Yang P, Holliday M, Trefely S, Fazakerley DJ, et al. 2015. Global phosphoproteomic analysis of human skeletal muscle reveals a network of exercise-regulated kinases and AMPK substrates. Cell Metab 22: 922–935. doi:10.1016/j.cmet.2015.09.001
      CrossRefMedlineGoogle Scholar
    22. ↵
      Hu Z-L, Park CA, Reecy JM. 2019. Building a livestock genetic and genomic information knowledgebase through integrative developments of Animal QTLdb and CorrDB. Nucleic Acids Res 47: D701–D710. doi:10.1093/nar/gky1084
      CrossRefMedlineGoogle Scholar
    23. ↵
      Kai Y, Li BE, Zhu M, Li GY, Chen F, Han Y, Cha HJ, Orkin SH, Cai W, Huang J, et al. 2021. Mapping the evolving landscape of super-enhancers during cell differentiation. Genome Biol 22: 269. doi:10.1186/s13059-021-02485-x
      CrossRefMedlineGoogle Scholar
    24. ↵
      Kaya-Okur HS, Wu SJ, Codomo CA, Pledger ES, Bryson TD, Henikoff JG, Ahmad K, Henikoff S. 2019. CUT&Tag for efficient epigenomic profiling of small samples and single cells. Nat Commun 10: 1930. doi:10.1038/s41467-019-09982-5
      CrossRefMedlineGoogle Scholar
    25. ↵
      Kern C, Wang Y, Xu X, Pan Z, Halstead M, Chanthavixay G, Saelao P, Waters S, Xiang R, Chamberlain A, et al. 2021. Functional annotations of three domestic animal genomes provide vital resources for comparative and agricultural research. Nat Commun 12: 1821. doi:10.1038/s41467-021-22100-8
      CrossRefMedlineGoogle Scholar
    26. ↵
      Kogelman LJA, Pant SD, Fredholm M, Kadarmideen HN. 2014. Systems genetics of obesity in an F2 pig model by genome-wide association, genetic network, and pathway analyses. Front Genet 5: 214. doi:10.3389/fgene.2014.00214
      CrossRefMedlineGoogle Scholar
    27. ↵
      Langmead B, Salzberg SL. 2012. Fast gapped-read alignment with Bowtie 2. Nat Methods 9: 357–359. doi:10.1038/nmeth.1923
      CrossRefMedlineGoogle Scholar
    28. ↵
      Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, 1000 Genome Project Data Processing Subgroup. 2009. The Sequence Alignment/Map format and SAMtools. Bioinformatics 25: 2078–2079. doi:10.1093/bioinformatics/btp352
      CrossRefMedlineGoogle Scholar
    29. ↵
      Li Q, Brown JB, Huang H, Bickel PJ. 2011. Measuring reproducibility of high-throughput experiments. Ann Appl Stat 5: 1752–1779. doi:10.1214/11-AOAS466
      CrossRefMedlineGoogle Scholar
    30. ↵
      Li J, Xiang Y, Zhang L, Qi X, Zheng Z, Zhou P, Tang Z, Jin Y, Zhao Q, Fu Y, et al. 2022. Enhancer-promoter interaction maps provide insights into skeletal muscle-related traits in pig genome. BMC Biol 20: 136. doi:10.1186/s12915-022-01322-2
      CrossRefMedlineGoogle Scholar
    31. ↵
      Ling Z, Li J, Jiang T, Zhang Z, Zhu Y, Zhou Z, Yang J, Tong X, Yang B, Huang L. 2024. Omics-based construction of regulatory variants can be applied to help decipher pig liver-related traits. Commun Biol 7: 381. doi:10.1038/s42003-024-06050-7
      CrossRefMedlineGoogle Scholar
    32. ↵
      Liu S, Liu Y, Zhang Q, Wu J, Liang J, Yu S, Wei G-H, White KP, Wang X. 2017a. Systematic identification of regulatory variants associated with cancer risk. Genome Biol 18: 194. doi:10.1186/s13059-017-1322-z
      CrossRefMedlineGoogle Scholar
    33. ↵
      Liu Y, Yu S, Dhiman VK, Brunetti T, Eckart H, White KP. 2017b. Functional assessment of human enhancer activities using whole-genome STARR-sequencing. Genome Biol 18: 219. doi:10.1186/s13059-017-1345-5
      CrossRefMedlineGoogle Scholar
    34. ↵
      Liu X, Zhang Y, Liu W, Li Y, Pan J, Pu Y, Han J, Orlando L, Ma Y, Jiang L. 2022. A single-nucleotide mutation within the TBX3 enhancer increased body size in Chinese horses. Curr Biol 32: 480–487.e6. doi:10.1016/j.cub.2021.11.052
      CrossRefMedlineGoogle Scholar
    35. ↵
      Lunney JK, Van Goor A, Walker KE, Hailstock T, Franklin J, Dai C. 2021. Importance of the pig as a human biomedical model. Sci Transl Med 13: eabd5758. doi:10.1126/scitranslmed.abd5758
      CrossRefMedlineGoogle Scholar
    36. ↵
      Moore LD, Le T, Fan G. 2013. DNA methylation and its basic function. Neuropsychopharmacology 38: 23–38. doi:10.1038/npp.2012.112
      CrossRefMedlineGoogle Scholar
    37. ↵
      Mukund K, Subramaniam S. 2020. Skeletal muscle: a review of molecular structure and function, in health and disease. Wiley Interdiscip Rev Syst Biol Med 12: e1462. doi:10.1002/wsbm.1462
      CrossRefMedlineGoogle Scholar
    38. ↵
      Muráni E, Murániová M, Ponsuksili S, Schellander K, Wimmers K. 2007. Identification of genes differentially expressed during prenatal development of skeletal muscle in two pig breeds differing in muscularity. BMC Dev Biol 7: 109. doi:10.1186/1471-213X-7-109
      CrossRefMedlineGoogle Scholar
    39. ↵
      Nègre N, Hennetin J, Sun LV, Lavrov S, Bellis M, White KP, Cavalli G. 2006. Chromosomal distribution of PcG proteins during Drosophila development. PLoS Biol 4: e170. doi:10.1371/journal.pbio.0040170
      CrossRefMedlineGoogle Scholar
    40. ↵
      Pan Z, Yao Y, Yin H, Cai Z, Wang Y, Bai L, Kern C, Halstead M, Chanthavixay G, Trakooljul N, et al. 2021. Pig genome functional annotation enhances the biological interpretation of complex traits and human disease. Nat Commun 12: 5848. doi:10.1038/s41467-021-26153-7
      CrossRefMedlineGoogle Scholar
    41. ↵
      Peng XL, So KK, He L, Zhao Y, Zhou J, Li Y, Yao M, Xu B, Zhang S, Yao H, et al. 2017. MyoD- and FoxO3-mediated hotspot interaction orchestrates super-enhancer activity during myogenic differentiation. Nucleic Acids Res 45: 8785–8805. doi:10.1093/nar/gkx488
      CrossRefMedlineGoogle Scholar
    42. ↵
      Peterson CM, Johannsen DL, Ravussin E. 2012. Skeletal muscle mitochondria and aging: a review. J Aging Res 2012: 194821. doi:10.1155/2012/194821
      CrossRefMedlineGoogle Scholar
    43. ↵
      Picard B, Lefaucheur L, Berri C, Duclos MJ. 2002. Muscle fibre ontogenesis in farm animal species. Reprod Nutr Dev 42: 415–431. doi:10.1051/rnd:2002035
      CrossRefMedlineGoogle Scholar
    44. ↵
      Qiu K, Xu D, Wang L, Zhang X, Jiao N, Gong L, Yin J. 2020. Association analysis of single-cell RNA sequencing and proteomics reveals a vital role of Ca2+ signaling in the determination of skeletal muscle development potential. Cells 9: 1045. doi:10.3390/cells9041045
      CrossRefGoogle Scholar
    45. ↵
      Quinlan AR, Hall IM. 2010. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26: 841–842. doi:10.1093/bioinformatics/btq033
      CrossRefMedlineGoogle Scholar
    46. ↵
      R Core Team. 2023. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna. https://www.R-project.org/.
      Google Scholar
    47. ↵
      Roadmap Epigenomics Consortium, Kundaje A, Meuleman W, Ernst J, Bilenky M, Yen A, Kheradpour P, Wang J, Whitaker JW, Schultz MD, et al. 2015. Integrative analysis of 111 reference human epigenomes. Nature 518: 317–330. doi:10.1038/nature14248
      CrossRefMedlineGoogle Scholar
    48. ↵
      Roberts BS, Anderson AG, Partridge EC, Cooper GM, Myers RM. 2024. Probabilistic association of differentially expressed genes with cis-regulatory elements. Genome Res 34: 620–632. doi:10.1101/gr.278598.123
      Abstract/FREE Full Text
    49. ↵
      Robinson MD, McCarthy DJ, Smyth GK. 2010. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26: 139–140. doi:10.1093/bioinformatics/btp616
      CrossRefMedlineGoogle Scholar
    50. ↵
      Scheffler TL, Scheffler JM, Park S, Kasten SC, Wu Y, McMillan RP, Hulver MW, Frisard MI, Gerrard DE. 2014. Fiber hypertrophy and increased oxidative capacity can occur simultaneously in pig glycolytic skeletal muscle. Am J Physiol Cell Physiol 306: C354–C363. doi:10.1152/ajpcell.00002.2013
      CrossRefMedlineGoogle Scholar
    51. ↵
      Sharifi-Zarchi A, Gerovska D, Adachi K, Totonchi M, Pezeshk H, Taft RJ, Schöler HR, Chitsaz H, Sadeghi M, Baharvand H, et al. 2017. DNA methylation regulates discrimination of enhancers from promoters through a H3K4me1-H3K4me3 seesaw mechanism. BMC Genomics 18: 964. doi:10.1186/s12864-017-4353-7
      CrossRefMedlineGoogle Scholar
    52. ↵
      Shen Y, Yue F, McCleary DF, Ye Z, Edsall L, Kuan S, Wagner U, Dixon J, Lee L, Lobanenkov VV, et al. 2012. A map of the cis-regulatory sequences in the mouse genome. Nature 488: 116–120. doi:10.1038/nature11243
      CrossRefMedlineGoogle Scholar
    53. ↵
      Soundararajan R, Hernández-Cuervo H, Stearns TM, Griswold AJ, Patil SS, Fukumoto J, Narala VR, Galam L, Lockey R, Kolliputi N. 2022. A-kinase Anchor Protein 1 deficiency causes mitochondrial dysfunction in mouse model of hyperoxia induced acute lung injury. Front Pharmacol 13: 980723. doi:10.3389/fphar.2022.980723
      CrossRefMedlineGoogle Scholar
    54. ↵
      Teng J, Gao Y, Yin H, Bai Z, Liu S, Zeng H, the PigGTEx Consortium, Bai L, Cai Z, Zhao B, et al. 2024. A compendium of genetic regulatory effects across pig tissues. Nat Genet 56: 112–123. doi:10.1038/s41588-023-01585-7
      CrossRefMedlineGoogle Scholar
    55. ↵
      Van Erp-Van der Kooij E, Kuijpers AH, Van Eerdenburg F, Tielen MJM. 2003. Coping characteristics and performance in fattening pigs. Livest Prod Sci 84: 31–38. doi:10.1016/S0301-6226(03)00072-1
      CrossRefGoogle Scholar
    56. ↵
      Wang C, Chen C, Lei B, Qin S, Zhang Y, Li K, Zhang S, Liu Y. 2024. Constructing eRNA-mediated gene regulatory networks to explore the genetic basis of muscle and fat-relevant traits in pigs. Genet Sel Evol 56: 28. doi:10.1186/s12711-024-00897-4
      CrossRefMedlineGoogle Scholar
    57. ↵
      Whyte WA, Orlando DA, Hnisz D, Abraham BJ, Lin CY, Kagey MH, Rahl PB, Lee TI, Young RA. 2013. Master transcription factors and mediator establish super-enhancers at key cell identity genes. Cell 153: 307–319. doi:10.1016/j.cell.2013.03.035
      CrossRefMedlineGoogle Scholar
    58. ↵
      Woolfe A, Goodson M, Goode DK, Snell P, McEwen GK, Vavouri T, Smith SF, North P, Callaway H, Kelly K, et al. 2005. Highly conserved non-coding sequences are associated with vertebrate development. PLoS Biol 3: e7. doi:10.1371/journal.pbio.0030007
      CrossRefMedlineGoogle Scholar
    59. ↵
      Yan W, Chen D, Schumacher J, Durantini D, Engelhorn J, Chen M, Carles CC, Kaufmann K. 2019. Dynamic control of enhancer activity drives stage-specific gene expression during flower morphogenesis. Nat Commun 10: 1705. doi:10.1038/s41467-019-09513-2
      CrossRefGoogle Scholar
    60. ↵
      Yang Y, Zhu M, Fan X, Yao Y, Yan J, Tang Y, Liu S, Li K, Tang Z. 2019. Developmental atlas of the RNA editome in Sus scrofa skeletal muscle. DNA Res 26: 261–272. doi:10.1093/dnares/dsz006
      CrossRefMedlineGoogle Scholar
    61. ↵
      Yang Y, Yan J, Fan X, Chen J, Wang Z, Liu X, Yi G, Liu Y, Niu Y, Zhang L, et al. 2021a. The genome variation and developmental transcriptome maps reveal genetic differentiation of skeletal muscle in pigs. PLoS Genet 17: e1009910. doi:10.1371/journal.pgen.1009910
      CrossRefMedlineGoogle Scholar
    62. ↵
      Yang Y, Fan X, Yan J, Chen M, Zhu M, Tang Y, Liu S, Tang Z. 2021b. A comprehensive epigenome atlas reveals DNA methylation regulating skeletal muscle development. Nucleic Acids Res 49: 1313–1329. doi:10.1093/nar/gkaa1203
      CrossRefMedlineGoogle Scholar
    63. ↵
      Yuan R, Zhang J, Wang Y, Zhu X, Hu S, Zeng J, Liang F, Tang Q, Chen Y, Chen L, et al. 2021. Reorganization of chromatin architecture during prenatal development of porcine skeletal muscle. DNA Res 28: dsab003. doi:10.1093/dnares/dsab003
      CrossRefMedlineGoogle Scholar
    64. ↵
      Yue F, Cheng Y, Breschi A, Vierstra J, Wu W, Ryba T, Sandstrom R, Ma Z, Davis C, Pope BD, et al. 2014. A comparative encyclopedia of DNA elements in the mouse genome. Nature 515: 355–364. doi:10.1038/nature13992
      CrossRefMedlineGoogle Scholar
    65. ↵
      Yue J, Hou X, Liu X, Wang L, Gao H, Zhao F, Shi L, Shi L, Yan H, Deng T, et al. 2021. The landscape of chromatin accessibility in skeletal muscle during embryonic development in pigs. J Anim Sci Biotechnol 12: 56. doi:10.1186/s40104-021-00577-z
      CrossRefGoogle Scholar
    66. ↵
      Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, Nusbaum C, Myers RM, Brown M, Li W, et al. 2008. Model-based analysis of ChIP-Seq (MACS). Genome Biol 9: R137. doi:10.1186/gb-2008-9-9-r137
      CrossRefMedlineGoogle Scholar
    67. ↵
      Zhang X, Chen M, Yan E, Wang Y, Ma C, Zhang P, Yin J. 2022. Dietary malic acid supplementation induces skeletal muscle fiber-type transition of weaned piglets and further improves meat quality of finishing pigs. Front Nutr 8: 825495. doi:10.3389/fnut.2021.825495
      CrossRefMedlineGoogle Scholar
    68. ↵
      Zhao X, Mo D, Li A, Gong W, Xiao S, Zhang Y, Qin L, Niu Y, Guo Y, Liu X, et al. 2011. Comparative analyses by sequencing of transcriptomes during skeletal muscle development between pig breeds differing in muscle growth rate and fatness. PLoS One 6: e19774. doi:10.1371/journal.pone.0019774
      CrossRefMedlineGoogle Scholar
    69. ↵
      Zhao Y, Li J, Liu H, Xi Y, Xue M, Liu W, Zhuang Z, Lei M. 2015. Dynamic transcriptome profiles of skeletal muscle tissue across 11 developmental stages for both Tongcheng and Yorkshire pigs. BMC Genomics 16: 377. doi:10.1186/s12864-015-1580-7
      CrossRefMedlineGoogle Scholar
    70. ↵
      Zhao Y, Hou Y, Xu Y, Luan Y, Zhou H, Qi X, Hu M, Wang D, Wang Z, Fu Y, et al. 2021. A compendium and comparative epigenomics analysis of cis-regulatory elements in the pig genome. Nat Commun 12: 2217. doi:10.1038/s41467-021-22448-x
      CrossRefMedlineGoogle Scholar
    71. ↵
      Zhou Z-Y, Li A, Otecko NO, Liu Y-H, Irwin DM, Wang L, Adeola AC, Zhang J, Xie H-B, Zhang Y-P. 2017. PigVar: a database of pig variations and positive selection signatures. Database 2017: bax048. doi:10.1093/database/bax048
      CrossRefMedlineGoogle Scholar
    72. ↵
      Zhu Y, Zhou Z, Huang T, Zhang Z, Li W, Ling Z, Jiang T, Yang J, Yang S, Xiao Y. 2022. Mapping and analysis of a spatiotemporal H3K27ac and gene expression spectrum in pigs. Sci China Life Sci 65: 1517–1534. doi:10.1007/s11427-021-2034-5
      CrossRefMedlineGoogle Scholar
    73. ↵
      Zuin J, Roth G, Zhan Y, Cramard J, Redolfi J, Piskadlo E, Mach P, Kryzhanovska M, Tihanyi G, Kohler H, et al. 2022. Nonlinear control of transcription through enhancer–promoter interactions. Nature 604: 571–577. doi:10.1038/s41586-022-04570-y
      CrossRefGoogle Scholar

    This article has not yet been cited by other articles.

    | Table of Contents

    This Article

    1. Published in Advance October 21, 2024, doi: 10.1101/gr.278344.123 Genome Res. 34: 2190-2202 © 2024 Zhang et al.; Published by Cold Spring Harbor Laboratory Press

    Article Category

    Share

    Preprint Server



    Navigate This Article

    Current Issue

    1. January 2026, 36 (1)
    1. Current Issue
    1. Alert me to new issues of Genome Research