Skip to main content
Download PDF

Genome-wide chromatin accessibility and selective signals of meat rabbits reveal key Cis-regulatory elements and variants during postnatal development of skeletal muscles in rabbits

BMC Genomics volume 26, Article number: 296 (2025) Cite this article

Abstract

Background

The development of skeletal muscles is intricately modulated by multiple genetic factors and significantly impacts the economic value of meat rabbits. However, our knowledge of epigenetics in rabbit skeletal muscles remains largely unknown.

Results

In this study, we collected leg skeletal muscles of rabbits and performed assays for transposase-accessible chromatin with high throughput sequencing (ATAC-seq) to detect open chromatin across three developmental stages: birth (D1), weaning (D35), and adulthood (D75). A total of 126,959 accessible chromatin regions (ACRs) were identified across samples, and a broad increase and decrease in chromatin accessibility were found from D1 to D35 and D35 to D75, respectively. Integrative analysis of chromatin accessibility and transcriptome data revealed ACRs that were nearly closed at D1 but highly accessible at D35 and D75 were significantly enriched in skeletal muscle development. Cis-regulation analysis further revealed that genes dominated by enhancers mainly play roles in the neuron development of rabbit skeletal muscles. Moreover, the detection of selection signals of meat rabbits and the footprinting analysis of transcription factor at open chromatin revealed that both base transversion (Chr13:12144967 A-> G) and the dynamics of chromatin accessibility at the PRDM1 binding site might regulate ZSWIM5 during the development of skeletal muscles in rabbits.

Conclusions

Our study provided a category of potential cis-regulatory elements for understanding the development of skeletal muscles at the tissue level and might facilitate potential insights into growth regulation in rabbits.

Peer Review reports

Introduction

With the increase in demand for healthy food in humans, rabbit (Oryctolagus cuniculus) meat has become increasingly popular due to the rich protein, low cholesterol levels, and low fat in the products [1]. The growth and meat quality of skeletal muscles are the core traits for rabbit breeding, especially in improving specialized meat rabbit strains [2]. Skeletal muscle development is generally divided into prenatal and postnatal development [3, 4]. During the prenatal stage, skeletal muscle mainly increases the number of muscle fibers in mammals [5]. During the postnatal stage, mammals experience an increase in the diameter and length of muscle fibers, as well as the acquisition of motor function and the deposition of flavor substances, which is a complex biological process [6]. Previous studies have revealed a sort of genetic factors regulating the prenatal development of skeletal muscles at the post-transcriptional layer of regulatory networks, such as microRNAs [7], lncRNAs [8], and circular RNAs [9]. At the genomic level, genome-wide association studies have revealed significant genetic variants affecting the growth trait [10]. Nevertheless, compared with other common livestock species, the genomic state during skeletal muscle development remains largely unknown in rabbits.

Recent studies have emphasized the importance of epigenetic dynamics in transcriptional regulatory networks and identified a large number of activated genomic sequences, including promoters and enhancers [11, 12]. One of the most exciting progresses is the illustration of the pivotal roles of chromatin accessibility during skeletal muscle development, which can directly reflect the effects of chromatin structural modification on gene transcription [13]. Chromatin accessibility influences the binding of transcription factors (TFs) and regulates the temporal and spatial expression patterns of target genes in muscle development, which provides an attractive opportunity to identify cis-regulatory elements [5, 14]. Previous studies have revealed that the variants residing in the cis-regulatory elements play critical roles in regulating gene expression [15]. Given its advantage in exploring key genomic regulatory elements, the Functional Annotation of Animal Genomes (FAANG) has begun to conduct large-scale chromatin accessibility analysis on muscle development in typical livestock, such as that in pigs [16], cattle [17], and chickens [18] and revealed that E2F6, OTX2, CTCF, SP1, MEF2C, EGR1, MyoD, AP-1, KLF, TEAD, and MEF2 play important roles in regulating myogenic genes expression [13, 19, 20]. Although significant progress has been made in these livestock, the landscape of chromatin accessibility of skeletal muscles remains largely unknown in rabbits.

In this study, we performed assays for transposase-accessible chromatin with high-throughput sequencing (ATAC-seq) of skeletal muscles at three key time points of growing meat rabbits. Then, the selection signals of the meat rabbit populations were analyzed to identify key genomic regions and TFs regulating the postnatal development of rabbit skeletal muscles. Our study provides a category of potential cis-regulatory elements for understanding the development of skeletal muscles and is expected to shed light on the molecular mechanisms regulating postnatal skeletal muscle development in rabbits.

Methods

Ethics approval

All experiments were performed in accordance with relevant guidelines and adhered to the ARRIVE guidelines (https://arriveguidelines.org/) for the reporting of animal experiments. This study was carried out in accordance with the principles of the Basel Declaration and recommendations of the Guide for the Care and Use of Laboratory Animals (http://grants1.nih.gov/grants/olaw/references/phspol.htm). All surgical procedures involving pigs were performed according to the approved protocols of the Biological Studies Animal Care and Use Committee, Sichuan Province, China. The protocol was approved by the ethics committee of Mianyang Normal University under permit No. SKY101368.

Experimental rabbits and histological analysis

In this study, the New Zealand white (NZW) rabbits were raised in the animal experiment center of Mianyang Normal University, Mianyang, China. These rabbits were breastfed for 35 days after birth and then fed a standard diet and water ad libitum, as described in our previous studies [21]. The skeletal muscle tissues were collected from the left hind thigh of six male rabbits at three representative time points of 1 day (D1), 35 days (D35), and 75 days (D75), representing the newborn, weaning, and adult rabbits, respectively. These skeletal muscle tissues were immediately isolated from rabbits under sterile conditions and then were snap-frozen in liquid nitrogen and stored at − 80℃ until cell isolation. The tissue samples used for histological assay were fixed using 4% paraformaldehyde at 4℃ overnight. Then, the fixed tissues were embedded in paraffin and sliced using a microtome (Leica, Bensheim, Germany). Haematoxylin and eosin (H&E) staining was performed after deparaffination according to standard protocols described in our previous study [22]. All slices of the H&E assay were observed and imaged using an Olympus BX-50 F light microscope (Olympus Optical, Tokyo, Japan). The diameters of muscle fibers were quantified by detecting the area of cross-sectional area of H&E slices using ImageJ [23] based on the scaled length in the images of H&E slices.

ATAC-seq library Preparation and sequencing

Two biological replicates of skeletal muscles were used to construct ATAC-seq libraries using previously published methods [24, 25]. Briefly, the skeletal muscle tissues were digested into cell suspension, which was re-suspended in the nuclear isolation buffer and washed repeatedly using nuclear wash buffer following a standard nuclear isolation protocol. Then, the DNA of approximately 60,000 nuclei was transposed using Tn5 transposase for 40 min at 37℃ using the Hyperactive ATAC-Seq Library Prep Kit for Illumina (Vazyme, Nanjing, China) according to the manufacturer’s instructions. The transposed DNA was purified and added sequencing adapters. After PCR amplification, the libraries were purified using a Qiagen MinElute PCR Purification Kit (Qiagen Inc., Valencia, CA, United States) following the manufacturer’s instructions. Finally, the qualified libraries were sequenced on an Illumina NovaSeq 6,000 platform, and paired-end raw reads were generated.

Data processing of ATAC-seq

The sequencing adapters of raw reads and low-quality reads (N base number > 5 and mean quality < 20) were removed using Fastp software [26]. The clean reads were mapped to the rabbit reference genome (OryCun 2.0, retrieved from https://asia.ensembl.org/index.html) using Burrows-Wheeler Aligner (BWA, v0.7.17-r1188) with “mem” mode [27]. The reads with mapping quality < 20 or reads mapped to mitochondrial DNA (due to mitochondrial DNA being exposed without chromatin) were removed using the Samtools [27]. The PCR duplication caused by PCR was removed using the GATK4 “MarkDuplicates” subprogram [28]. The accessible chromatin regions (ACRs, referred to as peaks) were called using Masc2 software (v2.2.6) [29] with parameters of “--nonmodel --shift − 100 --text size 200,” and the q value cutoff for peak calling was 0.05. The genome-wide ATAC-seq peaks were annotated using CHIPSeeker (v1.36.0) [30]. The ATAC-seq peaks in each sample were merged into a consensus region set using Diffbind2 (v3.10.0) [31], and the “DBA_SCORE_READS” in DiffBind2 was subjected to calculate normalized reads using the CPM method, which was calculated by dividing raw read counts by the count number of reads that mapped to all peaks. The differential peaks were identified using DESeq2, with the thresholds|logFC| > 1 and FDR < 0.05. The gene expression profiles of skeletal muscles during rabbit growth were retrieved from our previous study, in which RNA-seq was performed on the skeletal muscles at birth (D1), weaning (D35), and adult (D70) [9]. Because D70 and D75 were very close, they were considered to have the same genomic and transcription state in this study. The reads mapped to the ACRs and genes were normalized and converted into Bigwig files [32], and then they were visualized using IGV (v2.16.2) [33].

Footprinting analysis of TFs in ATAC-seq peaks

The TF footprinting analysis of the ATAC-seq peaks was performed using Tobias software (v0.16.1) [34], in which read alignment experienced a correction of Tn5 cutting preference, estimation of footprinting score, and detection of TF binding. Briefly, the vertebrate TF motifs were downloaded from the JASPAR database (release 2023). The Tn5 transposase sequence preference of cutting sites was estimated and corrected using the subprogram ‘ATACorrect’ of Tobias. The deletion of ATAC-seq signals given rise from protein binding and the neighboring signals around binding sites were calculated using the parameter of subprogram ‘FootprintScores’ of Tobias. Subsequently, the subprogram ‘BINDetect’ of Tobias was applied to detect the difference of TF binding based on the deletion of ATAC-seq signals of corresponding sites. All TFs with -log10 (p-value) above the 95% quantile or differential binding scores smaller/larger than the 5 and 95% quantiles (top 5% in each direction) were considered differential binding TFs.

Conservation analysis of ACRs

The conservation analysis of ACRs of rabbit skeletal muscles was conducted using UCSC LiftOver tools (https://genome.ucsc.edu/cgi-bin/hgLiftOver). We referred to the classification criteria from the study of Wang and colleagues [35], and an ACR with a value of the minimum ratio of bases that must remap (minMatch) of 0.1, 0.5, and 0.99 was considered highly, moderately, lowly conserved region between rabbits (OryCun2.0) and humans (hg19). While an ACR with a lower of 0.1 minMatch was considered a non-conserved region. The external database of the UCSC Genome Browser [36] for phastCons 100-vertebrate scores [37] was used to interpret the evolutionary and biochemical relevance of ACRs in different sequence conservation.

Selection signature of meat rabbits among different populations

To detect selection signals of meat rabbits, rabbit resequencing datasets were retrieved from our and other studies [38,39,40], which contain two specialized meat rabbit strains, Tianfu black rabbits (TB, n = 34) and New Zealand white rabbits (NZW, n = 10), and two indigenous breeds in China, Fujian black rabbits (FJB, n = 8) and Sichuan white rabbits (SCW, n = 40). The accession numbers of these populations are showed in Table S1. These datasets were subjected to the de novo SNP calling process using BWA [27] and GATK4 [41]. Briefly, raw sequences in Fastq format were first subjected to quality control by removing these low-quality reads (N base number > 5 and mean quality < 20) using Fastp [26], after which we obtained the clean reads. All clean reads were mapped to the rabbit reference genome using the BWA mem algorithm [27]. The PCR duplications of mapped reads were marked and removed using Sambamba [42]. Subsequently, GATK4 software was applied to SNP calling and individual genotyping according to recommendations of GATK Best Practices [41]. We further performed the hard filtering with an expression of “QD < 2.0|| FS > 60.0|| MQ < 40.0|| MQRankSum < -12.5|| ReadPosRankSum < -8.0”, and SNPs with an absence rate of all samples greater than 50% were removed. SNP genotypes were phased and imputed using the Beagle software [31]. Breeds were divided into two groups based on whether they were specialized meat rabbit strains. To avoid the influence of hair color genes, we conducted two comparative analyses, including Tianfu black rabbits (TFB) vs. Fujian black rabbits (FJB) and New Zealand white rabbits (NZW) vs. Sichuan white rabbits (SCW), with FJB and SCW as the reference populations and TFB and NZW rabbits as the target populations. The SNPs of TFB and NZW rabbits that were under selection were identified by using the population differentiation index (FST) in Vcftools [43]. The SNP with the top 5% FST values were selected as candidate selection sites.

Functional annotation and pathway analysis

Gene Ontology in biological process (GO-BP) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) [44]. The enriched GO-BP terms or KEGG pathways with a p-value < 0.05 were considered significant.

Statistical analysis

Statistical analyses, including T-test and One-way ANOVA, were conducted on R software. The p-value < 0.05 was considered significant. The “*” represents “p-value < 0.05” in a statistical test, and the “**” represents “p-value < 0.01” in a statistical test.

Results

Histological analysis and ATAC-seq of rabbit skeletal muscle

The analysis framework of this study is shown in Fig. 1A, which contains open chromatin detection, selection analysis of population, and transcriptome analysis. The skeletal muscle tissues were isolated from New Zealand white (NZW) rabbits at growth stages of 1 day (D1), 35 days (D35), and 75 days (D75) (Fig. 1A). H&E staining assay showed that the diameters of muscle fibers significantly increased from D1 to D35 and D35 to D75 (cell area D1 = 164.44 ± 13.67 µm2, D35 = 1145.24 ± 55.00 µm2, and D75 = 3072.96 ± 142.35 µm2, Kruskal Wallis Test p < 0.01, Fig. 1B). On the other hand, the shapes of the cross-section of muscle fibers gradually change from ellipses to relatively regular polygons (Fig. 1B). To determine the global epigenetics landscape during skeletal muscle development in rabbits, we carried out ATAC-seq to interrogate the chromatin accessibility of the samples. We obtained an average of 135.49 million quantified paired-end reads per sample, of which the Q20 ratio ranged from 96.86 to 97.58% and the Q30 ratio ranged from 86.24 to 93.33% (Table S3). When performing BWA-mem alignments to reference genome, the paired-mapping ratio ranged from 89.22 to 95.82% (Table S2). Library analysis showed that we generated typical ATAC-seq data flowing two characteristics. Firstly, the insert fragments of libraries consisted of major length enrichments in short (DNA sequences from intervening regions between two consecutive nucleosomes) and secondary length enrichment in longer (DNA sequences from spanned one or more nucleosomes) (Fig. 1C). Secondly, the ATAC-seq fragments were mainly enriched in the regions around the transcriptional start sites of genes (Fig. 1D). In total, 52,731, 53,127, and 45,791 accessible chromatin regions (ACRs) were identified in the D1, D35, and D75 samples after performing peak calling, respectively (Fig. 1E). The values of the fraction of reads in peaks (FRiP) ranged from 0.06 to 0.10 (Table S2). Analysis of these ACRs found that they were widely distributed throughout the genome, and their distribution was found to be positively related to gene density in the genome in general (Fig. 1E). The DiffBind2-based merging analysis re-centred the ACRs of different samples and established 126,959 non-redundant ACRs across all samples, of which major ACRs (75.50%) were annotated to the intergenic and intronic regions, and 19.45%, 4.24%, and 0.82% were annotated to the promoter, exon, and untranslated regions (UTRs), respectively (Fig. 1E). These data indicated that we obtained genome-wide high-quality chromatin accessibility data in the skeletal muscles of rabbits.

Fig. 1
figure 1

Histology and ATAC-seq analysis of skeletal muscles at three growth stages of rabbits. (A) Experimental design in this study. (B) H&E analysis of muscle tissues of rabbits at D1, D35, and D75. The figures show the cross- (upper panel) and longitudinal section (bottom panel) of muscle fibers. The scale bars represent 50 μm length size. The boxplot shows the cross-sectional area of muscle fibers, which was quantified using ImageJ software based on the scale bars. A total of 50 cells were randomly accounted in the analysis for each sample. (C) Inserted fragments of a representative ATAC-seq library. The fragments were expanded by the paired-end mode of mapped reads. (D) Enrichment of ATAC-seq signals around gene bodies. The heatmap shows the enrichment results of one chromosome (chromosome 1) of the rabbit genome. (E) Genome-wide chromatin accessibility of samples in D1, D35, D75. The bar plot shows the genomic annotation of accessible chromatin regions (ACRs)

Full size image

Conservation analysis of ACRs

To identify conserved cis-regulatory elements, we further mapped ACRs from the rabbit to the human genome. Our data showed that 17,734, 45,842, 3211, and 64,596 regions were highly, moderately, lowly, and non-conserved (HC, MC, LC, and NC) between rabbits and humans in sequence (Fig. 2A). Characterization of regions with different types of conservation found that ACRs with high conservation tended to be distributed near TSSs (Fig. 2B), of which only approximately 6% ACRs in LC while more than 22% ACRs in HC were annotated into TSS ± 1Kb regions. On the other hand, the HC ACRs exhibit high chromosomal synteny in the genomes of rabbits and humans, such as ACRs in Chr13-Chr1, Chr1-Chr11, and Chr19-Chr17 (rabbit-human) (Fig. 2C). To assess their evolutionary and biochemical relevance, we overlapped HC, MC, and LC ACRs with the 100-vertebrate phastCons scores. Reassuringly, HC regions were more evolutionarily conserved than MC regions (median increase: 0.11; Wilcoxon test p < 2.2 × 10 − 16), and MC regions were also more evolutionarily conserved than LC regions (median increase: 0.11; Wilcoxon test p < 2.2 × 10 − 16) (Fig. 2D). GO in the biological process (GO-BP) analysis of the nearby genes of ACRs in the HC found these conserved DNA sequences were significantly enriched in tissue development, regulation of cell communication, and nervous system development (Fig. 2E). On the other hand, KEGG pathway analysis showed that these conserved DNA sequences were significantly enriched in the cytoskeleton in muscle cells, signaling pathways regulating the pluripotency of stem cells, and MAPK signaling pathway (Fig. 2E).

Fig. 2
figure 2

Conservation analysis of ACRs. (A) Identification of highly (HC), moderately (MC), lowly (LC), and non-conserved (NC) ACRs between humans and rabbits. (B) Genomic annotation of ACRs with different conservation between humans and rabbits. (C) Mapping HC ACRs from rabbits to human genomes. (D) Analysis of 100-vertebrate phastCons scores of HC, MC, and LC ACRs. (E) GO-BP and KEGG pathway analysis of nearby genes of HC ACRs

Full size image

Dynamic analysis of chromatin accessibility during the development of skeletal muscles

The principal component analysis (PCA) based on the normalized read counts of the ACRs (CPM value) showed that the first (44% variance) and the second (25% variance) principal components explained proper variances, and the samples in D1, D35, and D75 were classified into expected groups, indicating excellent reproducibility between replicates and distinct chromatin state among different growth stages of skeletal muscles (Fig. 3A). Differential analysis showed that chromatin accessibility broadly increased from D1 to D35, with a total of 2550 and 115 ACRs significantly increased and decreased chromatin accessibility in D35 vs. D1, respectively. On the other hand, chromatin accessibility broadly decreased from D35 to D75, with a total of 261 and 2949 ACRs significantly increased and decreased chromatin accessibility in D75 vs. D35, respectively. A Few regions significantly changed chromatin accessibility in the comparison of D75 vs. D1, with 18 and 21 regions increasing and decreasing chromatin accessibility (Fig. 3B and Table S3, |log2(FC)| > 1 and p-adj < 0.05). To understand the function of differential ACRs, we further subjected them to k-means clustering analysis and found that 986, 790, 319, and 2893 ACRs were classified into the four groups in different dynamic patterns (group A, B, C, and D), respectively (Fig. 3C and Table S4). The functional analysis found that ACRs in different groups were responsible for different biological processes and signalling pathways. The chromatin of regions in group A was highly accessible at D1 and D35 and almost closed at D75. GO analysis of nearby genes of regions in group A found that they were mainly involved in metabolism and gene transcription-related functions. The top 5 significantly enriched terms in the GO-BP category were “positive regulation of metabolic process”, “regulation of gene expression”, “transcription, DNA-templated”, “positive regulation of macromolecule metabolic process”, and “positive regulation of cellular metabolic process” (Fig. 3D). KEGG pathway analysis of nearby genes of regions in group A found that they were significantly enriched in the AMPK signalling pathway (Fig. 3E). The chromatin of regions in group B was almost closed at D1 and highly accessible at D35 and D75. GO analysis of nearby genes of regions in group B found that they were mainly involved in skeletal muscle development-related functions. The skeletal muscle development-related terms in the top 10 significantly enriched GO-BP terms were “animal organ morphogenesis”, “muscle cell differentiation”, and “muscle structure development” (Fig. 3D). The top 5 significantly enriched KEGG pathways were “Gastric acid secretion”, “Oxytocin signalling pathway”, “cGMP-PKG signalling pathway”, “Calcium signalling pathway”, and “cAMP signalling pathway” (Fig. 3E). The chromatin of regions in group C decreased accessibility from D1 to D35 and increased accessibility from D35 to D75. GO analysis of nearby genes of regions in group C found that they were mainly involved in neuron development-related functions. The top 5 significantly enriched GO-BP terms were “neuron differentiation”, “neuron development”, “generation of neurons”, “neurogenesis”, and “nervous system development” (Fig. 3D). The chromatin of regions in group D increased accessibility from D1 to D35 and decreased chromatin accessibility from D35 to D75. GO analysis of nearby genes of regions in group D found that they were mainly involved in signal transduction-related functions. The top 5 significantly enriched GO-BP terms were “regulation of cellular component organization”, “regulation of cell communication”, “regulation of signaling”, “response to endogenous stimulus”, and “regulation of signal transduction” (Fig. 3D). KEGG analysis found that the top 5 significantly enriched signaling pathways were “Hippo signaling pathway”, “Pathways in cancer”, “Wnt signaling pathway”, “Calcium signaling pathway”, and “Cushing syndrome” (Fig. 3E).

Fig. 3
figure 3

Differential analysis of chromatin accessibility of skeletal muscles in rabbits. (A) PCA analysis of samples based on the normalized intensity of ATAC-seq signals in the ACRs. (B) Differential analysis of the intensity of ACRs. The red and blue points represent the significantly increased and decreased ACRs, respectively. (C) Clustering analysis of differential ACRs. The line plot in the zoomed annotation of the heatmap shows the change trend in the chromatin accessibility of corresponding groups. (D) GO-BP analysis of nearby genes of ACRs in different groups. The gene count was divided by ten. (E) KEGG pathway analysis of nearby genes of ACRs in different groups. The gene count was divided by ten

Full size image

Integrated analysis of ATAC-seq and transcriptome data revealed key cis-regulation patterns of genes

To identify key cis-regulatory elements for the corresponding gene, we integrated the ATAC-seq in this study and the RNA-seq data from our previous study [9]. All the ACRs were classified into promoters (ACRs in TSS ± 3 kb) and enhancers (ACRs not in promoters), and then identified different cis-regulation types of genes. The first cis-regulation type of genes is promoter-domination (PD), in which the chromatin accessibility of enhancer ACR(s) was not increased, but that of promoter ACR(s) was synchronously increased with the expression of a nearby gene. For example, the upregulated TSGA10 contains one increased ACR in the promoter region but contains 0 ACR in the enhancer regions (Fig. 4A). The second cis-regulation type of genes is the enhancer-domination (ED), in which the chromatin accessibility of enhancer ACR(s) was increased, but that of promoter ACR(s) was not synchronously increased with the expression of nearby gene(s). For example, the upregulated DTL contains ACRs in promoter and enhancer regions, but only one ACR in the enhancer region increased chromatin accessibility (Fig. 4A). The third cis-regulation type of genes is promoter and enhancer co-regulation (PEC), in which the chromatin accessibility of both promoter and enhancer ACR(s) was synchronously increased with the expression of nearby genes. For example, the chromatin accessibilities of both ACR in the promoter region and ACR in the enhancer regions of ENSOCUT00000058429 were synchronously increased with its expression (Fig. 4A). Based on the classification criteria, a total of 390 genes in ED, 21 genes in PD, and 1 gene in PEC were identified from D1 to D35 (Fig. 4B). On the other hand, a total of 292 genes in ED, 96 genes in PD, and 4 genes in PEC were identified from D35 to D75 (Fig. 4C). Because genes in ED and PD were the major cis-regulatory type, we further analyzed the functional difference between genes in ED and PD. Our data showed that the genes in ED mainly governed neuron development. The top 5 significantly enriched GO-BP terms were “head development”, “animal organ morphogenesis”, “central nervous system development”, “positive regulation of neurogenesis”, and “regulation of cellular component organization”. While the genes in PD mainly governed cell cycle-regulated functions. The top 5 significantly enriched GO-BP terms were “mitotic cell cycle process”, “regulation of cell cycle”, “cell cycle process”, “mitotic cell cycle”, and “transcription, DNA-templated” (Fig. 4C).

Fig. 4
figure 4

Integrated analysis of ATAC-seq and transcriptome data revealed key cis-regulatory elements of rabbit skeletal muscle. (A) Tracks of RNA-seq and ATAC-seq data. Tracks of representative genes in different types of cis-regulation, including promoter-domination, enhancer-domination, and promoter and enhancer co-regulation. Each track represents a single RNA-seq (blue) or ATAC-seq (red) sample. Increased ACRs were marked with grey shadows. (B) Summary of genes in different cis-regulatory types. (C) GO analysis of genes in ED and PD. EDG: enhancer-domination genes; PDG: promoter-domination genes

Full size image

Footprinting and selection signal analyses identified key TFs and mutations in ACRs

Chromatin accessibility is a key epigenetic factor that physically mediates transcription factors and DNA. We further performed footprinting analysis using the ATAC-seq after correcting Tn5 cutting preference, which provides evidence of direct occupancy of TF candidates on genomic DNA. For instance, the IRF7 that has not been reported in the skeletal muscle development in mammals was found to have deeper footprints and higher DNA accessibility, flanking its motif in D35 samples than in D1 and D75 samples (Fig. 5A). An average of 151 binding sites of TFs were found in the differential ACRs. A total of 17 TFs have > 1000 binding sites in the differential ACRs, which were mainly the members of SP and KLF families (Table S5 and Figure S1A). Statistics analysis of the subjected TFs found that 24 and 30 TFs significantly increased and decreased footprints from D1 to D35. On the other hand, 29 and 22 TFs significantly increased and decreased footprints from D35 to D75 (Fig. 5B). To identify key TFs and their target genes, we further performed function-based filtering of TFs. The top 5 TFs that have the most significantly increased footprints from D1 to D35 were ZNF354A, ZNF558, PHOX2B, HMGA1, and FOXD3, and the top 5 TFs that have the most significantly increased footprints from D35 to D75 were PATZ1, SP2, KLF10, ZNF148, and ZBEB4.

FST was used to select the top 5% of SNPs, with Fujian black rabbits (FJB) as the control group and Tianfu black rabbits (TFB) as the selection group, 10,793 SNPs with a selection signature were obtained (Fig. 5C). With Sichuan white rabbits (SCW) as the control group and New Zealand white rabbits (NZW) as the selection group, 110,747 genomic windows with a selection signature were obtained. An intersection analysis showed that a total of 2947 selection SNPs were found between the two FST comparisons (Figure S1B). We then integrated the results of the selection signals and TF footprinting analysis and found that 6 SNPs (Chr2:20258634, Chr2:150969157, Chr2:150969186, Chr3:54676164, Chr13:121449670, Chr20:30932071) with selection signals were located at the TF binding sites (Fig. 5D), among which one SNP (Chr13:121449670) was found as base transversion (A-> G) and located at a position with high positional weight of a key muscle TF PRDM1 [45,46,47]. Significantly, the “A” base of this position was found to have a high allele frequency (TFB, 82.4%; NZW, 90.0%) in both two specialized meat rabbit strains compared to the two indigenous breeds (FJB, 31.25%; SCW, 46.25%), while “G” base was found to have a low allele frequency (TFB, 17.6%; NZW, 10.0%) in both two specialized meat rabbit strains compared to the two indigenous breeds (FJB, 68.75%; SCW, 53.75%) (Fig. 5E). Combining the results of differential analysis of ATAC-seq peaks found that the region where this SNP is located maintains high chromatin accessibility at D35, while chromatin accessibility is lower at D1 and D75 and might affect nearby genes of TRL6 and ZSWIM5 (Fig. 5E).

Fig. 5
figure 5

TF footprinting analysis of ACRs. (A) Visualization of ATAC-seq footprint for motifs of an example TF (IRF7) in three developmental stages of rabbit muscles. (B) Differential binding analysis of TFs. The red and blue points showed TFs with increased and decreased footprints. The normal and small points represent mRNA-expressed and unexpressed TF genes. ATAC-seq signals across all motif binding sites in differential ACRs were aligned on the motifs and averaged. (C) Detection of selection signals of meat rabbits using the FST method. The upper panel indicates the selected SNPs by using Fujian black rabbits and Tianfu black rabbits as reference and target populations, respectively. The bottom panel indicates the selected SNPs by using Sichuan white rabbits and New Zealand white rabbits as reference and target populations, respectively. The scatters above the dashed line indicate the selected SNPs. (D) Summary of the SNPs located in the binding sites of TFs. (E) A key SNP located at the binding site of PRDM1. The orange tracks indicate the reads coverage of ATAC-seq. The variant A/G was marked using a red box. The allele frequencies of Chr13:121449670 in different populations are indicated using a bar plot

Full size image

Discussion

In this study, we performed histological analysis for the skeletal muscle of rabbits and found the diameters of muscle fibers continuously increased after rabbit birth, which was in line with other mammal livestock, such as cattle, pigs, and goats [5, 48, 49]. Newborn (D1), weaning (D35), and adulthood (D75) are the three key time points for skeletal muscle development in meat rabbits. ATAC-seq has rapidly become the preferred approach for investigating activated genomic cis-regulatory elements. Many previous studies have utilized the efficiency of ATAC-seq data to explore the regulation of muscle development in common livestock and identified key promoters and enhancers [5, 50, 51]. In this study, we carried out ATAC-seq technology to detect accessible chromatin regions in rabbit skeletal muscle tissues, which is the first comprehensive assessment of activated genomic regions in rabbit muscles. The distribution of length sizes of inserted fragments and enrichment characteristics of open chromatin signals were in line with previous ATAC-seq [52, 53], indicating the reliability of our ATAC-seq data. In general, genomic cis-regulatory elements were classified into two types, including promoters [54] and enhancers [55]. In our results, genomic annotation of identified ACRs showed that major activated genomic regions were located in the intergenic and intron sequences, suggesting that a large number of enhancers regulated skeletal muscle development in rabbits. Previous studies have revealed that dynamics of chromatin accessibility play important roles in different types of tissues of domestic animals, such as adipose tissue [56] and skeletal muscle tissues [20]. In this study, we found thousands of ACRs increased chromatin accessibility from D1 to D35 and decreased chromatin accessibility from D35 to D75, which indicated that chromatin accessibility is a key epigenetic factor regulating skeletal muscle development in rabbits. Furthermore, muscle development-related genes were found to increase chromatin accessibility in the promoter regions, suggesting that chromatin accessibility plays a crucial role in gene expression during skeletal muscle development in rabbits. Interestingly, our data showed that there were a large number of ACRs that increased chromatin accessibility from D1 to D35, but very few differential ACRs were found in the comparison of D75 vs. D1, which might suggest that the period from newborn to weaning is important for obtaining sufficient epigenetic regulation.

The clustering analysis found that the ACRs in different dynamic patterns play different roles during skeletal muscle development in rabbits. Previous studies have revealed that the neuron system develops coordinatedly with skeletal muscle in mammals [57]. Our data showed that ACRs that decreased their accessibility from D1 to D35 and increased accessibility at D75 were involved in nervous system development, which indicates chromatin accessibility is crucial for maintaining the function of the neuromuscular system. Previous studies have revealed that skeletal muscle is a complex heterogeneous tissue composed of various cell types [58], and single-cell ATAC-seq is used to analyze skeletal muscle development in humans and mice [59, 60]. Currently, single-cell ATAC-seq is gradually being utilized in model animals and major livestock, such as mice [61], pigs [13], and cattle [62]. However, single-cell ATAC-seq technology on rabbit muscles has not yet been applied. Although epigenetics in single-cell resolution has not been reported in rabbit muscle, considering that muscle fibers are the main cell type in muscle tissue, we obtained a group of ACRs that were involved in muscle development in this study by clustering analysis (group B). Therefore, the ACRs in group B might be key cis-regulatory elements for developing muscle fibres in rabbits. Integrating ATAC-seq and RNA-seq data to identify key cis-regulatory elements is an effective and scientific approach in many previous studies [13, 63, 64]. Uncovering the cis-regulatory elements that govern when and how much each gene is transcribed in a given genome and cellular state remains a central goal of biology [65]. ATAC-seq can reveal the open state of chromatin, indirectly reflecting the activity of regulatory elements; RNA-seq provides information on gene expression levels. Combining these two types of data can more accurately identify cis-regulatory elements closely related to gene expression. In this study, we provide a more detailed analysis of the ACRs surrounding genes. Our data showed that ED was the major type of cis-regulation, which might indicate that distal regulation is an important mode of muscle development, and further analysis of the three-dimensional spatial structure of chromatin may be necessary for precise analysis of cis-regulation in the future. On the other hand, we found genes in proximal regulation (PD) were involved in cell cycle-related function, and genes in distal regulation (ED) were involved in neuron-related functions, which might suggest that promoters mainly govern basic or non-cell-selective biological functions and enhancers mainly govern specific and more complex biological functions.

Epigenetic mechanisms governing chromatin organization and TF binding are critical components of transcriptional regulation [34]. The open chromatin physically provides an opportunity for the binding of transcription factors. In this context, DNA sequences are passive, while transcription factors are active, and the engagement of TFs and binding sites promotes transcriptional effects. The members of the SP family could improve muscle endurance [66] and were enriched in many binding sites during the postanal development of skeletal muscle in pigs in the previous study [20]. On the other hand, the members of the KLF family were key TFs during muscle development [67]. For example, KLF7 retaining muscle stem satellite cell quiescence is important for the maintenance of stem cell population and tissue regeneration [68]. KLF15 was reported to play an important role in skeletal muscle energy metabolism [69] and positively regulate skeletal muscle differentiation [70]. Our data showed that the members of the SP and KLF family have the most binding sites in the differential ACRs and were highly active in the adult stage in rabbit muscles, which might indicate the activities of the SP and KLF family play key roles in terminal development of skeletal muscle in rabbits. On the other hand, we screen out differential binding TFs by comparing muscles from different stages, which indicates a large number of TFs regulate the development of rabbit muscles through interacting with DNA sequences of chromatin accessibility regions. To further understand the mechanisms of chromatin accessibility governing muscle development, we performed an integrative analysis of population data and chromatin accessibility data. Previous studies have revealed that PRDM1 plays key roles in muscle development, including differentiation [45], fibre type commitment [47], and organogenesis [71]. In this study, the population analysis showed that a key selected SNP maker was located in the PRDM1 binding site and there is a significant tendency in allele frequency between specialized meat rabbit strains and indigenous breeds. On the other hand, the PRDM1 binding site maintains high accessibility during the growth stage in rabbits, suggesting that two mechanisms may affect the regulation of rabbit muscle development by PRDM1, including the single base transversion (A-> G) and the dynamic changes in chromatin accessibility of DNA sequence. A previous study has revealed that ZSWIM8 is a myogenic protein and regulates C2C12 differentiation [72]. Our study showed that ZSWIM5, which was from the protein family of ZSWIM8 is the nearby gene of the key selected SNP maker. Therefore, the key SNP (Chr13:12144967 A-> G), PRMD1 binding site, and the ZSWIM5 thus warrant further investigation.

Conclusions

In summary, we identify genome-wide accessible chromatin regions and their potential regulatory mechanisms in rabbit skeletal muscles at the tissue level. Our study provided a category of potential cis-regulatory elements for understanding the development of skeletal muscles at key growth stages in rabbits and thus might facilitate potential insights into meat rabbit breeding.

Data availability

The raw sequencing datasets generated in this study have been deposited in the Sequence Read Archive at NCBI under accession code PRJNA1171873 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA1171873).

References

  1. Li S, Zeng W, Li R, Hoffman LC, He Z, Sun Q, Li H. Rabbit meat production and processing in China. Meat Sci. 2018;145:320–8. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.meatsci.2018.06.037

    Article  PubMed  Google Scholar 

  2. Ebeid TA, Tůmová E, Al-Homidan IH, Ketta M, Chodová D. The potential role of feed restriction on productivity, carcass composition, meat quality, and muscle fibre properties of growing rabbits: A review. Meat Sci. 2022;191:108845. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.meatsci.2022.108845

    Article  CAS  PubMed  Google Scholar 

  3. Li CY, Li X, Liu Z, Ni W, Zhang X, Hazi W, Ma Q, Zhang Y, Cao Y, Qi J, et al. Identification and characterization of long non-coding RNA in prenatal and postnatal skeletal muscle of sheep. Genomics. 2019;111(2):133–41. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.ygeno.2018.01.009

    Article  CAS  PubMed  Google Scholar 

  4. Cisternas P, Henriquez JP, Brandan E, Inestrosa NC. Wnt signaling in skeletal muscle dynamics: myogenesis, neuromuscular synapse and fibrosis. Mol Neurobiol. 2014;49(1):574–89. https://doiorg.publicaciones.saludcastillayleon.es/10.1007/s12035-013-8540-5

    Article  CAS  PubMed  Google Scholar 

  5. Yue J, Hou X, Liu X, Wang L, Gao H, Zhao F, Shi L, Shi L, Yan H, Deng T, et al. The landscape of chromatin accessibility in skeletal muscle during embryonic development in pigs. J Anim Sci Biotechnol. 2021;12(1):56. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s40104-021-00577-z

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Brinegar AE, Xia Z, Loehr JA, Li W, Rodney GG, Cooper TA. Extensive alternative splicing transitions during postnatal skeletal muscle development are required for calcium handling functions. Elife. 2017;6:e27192. https://doiorg.publicaciones.saludcastillayleon.es/10.7554/eLife.27192

    Article  PubMed  PubMed Central  Google Scholar 

  7. Li Y, Wang J, Elzo MA, Gan M, Tang T, Shao J, Lai T, Ma Y, Jia X, Lai S. Multi-Omics analysis of key microRNA-mRNA metabolic regulatory networks in skeletal muscle of obese rabbits. Int J Mol Sci. 2021;22(8):4204. https://doiorg.publicaciones.saludcastillayleon.es/10.3390/ijms22084204

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Kuang L, Lei M, Li C, Zhang X, Ren Y, Zheng J, Guo Z, Zhang C, Yang C, Mei X, et al. Identification of long Non-Coding RNAs related to skeletal muscle development in two rabbit breeds with different growth rate. Int J Mol Sci. 2018;19(7):2046. https://doiorg.publicaciones.saludcastillayleon.es/10.3390/ijms19072046

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Du K, Zhao X, Li Y, Wu Z, Sun W, Wang J, Jia X, Chen S, Lai S. Genome-Wide identification and characterization of circular RNAs during skeletal muscle development in meat rabbits. Anim (Basel). 2022;12(17):2208. https://doiorg.publicaciones.saludcastillayleon.es/10.3390/ani12172208

    Article  Google Scholar 

  10. Yang X, Deng F, Wu Z, Chen SY, Shi Y, Jia X, Hu S, Wang J, Cao W, Lai SJ. A Genome-Wide association study identifying genetic variants associated with growth, carcass and meat quality traits in rabbits. Anim (Basel). 2020;10(6):1068. https://doiorg.publicaciones.saludcastillayleon.es/10.3390/ani10061068

    Article  Google Scholar 

  11. Yan F, Powell DR, Curtis DJ, Wong NC. From reads to insight: a Hitchhiker’s guide to ATAC-seq data analysis. Genome Biol. 2020;21(1):22. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s13059-020-1929-3

    Article  PubMed  PubMed Central  Google Scholar 

  12. Pan Z, Yao Y, Yin H, Cai Z, Wang Y, Bai L, Kern C, Halstead M, Chanthavixay G, Trakooljul N, et al. Pig genome functional annotation enhances the biological interpretation of complex traits and human disease. Nat Commun. 2021;12(1):5848. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/s41467-021-26153-7

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Cai S, Hu B, Wang X, Liu T, Lin Z, Tong X, Xu R, Chen M, Duo T, Zhu Q, et al. Integrative single-cell RNA-seq and ATAC-seq analysis of myogenic differentiation in pig. BMC Biol. 2023;21(1):19. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12915-023-01519-z

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Dong A, Preusch CB, So WK, Lin K, Luan S, Yi R, Wong JW, Wu Z, Cheung TH. A long noncoding RNA, LncMyoD, modulates chromatin accessibility to regulate muscle stem cell myogenic lineage progression. Proc Natl Acad Sci U S A. 2020;117(51):32464–75. https://doiorg.publicaciones.saludcastillayleon.es/10.1073/pnas.2005868117

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Miao Y, Zhao Y, Wan S, Mei Q, Wang H, Fu C, Li X, Zhao S, Xu X, Xiang T. Integrated analysis of genome-wide association studies and 3D epigenomic characteristics reveal the BMP2 gene regulating loin muscle depth in Yorkshire pigs. PLoS Genet. 2023;19(6):e1010820. https://doiorg.publicaciones.saludcastillayleon.es/10.1371/journal.pgen.1010820

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Andersson L, Archibald AL, Bottema CD, Brauning R, Burgess SC, Burt DW, Casas E, Cheng HH, Clarke L, Couldrey C, et al. Coordinated international action to accelerate genome-to-phenome with FAANG, the functional annotation of animal genomes project. Genome Biol. 2015;16(1):57. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s13059-015-0622-4

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Alexandre PA, Naval-Sánchez M, Menzies M, Nguyen LT, Porto-Neto LR, Fortes MRS, Reverter A. Chromatin accessibility and regulatory vocabulary across indicine cattle tissues. Genome Biol. 2021;22(1):273. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s13059-021-02489-7

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Zhu XN, Wang YZ, Li C, Wu HY, Zhang R, Hu XX. Chicken chromatin accessibility atlas accelerates epigenetic annotation of birds and gene fine-mapping associated with growth traits. Zool Res. 2023;44(1):53–62. https://doiorg.publicaciones.saludcastillayleon.es/10.24272/j.issn.2095-8137.2022.228

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Lyu P, Settlage RE, Jiang H. Genome-wide identification of enhancers and transcription factors regulating the myogenic differentiation of bovine satellite cells. BMC Genomics. 2021;22(1):901. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12864-021-08224-7

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Feng L, Si J, Yue J, Zhao M, Qi W, Zhu S, Mo J, Wang L, Lan G, Liang J. The landscape of accessible chromatin and developmental transcriptome maps reveal a genetic mechanism of skeletal muscle development in pigs. Int J Mol Sci. 2023;24(7):6413. https://doiorg.publicaciones.saludcastillayleon.es/10.3390/ijms24076413

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Wang G, Guo G, Tian X, Hu S, Du K, Zhang Q, Mao J, Jia X, Chen S, Wang J, et al. Screening and identification of MicroRNAs expressed in perirenal adipose tissue during rabbit growth. Lipids Health Dis. 2020;19(1):35. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12944-020-01219-5

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Wang G, Du K, Xie Z, Tang R, Jia X, Chen S, Lai S. Screening and identification of differentially expressed and adipose Growth-Related Protein-Coding genes during the deposition of perirenal adipose tissue in rabbits. Diabetes Metab Syndr Obes. 2020;13:4669–80. https://doiorg.publicaciones.saludcastillayleon.es/10.2147/dmso.S284246

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Schneider CA, Rasband WS, Eliceiri KW. NIH image to imageJ: 25 years of image analysis. Nat Methods. 2012;9(7):671–5. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/nmeth.2089

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Cao Y, Ai Y, Zhang X, Zhang J, Long X, Zhu Y, Wang L, Gu Q, Han H. Genome-wide epigenetic dynamics during postnatal skeletal muscle growth in Hu sheep. Commun Biol. 2023;6(1):1077. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/s42003-023-05439-0

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Grandi FC, Modi H, Kampman L, Corces MR. Chromatin accessibility profiling by ATAC-seq. Nat Protoc. 2022;17(6):1518–52. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/s41596-022-00692-9

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Chen S, Zhou Y, Chen Y, Gu J. Fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):i884–90. https://doiorg.publicaciones.saludcastillayleon.es/10.1093/bioinformatics/bty560

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25(14):1754–60. https://doiorg.publicaciones.saludcastillayleon.es/10.1093/bioinformatics/btp324

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Babadi M, Fu JM, Lee SK, Smirnov AN, Gauthier LD, Walker M, Benjamin DI, Zhao X, Karczewski KJ, Wong I, et al. GATK-gCNV enables the discovery of rare copy number variants from exome sequencing data. Nat Genet. 2023;55(9):1589–97. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/s41588-023-01449-0

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, Nusbaum C, Myers RM, Brown M, Li W, et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 2008;9(9):R137. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/gb-2008-9-9-r137

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Yu G, Wang LG, He QY. ChIPseeker: an R/Bioconductor package for chip peak annotation, comparison and visualization. Bioinformatics. 2015;31(14):2382–3. https://doiorg.publicaciones.saludcastillayleon.es/10.1093/bioinformatics/btv145

    Article  CAS  PubMed  Google Scholar 

  31. Ross-Innes CS, Stark R, Teschendorff AE, Holmes KA, Ali HR, Dunning MJ, Brown GD, Gojis O, Ellis IO, Green AR, et al. Differential oestrogen receptor binding is associated with clinical outcome in breast cancer. Nature. 2012;481(7381):389–93. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/nature10730

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Ramírez F, Dündar F, Diehl S, Grüning BA, Manke T. DeepTools: a flexible platform for exploring deep-sequencing data. Nucleic Acids Res. 2014;42(Web Server issue):W187–191. https://doiorg.publicaciones.saludcastillayleon.es/10.1093/nar/gku365

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Robinson JT, Thorvaldsdóttir H, Wenger AM, Zehir A, Mesirov JP. Variant review with the integrative genomics viewer. Cancer Res. 2017;77(21):e31–4. https://doiorg.publicaciones.saludcastillayleon.es/10.1158/0008-5472.Can-17-0337

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Bentsen M, Goymann P, Schultheis H, Klee K, Petrova A, Wiegandt R, Fust A, Preussner J, Kuenne C, Braun T, et al. ATAC-seq footprinting unravels kinetics of transcription factor binding during zygotic genome activation. Nat Commun. 2020;11(1):4267. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/s41467-020-18035-1

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Zhu X, Ma S, Wong WH. Genetic effects of sequence-conserved enhancer-like elements on human complex traits. Genome Biol. 2024;25(1):1. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s13059-023-03142-1

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Kuhn RM, Haussler D, Kent WJ. The UCSC genome browser and associated tools. Brief Bioinform. 2013;14(2):144–61. https://doiorg.publicaciones.saludcastillayleon.es/10.1093/bib/bbs038

    Article  CAS  PubMed  Google Scholar 

  37. Siepel A, Bejerano G, Pedersen JS, Hinrichs AS, Hou M, Rosenbloom K, Clawson H, Spieth J, Hillier LW, Richards S, et al. Evolutionarily conserved elements in vertebrate, insect, worm, and yeast genomes. Genome Res. 2005;15(8):1034–50. https://doiorg.publicaciones.saludcastillayleon.es/10.1101/gr.3715005

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Chen XW, Li L, Yin M, Wang Q, Luo WT, Ma Y, Pu ZH, Zhou JL. Cloning and molecular characterization of the ORF5 gene from a PRRSV-SN strain from Southwest China. Microb Pathog. 2017;112:295–302. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.micpath.2017.09.011

    Article  CAS  PubMed  Google Scholar 

  39. Liu C, Wang S, Dong X, Zhao J, Ye X, Gong R, Ren Z. Exploring the genomic resources and analysing the genetic diversity and population structure of Chinese Indigenous rabbit breeds by RAD-seq. BMC Genomics. 2021;22(1):573. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12864-021-07833-6

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Fatima N, Jia L, Liu B, Li L, Bai L, Wang W, Zhao S, Wang R, Liu E. A homozygous missense mutation in the fibroblast growth factor 5 gene is associated with the long-hair trait in Angora rabbits. BMC Genomics. 2023;24(1):298. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12864-023-09405-2

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. DePristo MA, Banks E, Poplin R, Garimella KV, Maguire JR, Hartl C, Philippakis AA, del Angel G, Rivas MA, Hanna M, et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat Genet. 2011;43(5):491–8. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/ng.806

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Tarasov A, Vilella AJ, Cuppen E, Nijman IJ, Prins P. Sambamba: fast processing of NGS alignment formats. Bioinformatics. 2015;31(12):2032–4. https://doiorg.publicaciones.saludcastillayleon.es/10.1093/bioinformatics/btv098

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, Handsaker RE, Lunter G, Marth GT, Sherry ST, et al. The variant call format and vcftools. Bioinformatics. 2011;27(15):2156–8. https://doiorg.publicaciones.saludcastillayleon.es/10.1093/bioinformatics/btr330

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/nprot.2008.211

    Article  CAS  PubMed  Google Scholar 

  45. Vincent SD, Mayeuf A, Niro C, Saitou M, Buckingham M. Non conservation of function for the evolutionarily conserved prdm1 protein in the control of the slow twitch myogenic program in the mouse embryo. Mol Biol Evol. 2012;29(10):3181–91. https://doiorg.publicaciones.saludcastillayleon.es/10.1093/molbev/mss125

    Article  CAS  PubMed  Google Scholar 

  46. Elworthy S, Hargrave M, Knight R, Mebus K, Ingham PW. Expression of multiple slow myosin heavy chain genes reveals a diversity of zebrafish slow twitch muscle fibres with differing requirements for Hedgehog and Prdm1 activity. Development. 2008;135(12):2115–26. https://doiorg.publicaciones.saludcastillayleon.es/10.1242/dev.015719

    Article  CAS  PubMed  Google Scholar 

  47. von Hofsten J, Elworthy S, Gilchrist MJ, Smith JC, Wardle FC, Ingham PW. Prdm1- and Sox6-mediated transcriptional repression specifies muscle fibre type in the zebrafish embryo. EMBO Rep. 2008;9(7):683–9. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/embor.2008.73

    Article  CAS  Google Scholar 

  48. Picard B, Gagaoua M. Muscle fiber properties in cattle and their relationships with meat qualities: an overview. J Agric Food Chem. 2020;68(22):6021–39. https://doiorg.publicaciones.saludcastillayleon.es/10.1021/acs.jafc.0c02086

    Article  CAS  PubMed  Google Scholar 

  49. Huang J, Jiao J, Tan ZL, He Z, Beauchemin KA, Forster R, Han XF, Tang SX, Kang J, Zhou C. Inferring the skeletal muscle developmental changes of grazing and Barn-Fed goats from gene expression data. J Agric Food Chem. 2016;64(36):6791–800. https://doiorg.publicaciones.saludcastillayleon.es/10.1021/acs.jafc.6b02708

    Article  CAS  PubMed  Google Scholar 

  50. Kaplow IM, Schäffer DE, Wirthlin ME, Lawler AJ, Brown AR, Kleyman M, Pfenning AR. Inferring mammalian tissue-specific regulatory conservation by predicting tissue-specific differences in open chromatin. BMC Genomics. 2022;23(1):291. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12864-022-08450-7

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Li Q, Wang Y, Hu X, Zhang Y, Li H, Zhang Q, Cai W, Wang Z, Zhu B, Xu L, et al. Transcriptional States and chromatin accessibility during bovine myoblasts proliferation and myogenic differentiation. Cell Prolif. 2022;55(5):e13219. https://doiorg.publicaciones.saludcastillayleon.es/10.1111/cpr.13219

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Zhao Z, Guo D, Wei Y, Li J, Jia X, Niu Y, Liu Z, Bai Y, Chen Z, Shi B, et al. Integrative ATAC-seq and RNA-seq analysis of the longissimus dorsi muscle of Gannan Yak and Jeryak. Int J Mol Sci. 2024;25(11):6029. https://doiorg.publicaciones.saludcastillayleon.es/10.3390/ijms25116029

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Xu Z, Wu J, Zhou J, Zhang Y, Qiao M, Sun H, Li Z, Li L, Chen N, Oyelami FO, et al. Integration of ATAC-seq and RNA-seq analysis identifies key genes affecting intramuscular fat content in pigs. Front Nutr. 2022;9:1016956. https://doiorg.publicaciones.saludcastillayleon.es/10.3389/fnut.2022.1016956

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Xu X, Mo Q, Cai Z, Jiang Q, Zhou D, Yi J, Promoters. Key Cis-Regulatory elements, and their potential applications in regulation of cadmium (Cd) in rice. Int J Mol Sci. 2024;25(24):13237. https://doiorg.publicaciones.saludcastillayleon.es/10.3390/ijms252413237

  55. Kreibich E, Kleinendorst R, Barzaghi G, Kaspar S, Krebs AR. Single-molecule footprinting identifies context-dependent regulation of enhancers by DNA methylation. Mol Cell. 2023;83(5):787–e802789. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.molcel.2023.01.017

    Article  CAS  PubMed  Google Scholar 

  56. Du K, Chen GH, Bai X, Chen L, Hu SQ, Li YH, Wang GZ, He JW, Lai SJ. Dynamics of transcriptome and chromatin accessibility revealed sequential regulation of potential transcription factors during the brown adipose tissue whitening in rabbits. Front Cell Dev Biol. 2022;10:981661. https://doiorg.publicaciones.saludcastillayleon.es/10.3389/fcell.2022.981661

    Article  PubMed  PubMed Central  Google Scholar 

  57. Zmojdzian M, Jagla K. The relationship between muscle stem cells and motor neurons. Cell Mol Life Sci. 2021;78(12):5043–9. https://doiorg.publicaciones.saludcastillayleon.es/10.1007/s00018-021-03838-2

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Cai C, Yue Y, Yue B. Single-cell RNA sequencing in skeletal muscle developmental biology. Biomed Pharmacother. 2023;162:114631. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.biopha.2023.114631

    Article  CAS  PubMed  Google Scholar 

  59. Sahinyan K, Blackburn DM, Simon MM, Lazure F, Kwan T, Bourque G, Soleimani VD. Application of ATAC-Seq for genome-wide analysis of the chromatin state at single myofiber resolution. Elife. 2022;11:e72792. https://doiorg.publicaciones.saludcastillayleon.es/10.7554/eLife.72792

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Orchard P, Manickam N, Ventresca C, Vadlamudi S, Varshney A, Rai V, Kaplan J, Lalancette C, Mohlke KL, Gallagher K, et al. Human and rat skeletal muscle single-nuclei multi-omic integrative analyses nominate causal cell types, regulatory elements, and SNPs for complex traits. Genome Res. 2021;31(12):2258–75. https://doiorg.publicaciones.saludcastillayleon.es/10.1101/gr.268482.120

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Su Q, Huang W, Huang Y, Dai R, Chang C, Li QY, Liu H, Li Z, Zhao Y, Wu Q, et al. Single-cell insights: pioneering an integrated atlas of chromatin accessibility and transcriptomic landscapes in diabetic cardiomyopathy. Cardiovasc Diabetol. 2024;23(1):139. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12933-024-02233-y

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Cai C, Wan P, Wang H, Cai X, Wang J, Chai Z, Wang J, Wang H, Zhang M, Yang N, et al. Transcriptional and open chromatin analysis of bovine skeletal muscle development by single-cell sequencing. Cell Prolif. 2023;56(9):e13430. https://doiorg.publicaciones.saludcastillayleon.es/10.1111/cpr.13430

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Ranzoni AM, Tangherloni A, Berest I, Riva SG, Myers B, Strzelecka PM, Xu J, Panada E, Mohorianu I, Zaugg JB, et al. Integrative Single-Cell RNA-Seq and ATAC-Seq analysis of human developmental hematopoiesis. Cell Stem Cell. 2021;28(3):472–e487477. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.stem.2020.11.015

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Ackermann AM, Wang Z, Schug J, Naji A, Kaestner KH. Integration of ATAC-seq and RNA-seq identifies human alpha cell and beta cell signature genes. Mol Metab. 2016;5(3):233–44. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.molmet.2016.01.002

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Kim S, Wysocka J. Deciphering the multi-scale, quantitative cis-regulatory code. Mol Cell. 2023;83(3):373–92. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.molcel.2022.12.032

    Article  CAS  PubMed  Google Scholar 

  66. Lv Z, Meng J, Yao S, Xiao F, Li S, Shi H, Cui C, Chen K, Luo X, Ye Y, et al. Naringenin improves muscle endurance via activation of the Sp1-ERRγ transcriptional axis. Cell Rep. 2023;42(11):113288. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.celrep.2023.113288

    Article  CAS  PubMed  Google Scholar 

  67. Zhuang ZH, Zhong Y, Chen YH, Zhang ZW. Research progress on the roles of Krüppel-like factors in muscle tissues. Yi Chuan. 2018;40(9):733–48. https://doiorg.publicaciones.saludcastillayleon.es/10.16288/j.yczz.18-095

    Article  PubMed  Google Scholar 

  68. Wang X, Shen QW, Wang J, Zhang Z, Feng F, Chen T, Zhang Y, Wei H, Li Z, Wang X, et al. KLF7 regulates satellite cell quiescence in response to extracellular signaling. Stem Cells. 2016;34(5):1310–20. https://doiorg.publicaciones.saludcastillayleon.es/10.1002/stem.2346

    Article  CAS  PubMed  Google Scholar 

  69. Zhao ZD, Zan LS, Li AN, Cheng G, Li SJ, Zhang YR, Wang XY, Zhang YY. Characterization of the promoter region of the bovine long-chain acyl-CoA synthetase 1 gene: roles of E2F1, Sp1, KLF15, and E2F4. Sci Rep. 2016;6:19661. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/srep19661

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Wang J, Chen T, Feng F, Wei H, Pang W, Yang G, Shen QW. KLF15 regulates slow myosin heavy chain expression through NFATc1 in C2C12 myotubes. Biochem Biophys Res Commun. 2014;446(4):1231–6. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.bbrc.2014.03.091

    Article  CAS  PubMed  Google Scholar 

  71. Wilm TP, Solnica-Krezel L. Essential roles of a zebrafish prdm1/blimp1 homolog in embryo patterning and organogenesis. Development. 2005;132(2):393–404. https://doiorg.publicaciones.saludcastillayleon.es/10.1242/dev.01572

    Article  CAS  PubMed  Google Scholar 

  72. Okumura F, Oki N, Fujiki Y, Ikuta R, Osaki K, Hamada S, Nakatsukasa K, Hisamoto N, Hara T, Kamura T. ZSWIM8 is a myogenic protein that partly prevents C2C12 differentiation. Sci Rep. 2021;11(1):20880. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/s41598-021-00306-6

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We thank the staff at our laboratory for their ongoing assistance.

Funding

This research was funded by the Sichuan Science and Technology Program (grant number: 2024NSFSC1167) and the Key R&D Projects of the Sichuan Provincial Department of Science and Technology (grant number: 2023YFN0045).

Author information

Author notes

  1. Kun Du and Dai-hua Wang contributed equally to this work.

Authors and Affiliations

  1. Technology Research Center of Modern Breeding Development, Mianyang Normal University, Mianyang, 621000, China

    Kun Du, Dai-hua Wang, Yu Xia, Qian Wu, Mao-qing Gu & Xi-wen Chen

  2. Farm Animal Genetic Resources Exploration and Innovation Key Laboratory of Sichuan Province, Sichuan Agricultural University, Chengdu, China

    Shen-qiang Hu

Authors
  1. Kun Du
    View author publications

    You can also search for this author inPubMed Google Scholar

  2. Dai-hua Wang
    View author publications

    You can also search for this author inPubMed Google Scholar

  3. Shen-qiang Hu
    View author publications

    You can also search for this author inPubMed Google Scholar

  4. Yu Xia
    View author publications

    You can also search for this author inPubMed Google Scholar

  5. Qian Wu
    View author publications

    You can also search for this author inPubMed Google Scholar

  6. Mao-qing Gu
    View author publications

    You can also search for this author inPubMed Google Scholar

  7. Xi-wen Chen
    View author publications

    You can also search for this author inPubMed Google Scholar

Contributions

Conceptualization, X.C. and K.D.; methodology, D.W.; software, S.H. and M.G.; validation, K.D. and Q.W.; formal analysis, K.D. and Y.X.; resources, D.W.; writing—original draft preparation, K.D.; writing—review and editing, X.C.; visualization, K.D.; supervision, X.C.; project administration, Y.X.; funding acquisition, K.D. and X.C. All authors have read and agreed to the published version of the manuscript.

Corresponding author

Correspondence to Xi-wen Chen.

Ethics declarations

Institutional review board statement

The study was conducted in accordance with the Declaration of Helsinki, and approved by the Institutional Review Board of the Biological Studies Animal Care and Use Committee, Sichuan Province, China (protocol code SKY101368, 4 November 2023).

Conflicts of interest

The authors declare no conflicts of interest.

Additional information

Publisher’s note

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

Electronic supplementary material

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Du, K., Wang, Dh., Hu, Sq. et al. Genome-wide chromatin accessibility and selective signals of meat rabbits reveal key Cis-regulatory elements and variants during postnatal development of skeletal muscles in rabbits. BMC Genomics 26, 296 (2025). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12864-025-11496-y

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12864-025-11496-y

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us