LncRNAs are regulated by chromatin states and affect the skeletal muscle cell differentiation

Abstract Objective This study aims to clarify the mechanisms underlying transcriptional regulation and regulatory roles of lncRNAs in skeletal muscle cell differentiation. Methods We analysed the expression patterns of lncRNAs via time‐course RNA‐seq. Then, we further combined the ATAC‐seq and ChIP‐seq to investigate the governing mechanisms of transcriptional regulation of differentially expressed (DE) lncRNAs. Weighted correlation network analysis and GO analysis were conducted to identify the transcription factor (TF)‐lncRNA pairs related to skeletal muscle cell differentiation. Results We identified 385 DE lncRNAs during C2C12 differentiation, the transcription of which is determined by chromatin states around their transcriptional start sites. The TF‐lncRNA correlation network showed substantially concordant changes in DE lncRNAs between C2C12 differentiation and satellite cell rapid growth stages. Moreover, the up‐regulated lncRNAs showed a significant decrease following the differentiation capacity of satellite cells, which gradually declines during skeletal muscle development. Notably, inhibition of the lncRNA Atcayos and Trp53cor1 led to the delayed differentiation of satellite cells. Those lncRNAs were significantly up‐regulated during the rapid growth stage of satellite cells (4‐6 weeks) and down‐regulated with reduced differentiation capacity (8‐12 weeks). It confirms that these lncRNAs are positively associated with myogenic differentiation of satellite cells during skeletal muscle development. Conclusions This study extends the understanding of mechanisms governing transcriptional regulation of lncRNAs and provides a foundation for exploring their functions in skeletal muscle cell differentiation.


| INTRODUC TI ON
The formation of myotubes is a necessary step in the complex and multi-stage process of muscle development, which is affected by many cytokines and self-proteins. 1 The vitality of satellite cells is essential for maintaining the stability of skeletal muscle tissue. 2 Moreover, both the capacity for differentiation and the number of satellite cells decrease with age during skeletal muscle development in mice. 3 However, the current understanding of gene regulation during the processes of skeletal muscle cell differentiation and functional decline in satellite cells remains incomplete, especially regarding the role of long non-coding RNAs (lncRNAs).
Transcription factors (TFs), such as the myogenic regulatory factors MyoD family, the MEF2 family (MEF2A-D) and others, have been shown to determine the characteristics of skeletal muscle cells. [4][5][6] The MyoD, MyoG and Mef2c are specifically involved in the regulation of myogenic differentiation, and their expression is significantly reduced and delayed in satellite cell senescence. 7,8 Although Heyl and Hey1 are considered potential effectors of the Notch pathway to inhibit myogenic differentiation, only constitutive expression of Hey1 blocked myogenesis. 9 Another study found that the weight and size of Heyl/1 double-knockout mice decreased, which indicates that Heyl has a positive effect on muscle development. 10 These previous findings demonstrate that the expression of TFs can be used as a molecular marker of cell differentiation state and function.
A growing body of evidence supports the finding of a close relationship between lncRNA function and skeletal myogenesis, and muscle diseases. 11,12 The Linc-MD1 RNA was reported to be involved in skeletal muscle differentiation through regulation of myogenic TFs. 13,14 Similarly, Linc-YY1 promotes myogenic differentiation via the interactions with YY1 and regulates satellite cell activation/proliferation by regulating Pax7 expression. 6 LncMyoD has been found to induce myogenic differentiation through disrupting the cell cycle. 15 The lncRNA Trp53cor1, also known as LincRNA-p21, has been shown to participate in repressing cell proliferation and smooth muscle cell apoptosis. 16,17 Collectively, these discoveries provide strong evidence of a contribution by lncRNAs to skeletal muscle development, although their transcriptional regulation and regulatory roles have not been well-studied at the whole genome level.
In this work, we systematically analysed expression patterns of ln-cRNAs using RNA-seq over a time course spanning from proliferation to differentiation of C2C12 myoblasts and mouse satellite cells across a range of ages to understand their regulatory functions. Moreover, the transcriptional regulation of lncRNAs was examined by combining ChIP-seq and ATAC-seq data. Our study provides new insights into the transcriptional regulation of lncRNAs and the mechanisms underlying their regulatory roles in skeletal muscle cell differentiation.

| Sample collection and cell culture
The C2C12 mouse myogenic cell line was acquired from the Cell Bank of Wuhan University and cultured in 1 × DMEM basic (Gibco-BRL) with 20% foetal bovine serum (Gibco 10099133) in 5% CO 2 at 37°C.
Skeletal muscle satellite cells were isolated from the hindlimb muscle of C57BL/6 mice at weeks 2, 4, 6, 8, 10 and 12, which were represented as W2, W4, W6, W8, W10 and W12, respectively. The detailed isolation process for satellite cells is described in our previous report. 3

| Analysis of differentially expressed lncRNAs
The mouse reference genome and GTF files of lncRNAs and protein-coding genes were downloaded from the GENCODE database version M11 (https://www.genco degen es.org/mouse/). The expression levels of protein-coding genes and lncRNAs were quantified by tophat v2.1.1 18 and the htseq-count script of htseq v0.6.0. 19 Differential expression analysis of C2C12 was performed on transcriptomic data collected between proliferation (cells cultured in growth medium, GM) and differentiation (cells cultured in differentiation medium, DM 12, 24, 48, 60 and 96 hours) by edger. 20 Differentially expressed (DE) genes were identified based on a fold change threshold value of (|log 2 FC| ≥ 1) and P value (P < .05).
Add 1 mL cold ATAC-buffer (10 mmol/L Tris-HCl, pH 7.4, 10 mmol/L NaCl, 3 mmol/L MgCl 2 and 0.1% Tween-20) was added to stop the lysis. The cells were further centrifuged at 500 RCF for 10 minutes at 4°C to collect the cell pellet after removing the supernatant. The following transposition reaction mix ( 21 Final library was electrophoresed in a 2% high-resolution agarose gel and 100-1000 bp fragments were cut out and sequenced by using Illumina HiSeq X Ten PE150 platform. The information of ATAC-seq data is shown in Table S1 and Figure S1.

| ATAC-seq and ChIP-seq data processing
Transcription factor (MyoD, MyoG, Cebpb, Usf1 and Max 23 ) and histone modification (H3K4me3 23 and H3K4me1 24 ) ChIP-seq data were downloaded from the previous studies. These data were obtained from C2C12 cells at proliferation and differentiation stages.
The ATAC-seq reads were mapped to the mouse M11 reference genome from GENCODE database using bowtie2 v2. 3 quantiles function of preprocesscore package v1.40.0 29 in r v3.5.1 was used to normalize the ATAC-seq signal. ChIP-seq data were subjected to a similar process as ATAC-seq reads. The motif analysis in this study was performed using the findMotifs.pl function of homer software. 30

| Correlation analysis
Weighted correlation network analysis (WGCNA) 31 was used to identify patterns of co-expression between significant DE genes and lncRNAs of C2C12 cell line RNA-seq data. The WGCNA package in the r v3.5.1 environment was applied to develop a weighted correlation network based on the normalized count matrix produced during differential expression analysis of lncRNAs using edger. 20 The Pearson correlation analysis was also used to identify significantly correlated DE lncRNA-gene pairs in C2C12 (P < .05). DAVID Bioinformatics Resources 6.8 (https://david.ncifc rf.gov/summa ry.jsp) was used for GO term functional annotation of significantly enriched genes.

| qPCR
Total RNA from C2C12 and satellite cells was extracted with TRIzol reagent (Invitrogen 15596026) according to the manufacturer's instructions. cDNA was obtained via reverse transcription of 1 μg of RNA using a High-Capacity cDNA Reverse Transcription kit (Thermo Fisher 4374967). THUNDERBIRD SYBR qPCR Mix (Toyobo) was used with a CFX384 real-time PCR system (Bio-Rad). All primer sequences in this study are listed in Table S2.  Table S3. . C, qPCR validation of RNA-seq differential expression analysis using four randomly selected DE lncRNAs in differentiating C2C12 myoblasts (n = 3). Tubulin was used as the internal control. '*' and '**' indicate P < .05 and P < .01, respectively. Fold change in relative expression was based on expression levels during the proliferation stage. D, Scatter plot of the log 2 FC correlation between RNA-seq (y-axis) and qPCR (x-axis) 4°C overnight with gentle shaking. Cells were further washed four times with 1 × PBS and then incubated with anti-mouse IgG (H+L), F (ab') 2 Fragment (Alexa Fluor™ 555 Conjugate; CST, 4409). A Nikon Eclipse TE2000-S microscope (Nikon) was used to observe the fluorescence.

| Western blot
Cells were lysed in RIPA buffer (Sigma R0278) with phosphatase and protease inhibitors on ice for 30 minutes. Total protein was electrophoresed in 10% SDS-PAGE gel and transferred onto a PVDF membrane (Millipore). The membranes were blocked in 5% skim milk for 2 hours at room temperature and then incubated overnight at

| Dynamic expression of lncRNAs during C2C12 myoblast differentiation
LncRNAs are involved in multiple biological processes in skeletal muscles. 11 Through the GENCODE project initiative, 9989 lncRNAs have been identified in the mouse genome. Among them, expression of 994 lncRNAs was detected in C2C12 proliferation and differentiation cells via the rRNA depletion strand-specific RNA-seq method used in our study. Here, 385 (38%) of DE lncRNAs and 3042 (25%) of DE protein-coding genes (PCGs) were identified via RNAseq over a time course spanning the C2C12 myoblast differentiation process ( Figure 1A). Visualization by heat map revealed dynamic changes in the expression of these DE lncRNAs during C2C12 myoblast differentiation ( Figure 1B). To verify the accuracy of our expression pattern analysis, we randomly selected four DE lncRNAs for qPCR validation, namely Gm10125, GM28653 and Malat1, which were up-regulated, and RP23-115O21.3, which was down-regulated ( Figure 1C). We found a high correlation between RNA-seq and qPCR (R = .89, P < .01; Figure 1D), confirming the reliability of our differential expression analysis. For the third category of DE lncRNAs, that is, associated with ATAC only, the difference between up-and down-regulated DE lncRNAs in the ATAC-seq signal was non-significant, although the distribution was similar to that of the H3K4me1-and H3K4me3associated categories ( Figure 2C). However, the TF binding results for this category were completely different in that the proportion of DE lncRNAs associated with MyoD or MyoG binding was substantially lower than the other two categories. In addition, recognition motifs of other TFs were significantly enriched in these ATAC-seq peak regions ( Figure 2G), and the expression levels of these TFs changed over the course of C2C12 myoblast differentiation ( Figure 2H).

| The correlations between DE lncRNAs and genes were revealed by WGCNA and Pearson's analysis
To explore the relationship between the expression of lncRNAs and protein-coding genes during myogenic differentiation, we  Figure 4D), which were also previously shown to provide key contributions to myogenic differentiation. [36][37][38] We thus developed a correlation network between these four TFs and 149 of their correlated DE lncRNAs to improve our understanding of the function of lncRNAs ( Figure 4E, Figure S3). Based on the changes in expression between proliferation and differentiation observed in this network, we divided the lncRNAs into an up-regulated group and a down-regulated group ( Figure 4F).

| LncRNAs are involved in the regulation of differentiation of satellite cell during skeletal muscle development
Satellite cells are essential components of skeletal muscle development and regeneration. To investigate the functional differences of satellite cells at different ages in mice, we followed the instruction from a previous study to isolate quiescent satellite cells from mouse skeletal muscle and adherent culture to obtain activated satellite cells. 3 Through expression analysis, we found that MyoG, MyoD and MyHC were all up-regulated in the satellite cells of 4-and 6-weekold mice compared with those in 2-week-old mice, indicating the rapid growth of skeletal muscle tissues ( Figure 5A). Thereafter, they were down-regulated, commensurate with age-related decreases in proliferation and the capacity for satellite cell differentiation. These results were consistent with our previous results generated from satellite cells of the same developmental stage which showed the capability of mouse muscle satellite cells to self-renew, as well as a significant decrease in their ability to differentiate after a rapid period of growth (4-6 weeks). 3 We further investigated whether or not the lncRNAs involved in our TF-lncRNA correlation network constructed with C2C12 differentiation data were involved in the regulation of the mouse satellite cell differentiation capacity. Our results showed that the up-regulation and down-regulation of the lncRNAs in the TF-lncRNA network were significantly concordant (P < .05) between the satellite cells in different rapid growth stages (4-6 weeks) and differentiated C2C12 via gene set enrichment analysis ( Figure 5B,C, Figure S4A).
Moreover, these lncRNA-correlated TFs exhibited the same changes as the lncRNAs ( Figure 5C,D, Figure S4). This finding thus indicated that these TFs had a regulatory function in the proliferation and differentiation of myogenic cells.
In the satellite cells from 8 to 12 week mice, the up-and down-regulation of TFs and lncRNAs identified in this correlation network were opposite to those in the satellite cells from 2-weekold mice ( Figure 5B,D, Figure S4B,C), and especially so for lncRNAs which were up-regulated in the rapid growth stages of satellite cells.
These results further suggested that these TFs and lncRNAs were also involved in the regulation of differentiation capacity decrease of satellite cells with age.  Figure 6A,B). Moreover, the expressions of these TFs and lncR-NAs, which peaked at week 6, were all positively correlated, in agreement with our data analysis, which also showed a significant  Figure 6A,B). Together, these results indicated that our expression analysis was reliable.

| Inhibition of Atcayos and Trp53cor reduced the differentiation of satellite cells
To further explore the contribution of lncRNAs in the capacity for differentiation among satellite cells, we chose the lncRNAs Atcayos and Trp53cor, which were up-regulated in rapid growth stages, for further loss of function assays. The expression of Atcayos and Trp53cor1 was successfully knocked down using RNAi ( Figure 7A,B). The Western blot results showed that myosin decreased when Atcayos or Trp53cor1 were inhibited ( Figure 7C,D), which was confirmed using immunofluorescence microscopy ( Figure 7E-H). These results showed that inhibition of the TFassociated lncRNAs Atcayos and Trp53cor1 led to the delayed differentiation of satellite cells, which also indicated that these two lncRNAs participated in positive regulation of the differentiation of satellite cells.

| D ISCUSS I ON
In the current study, we examined the differential expression patterns of lncRNAs during C2C12 myoblast differentiation to better Previous studies reported that the expression activation of ln-cRNA was associated with H3K4me3, 39 H3K4me1, 40 etc, histone modifications. In our study, a small fraction of DE lncRNAs (17.07%) related to myoblast differentiation did not associate with H3K4me3 and H3K4me1 but only associated with ATAC signals at their TSS regions. Previous study also reported that lncRNA expression is regulated by TFs in a similar way of protein-coding genes. 39,41 Thus, the regulation of TFs becomes particularly important for above only ATAC-associated DE lncRNAs. Moreover, the motifs of TFs Runx family (Runx1 and Runx2), ATF3, AP-1 family (Jund and Fosl1) and CTCF were significantly enriched in the ATAC-seq peak regions at TSS regions of the above DE lncRNAs. Previous studies showed that the TFs Runx1, 42 ATF3 43 and CTCF were related to muscle development. CTCF is essential in mediated chromatin loops, 44 delimiting enhancer-promoter interaction, 45 and is implicated in gene activation. 46 Thus, we speculated that these TFs play critical regulation roles in expression regulation of lncRNAs which were differentially Previous research has shown that Trp53cor1 functions in the regulation of vascular smooth muscle cell proliferation and apoptosis, 17 and also promotes migration of mesenchymal stem cells. 54 However, the lncRNA Atcayos remains functionally uncharacterized at present. We found, using siRNA knockdowns, that Trp53cor1 and Atcayos were both positively associated with the differentiation of satellite cells during rapid growth, thereby improving our current understanding of the regulatory function of lncRNAs in myogenic differentiation.
In conclusion, by combining ATAC-seq, ChIP-seq and RNA-seq data, we found that dynamic changes in the expression of lncRNAs between the developmental stages of myoblast proliferation and differentiation were closely related to chromatin states. In addition,

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

AUTH O R CO NTR I B UTI O N S
XQ and MH conceived and performed the experiments and explained the data. They also drafted, wrote and revised the manuscript and approved the version to be published. YX, YYX and YL analysed and helped explain the data. They also drafted and revised the manuscript and approved the version to be published. DW, YH, HZ, ZW and WZ assisted in the experiments and helped explain the data. They also drafted and revised the manuscript and approved the version to be published. SZ, XL and YZ designed the study and F I G U R E 6 qPCR relative expression analysis of transcription factors (TFs) and lncRNAs in C2C12 myoblasts and satellite cells. Three biological replicates were used for each analysis. Tubulin was used as the internal control. The '*' and '**' indicate P < .05 and P < .01, respectively. Scatter plots show correlations in expression between Heyl or Mef2c and lncRNAs in C2C12 myoblasts and satellite cells. A, Expression of TFs (Mef2c and Heyl) and lncRNAs (Atcayos, Trp53cor1, and GM10561) determined via qPCR in the C2C12 cell line. All comparisons were made with expression levels at the proliferation stage (GM). B, The results of qPCR and correlation of relative expression between TFs (Mef2c and Heyl) and lncRNAs (Atcayos, Trp53cor1, and GM10561) in satellite cells. Fold change is based on comparisons with expression in satellite cells from 2-week-old mice F I G U R E 7 Inhibition of Atcayos and Trp53cor1 in the differentiated satellite cells. A, B, qPCR results of Atcayos and Trp53cor1 expression in differentiated satellite cells at 24 h when Atcayos or Trp53cor1 were inhibited using siRNA, '**' indicate P < .01. C, D, Western blots of Atcayos and Trp53cor1 in siRNA knockdown differentiated satellite cells at 24 h. E, F, Immunofluorescence staining of myosin (red) in the differentiated satellite cells at 24 h when Atcayos or Trp53cor1 were inhibited using siRNA. Nuclei were stained with DAPI (blue) Scale bars: 100 μm. Magnification: 100×. G, H, Differentiation index (no. of nuclei in myosin + cell/total nuclei) for satellite cells in (E) and (F) explained the results. They also revised it critically for important intellectual content and approved the version to be published.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data are openly available in a public repository that issues data sets with DOIs: ChIP-seq data of H3K4me3, MyoD, MyoG, Usf1, Max and Cebpb that support the results of this study are openly available in NCBI at https://doi.org/10.1038/natur e13992. 23 ChIP-seq data of H3K4me1 that support the findings of this study are openly available in NCBI at https://doi.org/10.1016/j.molcel.2013.07.022. 24 Satellite cell RNA-seq data that support the findings of this study are openly available in NCBI at https://doi.