Skip to main content

Comprehensive analysis of lncRNAs involved in skeletal muscle development in ZBED6-knockout Bama pigs

Abstract

Background

The mutation of insulin-like growth factor 2 (IGF2 mutation) that a single-nucleotide substitution (G→A) in the third intron of IGF2 abrogates the interaction with zinc finger BED-type containing 6 (ZBED6) and leads to increased muscle mass in pigs. IGF2 mutation knock-in (IGF2 KI) and ZBED6 knockout (ZBED6 KO) lead to changes in IGF2 expression and increase muscle mass in mice and pigs. Long noncoding RNAs (lncRNAs) may participate in numerous biological processes, including skeletal muscle development. However, the role of the ZBED6-lncRNA axis in skeletal muscle development is poorly characterized.

Results

In this study, we assembled transcriptomes using RNA-seq data published in previous studies by our group and identified 11,408 known lncRNAs and 2269 potential lncRNAs in seven tissues, heart, longissimus dorsi, gastrocnemius muscle, liver, spleen, lung and kidney, of ZBED6 KO (lean mass model) and WT Bama pigs. ZBED6 affected the expression of 1570 lncRNAs (differentially expressed lncRNAs [DE-lncRNAs]; log2-fold change ≥ 1, nominal p-value ≤ 0.05) in the seven examined tissues. The expressed lncRNAs (FPKM > 0.1) exhibited tissue-specific patterns in WT pigs. Specifically, 3410 lncRNAs were expressed exclusively in only one tissue. Potential functions of lncRNAs were indirectly predicted by searching their target cis- and trans-regulated protein-coding genes. LncRNAs with tissue-specific expression influence numerous genes related to tissue functions. Weighted gene coexpression network analysis (WGCNA) of 1570 DE-lncRNAs between WT and ZBED6 KO pigs was used to define the following six lncRNA modules specific to different tissues: skeletal muscle, heart, lung, spleen, kidney and liver modules. Furthermore, by conjoint analysis of longissimus dorsi data (tissue-specific expression, muscle module and DE-lncRNAs) and ChIP-PCR revealed NONSUSG002145.1 (adjusted p-values = 0.044), which is coexpressed with the IGF2 gene and binding with ZBED6, may play important roles in ZBED6 KO pig skeletal muscle development.

Conclusions

These findings indicate that the identified lncRNAs may play essential roles in tissue function and regulate the mechanism of ZBED6 action in skeletal muscle development in pigs. To our knowledge, this is the first study describing lncRNAs in ZBED6 KO pigs. These results may open new research directions leading to a better understanding of the global functions of ZBED6 and of lncRNA functions in skeletal muscle development in pigs.

Peer Review reports

Background

With the rapid development of RNA sequencing (RNA-seq) technologies, tens of thousands of noncoding RNAs have been discovered [1, 2]. Long noncoding RNAs (lncRNAs) are defined as autonomously transcribed noncoding RNAs greater than 200 nt in length that do not overlap with annotated coding genes [3]. Previously, lncRNAs, which have lower expression and protein-coding potential than mRNAs, had long been regarded as transcription junk [4]. However, in recent decades, increasing evidence has indicated that lncRNAs play important roles in many biological processes, such as the regulation of skeletal muscle development [5, 6], cell fate decisions [7], and subcutaneous fat deposition [8].

Zinc finger BED domain containing protein 6 (ZBED6), which regulates skeletal muscle development through insulin-like growth factor 2 (IGF2), is a transcriptional repressor [9,10,11,12]. A single nucleotide transition from G to A in intron 3 of IGF2, a paternally expressed quantitative trait locus (QTL) in pigs, abrogates ZBED6-IGF2 binding and results in 3-fold greater postnatal expression of IGF2 mRNA in skeletal muscle, leading to increased muscle mass and heart size and reduced fat deposition in pigs. Recently, reports showed that ZBED6 regulates IGF2 mRNA expression in several human and murine cell lines and in mice [13,14,15,16,17]. Pigs are among the most essential commercially farmed animals and have attracted substantial attention in the context of muscle development. With long-term selection for lean meat content, the QTL in IGF2 is fixed in modern commercial pigs, such as Large White, Landrace and Hampshire [11], but most indigenous Chinese breeds homozygous for the IGF2 wild-type allele have a lower rate of lean meat [18]. IGF2 knock-in (IGF2 KI) and ZBED6 knockout (ZBED6 KO) leads to increased IGF2 expression and muscle mass in mice and pigs [19,20,21,22]. The functional role of ZBED6, in addition to its important role in regulating IGF2 expression, remains unclear. Hence, it is essential to elucidate the function of ZBED6 in pigs.

Recently, we obtained ZBED6 KO pigs that showed more muscle mass than wild-type (WT) pigs and analysed the effect of ZBED6 on gene expression [22], but the relationship of ZBED6 and lncRNAs was unclear. The ZBED6 KO pig provides a good lean mass pig model to illustrate expression profiles and the functions of lncRNAs in skeletal muscle development. Here, we performed RNA-seq analysis on seven different tissues of ZBED6 KO and WT pigs to systematically identify lncRNAs. Differentially expressed lncRNAs (DE-lncRNAs), tissue expression patterns and potential functions were elucidated. Our study will be of great use in future explorations of skeletal muscle lncRNA and ZBED6 functions.

Results

Identification and characteristics of lncRNAs expressed in different tissues of WT and ZBED6 KO pigs

To comprehensively investigate the lncRNAs in WT and ZBED6 KO pigs, the Illumina HiSeq 2500 platform was used to obtain a comprehensive view of lncRNAs in seven tissues, namely, muscle (gastrocnemius muscle and longissimus dorsi), heart, liver, spleen, lung and kidney, between ZBED6 KO and WT pigs, with three biological replicates for each group (a total of 42 samples). We generated approximately 4.2 to 8.74 G of raw reads, and 70 %~87.3 % of the clean reads mapped to the genome of the paired reads for each sample (Supplementary Table 2). The GC content of the clean data ranged from 47.05 to 57.92 %, demonstrating that the reliability and quality of the sequencing data were adequate for further analysis. To identify lncRNAs, the transcripts were assembled and reconstructed into a total of 210,858 transcripts using Cufflinks software. Ultimately, in addition to 11,408 known lncRNAs, 4352 potential lncRNA transcripts that overlapped 2269 lncRNA genes were identified by using CNCI software (Fig. 1 A-B). It is important to maximize our understanding of lncRNAs through multiple structural features. We found that 48.59 % of the lncRNAs were more than 2000 bp (Fig. 1 C), whereas 85.39 % spanned two or more exons (Fig. 1D). LncRNAs tended to have lower expression levels than mRNAs (Fig. 1E). The features of lncRNAs were consistent with those of lncRNAs reported in other studies.

Fig. 1
figure 1

Statistics of lncRNAs identified in this study. (A) Percentage of known lncRNA transcripts. (B) Percentage of potentially novel transcripts. (C) Length distribution of lncRNAs. (D) Exon number of lncRNAs. (E) Expression levels of mRNAs and lncRNAs

Principal component analysis (PCA) of these 42 samples clearly showed that the lncRNA expression profiles of seven tissues were clustered together separately, and the profiles were similar in skeletal muscle (gastrocnemius muscle and longissimus dorsi) (Fig. 2 A). We detected 13,677 expressed lncRNAs in this experiment Using log2-fold chang ≥ 1, nominal p-value ≤ 0.05 as the filter criteria,we found 127 to 295 DE-lncRNAs across seven tissues between the WT and ZBED6 KO samples (Fig. 2B), of which 177 and 293 DE-lncRNAs were identified in gastrocnemius muscle and longissimus dorsi, respectively (Fig. 2 C-D). Using log2-fold change ≥ 1, adjusted p-values ≤ 0.05 as the filter criteria, 82 to 211 DE-lncRNAs were detected, of which 116 and 211 DE-lncRNAs were found in gastrocnemius muscle and longissimus dorsi, respectively (Supplementary Table 3). The details of DE-lncRNAs in the seven tissues are shown in Supplementary Table 3. These results suggest that ZBED6 regulates the expression of lncRNAs in pigs.

Fig. 2
figure 2

Differentially expressed lncRNAs. (A) PCA of the expressed lncRNAs in multiple tissues in WT and ZBED6 −/− pigs (n = 3). (B) The number of DE-lncRNAs in each tissue. (C-D) Volcano plot showing 293 and 177 DE-lncRNAs of longissimus dorsi (LD) and gastrocnemius muscle (GM) between WT and ZBED6−/− pigs (n = 3)

Tissue-specific expression of lncRNAs in WT pig

To elucidate the functional roles of lncRNAs in WT pigs, we used lncRNA expression data (FPKM > 0.1) from the seven different tissues from three WT pigs to characterize tissue-specific expression patterns. A heatmap based on the expression patterns in all tissues was generated across all expressed lncRNAs, which indicated that the majority exhibited a preferential tissue expression pattern (Fig. 3 A). A total of 488 commonly expressed lncRNAs were found in all tissues, and 3410 lncRNAs displayed exclusive expression in one tissue (Supplementary Table 4), with the most tissue-specific lncRNAs present in the longissimus dorsi (N = 1177) and the fewest tissue-specific lncRNAs present in the gastrocnemius muscle (N = 258) (Fig. 3B). To analyse the possible roles of these lncRNAs in pig tissues, potential cis- and trans-target genes for lncRNAs were predicted. As a result, a total of 1815 cis- and 1804 trans-target genes for common and tissue-specific expressed lncRNAs were detected (Supplementary Table 5). Target genes of the 488 commonly expressed lncRNAs, including adiponectin (ADIPOQ), peroxisomal proliferating-activated receptor gamma (PPARG), anabolic-androgenic steroids (AASS) and WNT9B, play central roles in the regulation of growth and fat development in pigs and were enriched in 20 KEGG pathways, including biosynthesis of antibiotics, the peroxisome proliferator-activated receptor (PPAR) signalling pathway and the Wnt signalling pathway (Fig. 3 C). Additionally, tissue-specific lncRNAs of the seven tissues were enriched in metabolic pathways, the intestinal immune network for IgA production, fatty acid metabolism and transcriptional misregulation in cancer (Supplementary Table 6).

Fig. 3
figure 3

Tissue-specific expression of lncRNAs in WT pig. (A) Heatmap of expressed lncRNAs across 7 tissues. Each row represents the expression levels of all detected lncRNAs, and each column contains all expressed transcripts. We transformed the FPKM values into log2(FPKM + 1) values and then calculated the Z-score for every log 2 (FPKM + 1) value within each tissue. (B) Venn diagram of expressed lncRNAs in seven tissues. (C) Enriched gene pathway terms of targeted genes of 488 commonly expressed lncRNAs across 7 tissues [23]

Tissue expression patterns of lncRNAs in WT and ZBED6 KO pigs

To demonstrate the dynamics of lncRNA expression in different tissues of WT and ZBED6 KO pigs, the expression patterns of DE-lncRNAs (nominal p-value ≤ 0.05) were clustered using weighted gene coexpression network analysis (WGCNA). A total of 7 lncRNA transcriptional modules were identified in WT and ZBED6 KO pig tissues, with six modules showing a strong correlation with specific tissues (Fig. 4 A-B). We focused on modules that were highly positively correlated with tissues (correlation > 0.3, p-value < 0.05) and defined six tissue-specific modules (Supplementary Table 7). The magenta module, which included 163 lncRNAs strongly related to the four muscle tissues (gastrocnemius muscle and longissimus dorsi of WT and ZBED6 KO pigs), was identified as a skeletal muscle module, whereas grey indicated weak correlation with all tissues (Fig. 4B). The yellow, brown, turquoise, blue and black modules showed continuously positive correlations with heart, lung, spleen, kidney and liver tissues in WT and ZBED6 KO pigs, respectively.

Fig. 4
figure 4

Expression modules of lncRNAs in WT and ZBED6 KO pigs determined by WGCNA. (A) Hierarchical cluster tree of all differentially expressed lncRNA modules. Modules correspond to the branch and are marked by a colour, as indicated by colour strips under the tree. (B) Correlation analysis between lncRNA modules and tissues

LncRNAs play a potential role in skeletal muscle development in ZBED6 KO pigs

To identify the function of lncRNAs in skeletal muscle development in ZBED6 KO pigs, we overlapped the lncRNAs from three analyses: longissimus dorsi tissue-specific expression, skeletal muscle module (magenta module) and DE-lncRNAs (nominal p-value ≤ 0.05) of the longissimus dorsi. Thirty common lncRNAs were revealed (Fig. 5 A), of which eleven lncRNAs were novel and nineteen lncRNAs were known. However, no common lncRNAs of gastrocnemius muscle were revealed with a similar analysis. Therefore, the study mostly was focused on longissimus dorsi. Heatmap analysis further showed that the expression of ten of the thirty common lncRNAs was upregulated in the longissimus dorsi of ZBED6 KO pigs (Fig. 5B). These results suggested that lncRNAs play an important role in skeletal muscle development in ZBED6 KO pigs. To explore how lncRNAs work in concert with their target genes (mRNAs) to regulate skeletal muscle development and to identify the key molecules, possible regulatory networks of interactions between lncRNAs and their target genes (mRNAs) were constructed. In the present study, we predicted potential cis- and trans-target genes of DE-lncRNAs and further constructed a lncRNA-gene interaction network between the DE-lncRNAs and their corresponding differentially expressed cis- and trans-target genes by using Cytoscape. For the thirty common lncRNAs in longissimus dorsi, we detected 21 cis and 268 trans target genes (Supplementary Table 8) and performed KEGG enrichment analysis to reveal the overrepresented biological processes of the 289 target genes using DAVID (Fig. 5 C). The KEGG results showed that the genes were related to the myocyte proliferation or differentiation pathway, the AMPK signalling pathway, the FoxO signalling pathway, factor-regulated calcium reabsorption and others (Supplementary Table 9). In addition, metabolic pathways, the citrate cycle (TCA cycle), glucagon signalling pathway and pyruvate metabolism were overrepresented. The lncRNA-mRNA and ChIP-PCR analysis showed that among the thirty common lncRNAs, NONSUSG002145.1 that was upregulated in the longissimus dorsi of ZBED6 KO pigs (adjusted p-values = 0.044), is the target of ZBED6 and perfectly coexpressed with a crucial regulator of skeletal muscle development, IGF2 (Fig. 6 A-B), suggesting that NONSUSG002145.1 could be potential regulators of IGF2 expression during skeletal muscle development.

Fig. 5
figure 5

LncRNAs in longissimus dorsi (LD) development in ZBED6 KO pigs. (A) Thirty common lncRNAs in longissimus dorsi tissue-specific expression (tissue-specific), the skeletal muscle module (magenta module) and DE-lncRNAs of longissimus dorsi (LD_DE-lncRNAs). (B) Heatmap with log2(fold-change) between WT and ZBED6 KO pig longissimus dorsi of 30 common lncRNAs. (C) Enriched gene pathway terms of 30 common lncRNA-targeted genes [23]

Fig. 6
figure 6

Functional network of the key lncRNA in longissimus dorsi (LD) development in ZBED6 KO pigs. (A) Function network of NONSUSG002145.1. (B) ChIP-PCR analysis in longissimus dorsi of ZBED6 occupancy at the binding sites of NONSUSG002145.1. The red triangle represents lncRNA, the circle represents target mRNA of lncRNA, the blue circle represents target mRNA IGF2 related to muscle development. (C) qRT-PCR validation

RT-qPCR Validation of DE-lncRNAs

To validate the RNA sequencing data, 9 lncRNAs were randomly selected from DE-lncRNAs, and their expression levels in different tissues of WT and ZBED6 KO pigs were examined by RT-qPCR. The results showed that the expression of the 9 lncRNAs was consistent with the expression trends calculated from the RNA-seq data (Fig. 6 C).

Discussion

Thousands of lncRNAs have been identified and implicated in reprogramming [24], embryonic development [25, 26], imprinting contro l[27], cell differentiation and development [28, 29], and skeletal muscle growth [30, 31] using RNA-seq analyses, suggesting that RNA-seq is highly effective for the discovery of lincRNAs. Studies targeting lncRNAs in humans and mice have become a research hotspot in recent years. There are many excellent reviews with a specific focus on these functional and regulatory issues [3, 28, 32,33,34]. Compared to research on lncRNAs in humans and mice, studies on lncRNAs are still in their infancy in pigs.

ZBED6 is a transcriptional repressor and involved in skeletal muscle growth by regulating IGF2 expression [9,10,11,12,13,14,15,16]. Recently, pigs and mice with genetically modified IGF2 and ZBED6 showed increased muscle mass and IGF2 expression [17, 20,21,22]. The functional role of ZBED6 in relation to lncRNAs is still poorly characterized. Therefore, we obtained a lean mass model pig by knocking out ZBED6 and generated 545 G of RNA-seq data, which were made publicly available [22]. Using RNA-seq data we identified more than 13,000 lncRNAs in seven tissues. These lncRNAs showed a shorter length, fewer exons and lower expression levels than mRNAs, which is consistent with findings in humans [35], mice [36], sheep [29, 30] and horses [37]. Almost 83.41 % of the identified lncRNAs (11,408) in our study were previously reported in a pig RNA-seq study [5]. The identification of potentially novel lncRNAs might be due to differences in various aspects of these studies, such as types of pigs and sampled tissue. In addition, the RNA-seq analysis and RT-qPCR validation indicated that the lncRNAs are bona fide transcripts rather than random products from transcriptional noise.

Several studies have already revealed that lncRNAs tend to show tissue-specific expression [38, 39]. Studies in humans [35] and sheep [29] have shown that tissue-specific lncRNAs are involved in certain biological processes in specific tissues. Target gene analysis of lncRNAs indicated high correlation between lncRNAs and mRNAs, consistent with previous findings showing that lncRNAs were likely to be functionally associated with their nearest neighbouring mRNAs [40]. That is, the functions of lncRNAs can be predicted through their target mRNAs. We found 488 commonly expressed lncRNAs across the seven tissues in WT pigs, including 1213 lncRNA-mRNA pairs that regulate growth and fat development. ADIPOQ and PPARG, which are involved in fat metabolism pathways (PPAR signalling pathway and adipocytokine signalling pathway), play a crucial role in adipogenic induction [19, 41]. WNT9B, a secreted Wnt growth factor enriched in the development and growth pathway (Wnt signalling pathway), is a developmental regulator [42]. This suggested that lncRNAs have various roles in many biological processes in pigs, such as muscle development, metabolism and immunity, consistent with previous studies [26, 43,44,45,46]. These data provide a high-quality lncRNA resource for pigs.

293 DE-lncRNAs (nominal p-value ≤ 0.05) were found in the muscles of WT and ZBED6 KO pigs, of which the adjusted p-values of 211 DE-lncRNAs was less than 0.05. And 163 specific lncRNAs were observed in the muscles of WT and ZBED6 KO pigs, suggesting that ZBED6 is associated with the growth of pig tissues by regulating lncRNAs. Interestingly, we identified NONSUSG002145.1 (adjusted p-values = 0.044) which passed multiple testing or FDR correction, coexpressed with IGF2, could bind with ZBED6 directly. Due to the central role of IGF2 in skeletal muscle development in mice and pigs [9, 10, 47], NONSUSG002145.1 that was reported rarely may have regulatory roles in the expression of IGF2 during the development of muscle in ZBED6 KO pigs. Therefore, additional studies, such as loss-of-function experiments, are needed in the future to provide further insights into the regulatory functions of NONSUSG002145.1 in skeletal muscle development.

Conclusions

In the current study, we identified and characterized lncRNAs in WT pigs. LncRNAs with tissue-specific expression across seven tissues in WT pigs were discovered. Some lncRNAs expressed across all tissues tested may influence the expression of numerous genes, including those involved in (1) fat metabolism factors (e.g., ADIPOQ and PPARG) and (2) development and growth factors (e.g., WNT9B). We also examined the effects of ZBED6 on the lncRNA expression profile and identified DE-lncRNAs between WT and ZBED6 KO pigs. WGCNA revealed lncRNA tissue expression patterns in WT and ZBED6 KO pigs. NONSUSG002145.1 (adjusted p-values = 0.044), regulating the expression of IGF2, were identified as key molecules related to muscle development in ZBED6 KO pigs.

Methods

Ethics statement

All animal experiments were performed according to protocols and guidelines approved by the Institutional Animal Care and Use Committee of the Beijing Academy of Agricultural Sciences.

Datasets

All RNA-seq data of female ZBED6 KO and female WT Bama miniature pigs used in this study were from previous studies in our laboratory that have been deposited in the Sequence Read Archive (https://www.ncbi.nlm.nih.gov) with the indicated accession codes (BioProject ID: PRJNA663759). The RNA-seq data included seven tissues (heart, longissimus dorsi, gastrocnemius muscle, liver, spleen, lung and kidney) from three ZBED6 KO and three WT Bama pigs (Supplementary Table 10).

Sequence analysis

Raw reads were cleaned by removing adapter sequences, reads with over 10 % N sequences, and low-quality reads in which the number of bases with a quality value Q ≤ 10 was more than 50 %, and then clean reads were obtained. The Q20, Q30, and GC contents of the clean data were calculated. The clean reads were mapped to the pig reference genome (https://www.ncbi.nlm.nih.gov/genome/84?genome_assembly_id=317145) using TopHat (v2.0.11) with default parameters [48]. The BAM files generated from TopHat were assembled, and transcripts were constructed by Cufflinks (v2.2.1) [48]. Then, lncRNA analyses were performed using the lncRNA calling protocol and CNCI [49, 50]. RNA length ≥ 200 nt, exon number = > 2, CNCI score = > 0 were used to evaluate the coding potential of transcripts. We applied Cuffdiff (part of Cufflinks) to detect DElncRNAs between the ZBED6 KO and WT pig. Using the negative binomial distribution, p-value was calculated to estimate changes in a transcript’s fragment count. Then p-value was corrected multiply to calculate q-value (adjusted p-values) by false discovery rate (FDR) and the Storey’s q-value procedure [51]. Transcripts with fold changes ≥ 2.0 and nominal p-value ≤ 0.05 were identified as DElncRNAs. Cufflinks and Cuffdiff were used as described by Trapnell et al [48].

PCA and functional enrichment analysis

The FPKM values for all of the annotated transcripts from the seven tissue transcriptomes were used to perform PCA, which was implemented by gmodels in R (version 3.1.3, http://cran.r-project.org/). GO enrichment analysis of differentially expressed genes was implemented using DAVID (https://david.ncifcrf.gov/). Pig genes were used as the background gene set when performing the KEGG enrichment analysis [23]. P values ≤ 0.05 were considered significant.

Tissue-specific identification of lncRNAs

We used RNA-seq data of seven tissues (heart, longissimus dorsi, gastrocnemius muscle, liver, spleen, lung and kidney) from WT pigs to characterize the expression patterns of the lncRNA genes. The heatmap was generated using all expressed lncRNA genes in seven tissues using the pheatmap package in R based on the FPKM values (FPKM > 0.1). The level of gene expression is visualized by using a colour gradient from blue to red. For seven tissues from WT pigs, lncRNAs were classified into tissue-specific lncRNAs that were observed in only one tissue.

WGCNA Network Analysis

The genes that were differentially expressed between all stages of seven tissues from WT and KO pigs were selected to identify key genes in tissue development using the R package WGCNA [52]. A signed weighted correlation network was constructed by first creating a matrix of Pearson correlation coefficients between all pairs of genes across the measured samples. The adjacency matrix was then transformed into a topological overlap matrix (TOM) to minimize the effects of noise and spurious associations. To define modules as branches, we employed the dynamic tree cut algorithm with default parameters to cut the hierarchal clustering tree [53].

Target Gene Prediction

A previous study showed that lncRNAs exert cis-regulatory effects on their colocalized genes [54]. The colocalization role is lncRNA acting on neighbouring target genes. The 10 k coding genes upstream and downstream of lncRNAs were searched to identify their cis-acting effects. Additionally, lncRNAs recognize targeted genes at the expression level. The Pearson correlation coefficient method was used to analyse the correlation between lncRNAs and mRNAs. mRNAs with an absolute correlation value greater than 0.9 and a p-value less than 0.05 were identified as trans target genes.

Chromatin immunoprecipitation (ChIP)-PCR analysis

ChIP was conducted previously described [55]. Briefly, freshly longissimus dorsi were treated with 1 % formaldehyde for 18 min and neutralizing with glycine (AMRESCO, USA) for 5 min at room temperature. Then longissimus dorsi was smashed and resuspended in SDS lysis buffer (Beyotime, China). After incubation for 20 min at 4 °C, the lysates were sonicated 24 times (30 s each) (Bioruptor® Sonication System, USA). Total solution and Protein G agarose beads (ThermoFisher Scientific, USA) were incubated at 4℃ overnight with 10ug of the following antibodies: anti-ZBED6 antibody (HPA068807, ATLAS) and normal mouse IgG antibody (2729 S, CST). The beads were washed to get ChIP Elution Buffer. Add 20 µl 5 M NaCl to the combined eluates and reverse histone-DNA crosslinks by heating at °65 for 4 h. RNase A (TIANGEN, China) and Proteinase K (TIANGEN, China) were added to ChIP Elution Buffer to remove RNA and proteins. ChIP DNA Clean &Concentrator purification spin column (ZYMO, Irvine, California, USA) was used to purify coprecipitated DNA. The ZBED6 binding site, NONSUSG002145.1 was evaluated using PCR and normalized by total chromatin (input). Normal mouse IgG was used as the negative control, the primers are described in supplementary Table 1. The PCR conditions were: 1 cycle at 95℃ for 5 min followed by 33 cycles at 95℃ for 30 s, 58–59℃ for 30 s, and 72℃ for 30 s. The PCR products were then electrophoresed on 1.5 % agarose gels stained with GelGreen (TIANGEN, China).

Quantitative PCR

For the qRT-PCR analysis of genes, single-strand cDNA was synthesized from RNA using the PrimeScript RT reagent Kit With gDNA Eraser (Takara, China). The qPCR reactions was performed in ABI MicroAmp Optical 96-well Reaction plates on an ABI 7500 real-time PCR instrument (USA), using SYBR Premix Ex Taq (Tli RNaseH Plus; Takara). The primer sequences, designed using Premier Primer 5.0 software, are listed in Supplementary Table 1. Relative gene expression was calculated using the comparative cycle threshold (2−ΔΔCt) method.

Statistical analysis

Statistical analyses of the qRT-PCR results and graphs were carried out in GraphPad Prism 7 (GraphPad Software Inc, La Jolla, CA, USA). Statistical significance of the data was tested by performing paired t-tests. The P values were determined by Student’s t test (unpaired). Error bars indicate SEM. ***, P < 0.001; **, P < 0.01; *, P < 0.05.

Availability of data and materials

All the RNA-seq reads have been deposited in the Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra) with the accession codes (BioProject ID: PRJNA663759).

Abbreviations

RNA-Seq:

RNA sequencing

lncRNA:

Long non-coding RNA

ZBED6:

Zinc finger BED domain containing protein 6

IGF2:

Insulin-like growth factor 2

QTL:

Quantitative trait locus

DE-lncRNAs:

differentially expressed lncRNAs

ZBED6 KO:

ZBED6 knockout

WT:

Wildtype

WGCNA:

Weighted gene co-expression network analysis

GM:

Gastrocnemius muscle

LD:

Longissimus dorsi

TOM:

Topological overlap matrix

PCA:

Principal component analysis

ADIPOQ:

Adiponectin

PPARG:

Peroxisomal proliferating-activated receptor gamma

AASS:

Anabolic-androgenic steroids

PPAR:

Peroxisome proliferators-activated receptors

ChIP:

Chromatin immunoprecipitation

References

  1. Harrow J, Frankish A, Gonzalez JM, Tapanari E, Diekhans M, Kokocinski F, Aken BL, Barrell D, Zadissa A, Searle S et al: GENCODE: the reference human genome annotation for The ENCODE Project. Genome research 2012, 22(9):1760–1774.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  2. Guttman M, Amit I, Garber M, French C, Lin MF, Feldser D, Huarte M, Zuk O, Carey BW, Cassady JP et al: Chromatin signature reveals over a thousand highly conserved large non-coding RNAs in mammals. Nature 2009, 458(7235):223–227.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  3. Quinn JJ, Chang HY: Unique features of long non-coding RNA biogenesis and function. Nat Rev Genet 2016, 17(1):47–62.

    CAS  PubMed  Article  Google Scholar 

  4. Huttenhofer A, Schattner P, Polacek N: Non-coding RNAs: hope or hype? Trends Genet 2005, 21(5):289–297.

    PubMed  Article  CAS  Google Scholar 

  5. Zhao W, Mu Y, Ma L, Wang C, Tang Z, Yang S, Zhou R, Hu X, Li MH, Li K: Systematic identification and characterization of long intergenic non-coding RNAs in fetal porcine skeletal muscle development. Sci Rep 2015, 5:8957.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  6. Yang Y, Liang G, Niu G, Zhang Y, Zhou R, Wang Y, Mu Y, Tang Z, Li K: Comparative analysis of DNA methylome and transcriptome of skeletal muscle in lean-, obese-, and mini-type pigs. Sci Rep 2017, 7:39883.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  7. Joo MS, Shin SB, Kim EJ, Koo JH, Yim H, Kim SG: Nrf2-lncRNA controls cell fate by modulating p53-dependent Nrf2 activation as an miRNA sponge for Plk2 and p21(cip1). FASEB J 2019, 33(7):7953–7969.

  8. Liu X, Liu K, Shan B, Wei S, Li D, Han H, Wei W, Chen J, Liu H, Zhang L: A genome-wide landscape of mRNAs, lncRNAs, and circRNAs during subcutaneous adipogenesis in pigs. J Anim Sci Biotechnol 2018, 9:76.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  9. Jeon JT, Carlborg O, Tornsten A, Giuffra E, Amarger V, Chardon P, Andersson-Eklund L, Andersson K, Hansson I, Lundstrom K et al: A paternally expressed QTL affecting skeletal and cardiac muscle mass in pigs maps to the IGF2 locus. Nat Genet 1999, 21(2):157–158.

    CAS  PubMed  Article  Google Scholar 

  10. Nezer C, Moreau L, Brouwers B, Coppieters W, Detilleux J, Hanset R, Karim L, Kvasz A, Leroy P, Georges M: An imprinted QTL with major effect on muscle mass and fat deposition maps to the IGF2 locus in pigs. Nat Genet 1999, 21(2):155–156.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  11. Van Laere AS, Nguyen M, Braunschweig M, Nezer C, Collette C, Moreau L, Archibald AL, Haley CS, Buys N, Tally M et al: A regulatory mutation in IGF2 causes a major QTL effect on muscle growth in the pig. Nature 2003, 425(6960):832–836.

    PubMed  Article  CAS  PubMed Central  Google Scholar 

  12. Markljung E, Jiang L, Jaffe JD, Mikkelsen TS, Wallerman O, Larhammar M, Zhang X, Wang L, Saenz-Vash V, Gnirke A et al: ZBED6, a novel transcription factor derived from a domesticated DNA transposon regulates IGF2 expression and muscle growth. PLoS Biol 2009, 7(12):e1000256.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  13. Jiang L, Wallerman O, Younis S, Rubin CJ, Gilbert ER, Sundstrom E, Ghazal A, Zhang X, Wang L, Mikkelsen TS et al: ZBED6 modulates the transcription of myogenic genes in mouse myoblast cells. PLoS One 2014, 9(4):e94187.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  14. Akhtar Ali M, Younis S, Wallerman O, Gupta R, Andersson L, Sjoblom T: Transcriptional modulator ZBED6 affects cell cycle and growth of human colorectal cancer cells. Proc Natl Acad Sci U S A 2015, 112(25):7743–7748.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  15. Wang X, Jiang L, Wallerman O, Engstrom U, Ameur A, Gupta RK, Qi Y, Andersson L, Welsh N: Transcription factor ZBED6 affects gene expression, proliferation, and cell death in pancreatic beta cells. Proc Natl Acad Sci U S A 2013, 110(40):15997–16002.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  16. Wang X, Jiang L, Wallerman O, Younis S, Yu Q, Klaesson A, Tengholm A, Welsh N, Andersson L: ZBED6 negatively regulates insulin production, neuronal differentiation, and cell aggregation in MIN6 cells. FASEB J 2019, 33(1):88–100.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  17. Younis S, Schönke M, Massart J, Hjortebjerg R, Sundström E, Gustafson U, Björnholm M, Krook A, Frystyk J, Zierath JR et al: The ZBED6-IGF2 axis has a major effect on growth of skeletal muscle and internal organs in placental mammals. Proc Natl Acad Sci U S A 2018, 115(9):E2048-e2057.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  18. Yang GC, Ren J, Guo YM, Ding NS, Chen CY, Huang LS: Genetic evidence for the origin of an IGF2 quantitative trait nucleotide in Chinese pigs. Anim Genet 2006, 37(2):179–180.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  19. Achari AE, Jain SK: Adiponectin, a Therapeutic Target for Obesity, Diabetes, and Endothelial Dysfunction. Int J Mol Sci 2017, 18(6).

  20. Xiang G, Ren J, Hai T, Fu R, Yu D, Wang J, Li W, Wang H, Zhou Q: Editing porcine IGF2 regulatory element improved meat production in Chinese Bama pigs. Cell Mol Life Sci 2018, 75(24):4619–4628.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  21. Liu X, Liu H, Wang M, Li R, Zeng J, Mo D, Cong P, Liu X, Chen Y, He Z: Disruption of the ZBED6 binding site in intron 3 of IGF2 by CRISPR/Cas9 leads to enhanced muscle development in Liang Guang Small Spotted pigs. Transgenic Res 2019, 28(1):141–150.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  22. Wang D: Studying the Molecular Mechanism Underlying Lean Meat Trait by Using ZBED6-Knockout Pig. In: The 3rd China Animal Husbandry Genome Industry Transformation Summit Forum Paper Abstracts. Beijing, China; 2017.

  23. Kanehisa M, Furumichi M, Sato Y, Ishiguro-Watanabe M, Tanabe M: KEGG: integrating viruses and cellular organisms. Nucleic Acids Res 2021, 49(D1):D545-D551.

    PubMed  Article  CAS  PubMed Central  Google Scholar 

  24. Loewer S, Cabili MN, Guttman M, Loh YH, Thomas K, Park IH, Garber M, Curran M, Onder T, Agarwal S et al: Large intergenic non-coding RNA-RoR modulates reprogramming of human induced pluripotent stem cells. Nat Genet 2010, 42(12):1113–1117.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  25. Pauli A, Valen E, Lin MF, Garber M, Vastenhouw NL, Levin JZ, Fan L, Sandelin A, Rinn JL, Regev A et al: Systematic identification of long noncoding RNAs expressed during zebrafish embryogenesis. Genome research 2012, 22(3):577–591.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  26. Ulitsky I, Shkumatava A, Jan CH, Sive H, Bartel DP: Conserved function of lincRNAs in vertebrate embryonic development despite rapid sequence evolution. Cell 2011, 147(7):1537–1550.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  27. Lee JT, Bartolomei MS: X-inactivation, imprinting, and long noncoding RNAs in health and disease. Cell 2013, 152(6):1308–1323.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  28. Wapinski O, Chang HY: Long noncoding RNAs and human disease. Trends Cell Biol 2011, 21(6):354–361.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  29. Song S, Yang M, Li Y, Rouzi M, Zhao Q, Pu Y, He X, Mwacharo JM, Yang N, Ma Y et al: Genome-wide discovery of lincRNAs with spatiotemporal expression patterns in the skin of goat during the cashmere growth cycle. BMC Genomics 2018, 19(1):495.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  30. Ling Y, Zheng Q, Sui M, Zhu L, Xu L, Zhang Y, Liu Y, Fang F, Chu M, Ma Y et al: Comprehensive Analysis of LncRNA Reveals the Temporal-Specific Module of Goat Skeletal Muscle Development. Int J Mol Sci 2019, 20(16).

  31. Kallen AN, Zhou XB, Xu J, Qiao C, Ma J, Yan L, Lu L, Liu C, Yi JS, Zhang H et al: The imprinted H19 lncRNA antagonizes let-7 microRNAs. Mol Cell 2013, 52(1):101–112.

    CAS  PubMed  Article  Google Scholar 

  32. Rinn JL, Chang HY: Genome regulation by long noncoding RNAs. Annu Rev Biochem 2012, 81:145–166.

    CAS  PubMed  Article  Google Scholar 

  33. Heward JA, Lindsay MA: Long non-coding RNAs in the regulation of the immune response. Trends Immunol 2014, 35(9):408–419.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  34. Gloss BS, Dinger ME: The specificity of long noncoding RNA expression. Biochim Biophys Acta 2016, 1859(1):16–22.

    CAS  PubMed  Article  Google Scholar 

  35. Verma A, Jiang Y, Du W, Fairchild L, Melnick A, Elemento O: Transcriptome sequencing reveals thousands of novel long non-coding RNAs in B cell lymphoma. Genome Med 2015, 7:110.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  36. Lv J, Huang Z, Liu H, Liu H, Cui W, Li B, He H, Guo J, Liu Q, Zhang Y et al: Identification and characterization of long intergenic non-coding RNAs related to mouse liver development. Mol Genet Genomics 2014, 289(6):1225–1235.

    CAS  PubMed  Article  Google Scholar 

  37. Pu Y, Zhang Y, Zhang T, Han J, Ma Y, Liu X: Identification of Novel lncRNAs Differentially Expressed in Placentas of Chinese Ningqiang Pony and Yili Horse Breeds. Animals (Basel) 2020, 10(1).

  38. Zhou ZY, Li AM, Adeola AC, Liu YH, Irwin DM, Xie HB, Zhang YP: Genome-wide identification of long intergenic noncoding RNA genes and their potential association with domestication in pigs. Genome Biol Evol 2014, 6(6):1387–1392.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  39. Cabili MN, Trapnell C, Goff L, Koziol M, Tazon-Vega B, Regev A, Rinn JL: Integrative annotation of human large intergenic noncoding RNAs reveals global properties and specific subclasses. Genes Dev 2011, 25(18):1915–1927.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  40. Huarte M, Guttman M, Feldser D, Garber M, Koziol MJ, Kenzelmann-Broz D, Khalil AM, Zuk O, Amit I, Rabani M et al: A large intergenic noncoding RNA induced by p53 mediates global gene repression in the p53 response. Cell 2010, 142(3):409–419.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  41. Rosen ED, Walkey CJ, Puigserver P, Spiegelman BM: Transcriptional regulation of adipogenesis. Genes Dev 2000, 14(11):1293–1307.

    CAS  PubMed  PubMed Central  Google Scholar 

  42. Qian J, Jiang Z, Li M, Heaphy P, Liu YH, Shackleford GM: Mouse Wnt9b transforming activity, tissue-specific expression, and evolution. Genomics 2003, 81(1):34–46.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  43. Wang Y, Xue S, Liu X, Liu H, Hu T, Qiu X, Zhang J, Lei M: Analyses of Long Non-Coding RNA and mRNA profiling using RNA sequencing during the pre-implantation phases in pig endometrium. Sci Rep 2016, 6:20238.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  44. Zhao Z, Bai J, Wu A, Wang Y, Zhang J, Wang Z, Li Y, Xu J, Li X: Co-LncRNA: investigating the lncRNA combinatorial effects in GO annotations and KEGG pathways based on human RNA-Seq data. Database (Oxford) 2015, 2015.

  45. Satpathy AT, Chang HY: Long noncoding RNA in hematopoiesis and immunity. Immunity 2015, 42(5):792–804.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  46. Sohi G, Dilworth FJ: Noncoding RNAs as epigenetic mediators of skeletal muscle regeneration. Febs j 2015, 282(9):1630–1646.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  47. Zanou N, Gailly P: Skeletal muscle hypertrophy and regeneration: interplay between the myogenic regulatory factors (MRFs) and insulin-like growth factors (IGFs) pathways. Cell Mol Life Sci 2013, 70(21):4117–4130.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  48. Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L: Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc 2012, 7(3):562–578.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  49. Luo H, Bu D, Sun L, Fang S, Liu Z, Zhao Y: Identification and function annotation of long intervening noncoding RNAs. Brief Bioinform 2017, 18(5):789–797.

    CAS  PubMed  PubMed Central  Google Scholar 

  50. Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, Liu Y, Chen R, Zhao Y: Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res 2013, 41(17):e166.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  51. Storey JD, Tibshirani R: Statistical significance for genomewide studies. Proc Natl Acad Sci U S A 2003, 100(16):9440–9445.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  52. Langfelder P, Horvath S: WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 2008, 9:559.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  53. Langfelder P, Zhang B, Horvath S: Defining clusters from a hierarchical cluster tree: the Dynamic Tree Cut package for R. Bioinformatics 2008, 24(5):719–720.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  54. Bonasio R, Shiekhattar R: Regulation of transcription by long noncoding RNAs. Annu Rev Genet 2014, 48:433–455.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  55. Lindroth AM, Park YJ, McLean CM, Dokshin GA, Persson JM, Herman H, Pasini D, Miro X, Donohoe ME, Lee JT et al: Antagonism between DNA and H3K27 methylation at the imprinted Rasgrf1 locus. PLoS Genet 2008, 4(8):e1000145.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

Download references

Funding

The project was supported by the National Key R & D Program of China (2020YFA0509500), the Agricultural Science and Technology Innovation Program of China (ASTIP-IAS01) and the National High-tech R&D Program of China (2017YFC1103702, 2017YFC1103700).

Author information

Authors and Affiliations

Authors

Contributions

LJ conceived the project and designed the research; DP generated and expanded the KO pig populations; DW, SW and WT collected and prepared all the tissues; LJ, DW, and YL contributed to functional validation and data analysis; DW, LJ, YP, and YM wrote the paper based on input from all the authors. The author(s) read and approved the final manuscript.

Corresponding authors

Correspondence to Yuehui Ma or Lin Jiang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

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

Supplementary Information

Additional file 1.

Supplementary Table S1. Primer pairs of DEGs used for qRT-PCR validation.

Additional file 2.

Supplementary Table S2. Quality analysis and genome mapping analysis of transcriptome sequencing.

Additional file 3.

Supplementary Table S3. List of DE-lncRNAs in 7 tissues (gastrocnemius muscle, longissimus dorsi, heart, liver, spleen, lung and kidney) between WT and ZBED6 KO pigs.

Additional file 4.

Supplementary Table S4. Tissue-specific and common expression of lncRNAs across 7 tissues (gastrocnemius muscle, longissimus dorsi, heart, liver, spleen, lung and kidney) in WT pigs.

Additional file 5.

Supplementary Table S5. Cis and trans target genes predicted to be regulated by tissue-specific and commonly expressed lncRNAs in 7 tissues (gastrocnemius muscle, longissimus dorsi, heart, liver, spleen, lung and kidney) in WT pigs.

Additional file 6.

Supplementary Table S6. KEGG of targeted genes of tissue-specific lncRNAs across 7 tissues in WT pigs.

Additional file 7.

Supplementary Table S7. List of lncRNA modules.

Additional file 8.

Supplementary Table S8. Cis- and trans-target genes predicted to be regulated by 30 commonly expressed lncRNAs in longissimus dorsi tissue-specific expression (tissue-specific), skeletal muscle module (magenta module) and DE-lncRNAs of longissimus dorsi (LD_DE-lncRNAs).

Additional file 9.

Supplementary Table S9. KEGG enrichment of 30 commonly expressed lncRNAs.

Additional file 10.

Supplementary Table S10. The detail of RNA-seq data.

Rights and permissions

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

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Wang, D., Pu, Y., Li, Y. et al. Comprehensive analysis of lncRNAs involved in skeletal muscle development in ZBED6-knockout Bama pigs. BMC Genomics 22, 593 (2021). https://doi.org/10.1186/s12864-021-07881-y

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-021-07881-y

Keywords

  • ZBED6 KO pig
  • lncRNAs
  • RNA-Seq
  • Skeletal muscle development