Characterization of Long Non-coding RNAs Modified by m6A RNA Methylation in Skeletal Myogenesis
Shu-Juan Xie1,2† 
Shuang Tao2† 
Li-Ting Diao2† Pan-Long Li2 Wei-Cai Chen2 Zhi-Gang Zhou3 Yan-Xia Hu2 Ya-Rui Hou2 Hang Lei2 Wan-Yi Xu2 Wen-Jie Chen2 
Yan-Wen Peng2* 
Qi Zhang1,2* 
Zhen-Dong Xiao2*- 1Vaccine Research Institute of Sun Yat-sen University, The Third Affiliated Hospital of Sun Yat-sen University, Guangzhou, China
- 2Biotherapy Center, The Third Affiliated Hospital of Sun Yat-sen University, Guangzhou, China
- 3Department of Orthopedics, First Affiliated Hospital, Jinan University, Guangzhou, China
Proper development of mammalian skeletal muscle relies on precise gene expression regulation. Our previous studies revealed that muscle development is regulated by both mRNA and long non-coding RNAs (lncRNAs). Accumulating evidence has demonstrated that N6-methyladenosine (m6A) plays important roles in various biological processes, making it essential to profile m6A modification on a transcriptome-wide scale in developing muscle. Patterns of m6A methylation in lncRNAs in developing muscle have not been uncovered. Here, we reveal differentially expressed lncRNAs and report temporal m6A methylation patterns in lncRNAs expressed in mouse myoblasts and myotubes by RNA-seq and methylated RNA immunoprecipitation (MeRIP) sequencing. Many lncRNAs exhibit temporal differential expression, and m6A-lncRNAs harbor the consensus m6A motif “DRACH” along lncRNA transcripts. Interestingly, we found that m6A methylation levels of lncRNAs are positively correlated with the transcript abundance of lncRNAs. Overexpression or knockdown of m6A methyltransferase METTL3 alters the expression levels of these lncRNAs. Furthermore, we highlight that the function of m6A genic lncRNAs might correlate to their nearby mRNAs. Our work reveals a fundamental expression reference of m6A-mediated epitranscriptomic modifications in lncRNAs that are temporally expressed in developing muscle.
Introduction
Skeletal muscle plays critical roles in the regulation of wider metabolism as well as driving locomotion (Cong et al., 2020). Myogenesis, the development of muscle, is a complex biological process regulated by multiple transcription factors and specific signaling pathways (Bryson-Richardson and Currie, 2008; Bentzinger et al., 2012). Our previous studies showed that non-coding RNAs, including miRNAs and lncRNAs (long non-coding RNAs), play essential roles in skeletal muscle development (Xie et al., 2018; Liu et al., 2021; Tan et al., 2021). LncRNAs are a class of non-coding RNAs greater than 200 nucleotides in length with limited or no protein-coding capacity that possess complex spatial structures and diverse functions (Chen and Carmichael, 2010; Xie et al., 2021a). Numerous studies have shown that lncRNAs play a significant role in biological functions, such as epigenetic modification, mRNA transcription, splicing, stability and translation (Lan et al., 2021). Functionally, lncRNAs can either act in cis by regulating the expression of neighboring genes or in trans by regulating the expression of distant genes (Ulitsky and Bartel, 2013). Increasing studies have shown that lncRNAs participate in myogenesis. For example, H19, one of the earliest known imprinted lncRNAs, is strongly repressed after birth in all mouse tissues, but it remains expressed in the skeletal muscle and heart in adults (Milligan et al., 2000), controlling reactivation of the imprinted gene network and alleviating muscular dystrophy (Borensztein et al., 2013; Martinet et al., 2016; Zhang Y. et al., 2020). LncRNA MALAT1 interacts with miRNAs or mRNAs to regulate skeletal muscle maintenance (Yong et al., 2020; Liu et al., 2021). Additionally, several lncRNAs have been reported to shape muscle (Ro et al., 2018; Sweta et al., 2019; Martone et al., 2020), including linc-MD1 (Cesana et al., 2011), lincYY1 (Lu et al., 2013; Zhou et al., 2015), lncRNA Dum (Wang et al., 2015), linc-RAM (Yu et al., 2017; Zhao et al., 2018), muscle-specific lncR-Irm (Sui et al., 2019), lncR-Myoparr (Hitachi et al., 2019), lnc-MyoD (Gong et al., 2015; Lim et al., 2020), and lncMGPF (Lv et al., 2020).
RNA chemical modifications in coding and non-coding RNAs can regulate gene expression without changing the sequence of the RNA molecules via a process referred as “epitranscriptomics” (Helm and Motorin, 2017; Roundtree et al., 2017; Fazi and Fatica, 2019). Greater than 150 RNA modifications have been identified as posttranscriptional regulatory marks in multiple RNA species, including mRNAs, tRNAs, rRNAs, small non-coding RNAs, and lncRNAs (Yang Y. et al., 2018). Among these modifications, N6-methyladenosine (m6A) is the most common modification in mammalian mRNAs and lncRNAs (Pan, 2013). The effectors in m6A pathways include “writers” and “erasers” that install and remove the methylation, respectively, and “readers” that recognize it (Fu et al., 2014; Meyer and Jaffrey, 2017; Shi et al., 2019). The m6A writer machinery is a methyltransferase complex composed of multiple subunits with a stable core complex formed between methyltransferase-like 3 (METTL3) and methyltransferase-like 14 (METTL14) (Ping et al., 2014). Benefitting from deep sequencing, m6A patterns have been demonstrated to occur in a cell type- and cell state-dependent manner, and the landscape of the m6A methylome has been identified in human and mouse tissues (Liu et al., 2020; Zhang H. et al., 2020).
Given the important function of m6A modification in gene expression, emerging evidence has revealed the critical role of m6A in skeletal muscle regulation (Li et al., 2021). Our recent study elaborated the dynamic m6A methylation of mRNA during skeletal muscle differentiation and revealed the role of METTL3/14-m6A-MNK2-ERK signaling axis in skeletal muscle differentiation and regeneration (Xie et al., 2021b). Besides, another work of us confirmed that the m6A key methyltransferase METTL3 is involved in the biogenesis of muscle-specific miRNAs (Diao et al., 2021a). This finding is consistent with a previous report indicating that METTL3 is sufficient to enhance miRNA maturation in a global and non-cell-type specific manner (Alarcón et al., 2015). Other studies focused on the modification of mRNA by m6A methylation, which is involved in muscle formation, maintaining muscle homeostasis, and musculoskeletal disorders (Zhang W. et al., 2020). Recent studies revealed that METTL3-mediated m6A methylation is essential for muscle stem cell self-renewal (Lin et al., 2020), muscle regeneration (Liang et al., 2021), and muscle stem cell/myoblast state transitions (Gheller et al., 2020). Interestingly, myogenic potential is maintained partly by the Mettl3-mediated stabilization of processed MyoD mRNA through m6A modification of the 5′ untranslated regions (UTR) during proliferative phases (Kudou et al., 2017), and depletion of m6A “eraser” FTO in myoblasts leads to impaired skeletal muscle development (Wang et al., 2017). These data suggested that m6A modification could mediate muscle progenitor cell proliferation and differentiation. It has been demonstrated that m6A lncRNA modification plays roles in different biological processes (Patil et al., 2016; Yang D. et al., 2018; Fazi and Fatica, 2019; Ma et al., 2019; He et al., 2020; Lan et al., 2021). However, little is known about the m6A methylation status of lncRNAs involved in developing skeletal muscle.
In this study, we characterized lncRNAs in mouse myoblasts and differentiated myotubes using RNA sequencing (RNA-seq) and further uncovered abundant m6A sites and specific m6A patterns in these lncRNAs using methylated RNA immunoprecipitation sequencing (MeRIP-seq). Our results reveal the temporal expression profile and m6A methylation status of lncRNAs during skeletal myogenesis. We found that the m6A methylation levels of lncRNAs were positively correlated with their transcriptional abundance. Our data will provide a fundamental reference for further study on the function of lncRNAs.
Materials and Methods
Cell Culture
The C2C12 mouse myoblast cell line and HEK-293T cells were purchased from the Cellular Library of the National Collection of Authenticated Cell Cultures (Shanghai, China). The cells were cultured in growth medium (GM)-Dulbecco’s modified Eagle’s medium (DMEM, Gibco) with 10% fetal bovine serum (FBS, Gibco), 100 U/ml penicillin, and 100 μg/ml streptomycin (1 × penicillin–streptomycin) at 37°C in a humidified chamber supplemented with 5% CO2. When C2C12 cells reached about 90% confluency, the GM was replaced with differentiation medium (DM)-DMEM containing 2% horse serum (HyClone).
Stable Cell Generation
For METTL3 overexpression, cDNA of mouse METTL3 was cloned into the pKD-CMV-MCS-EF1-PURO (pKD) vector by Gibson Assembly, and pKD-GFP was used as a negative control. For METTL3 knockdown, the gRNAs downstream of the transcription start sites were used to guide the fusion of inactive Cas9 (dCas9) to the Krüppel-associated box (KRAB) repressor. To generate stable cells, lentiviruses were produced in HEK-293T cells by transfecting vectors together with psPAX.2 and pMD2.G. Lentiviruses were collected and filtered 48 h after transfection and then used to infect target cell lines. Stable cells were selected using the antibiotic puromycin.
All sequences of clone primers used in this study are listed in Supplementary Table 1.
Gene Knockdown
To knock down lncRNAs, custom designed siRNAs targeting selected lncRNAs and control siRNAs were synthesized by Shanghai GenePharma Co., Ltd. C2C12 cells were seeded in 12-well plates and transfected with siRNAs using Lipofectamine 2000 (Invitrogen) after the cells reached 30–40% confluency, according to the manufacturer’s instructions. The transfected cells were maintained in growth medium for two days, and then cells were harvested for analysis.
All siRNA sequences used in this study are listed in Supplementary Table 2.
m6A Dot Blot Assay
RNA dot blotting was performed as previously described with modifications (Chen et al., 2015). Cells were harvested carefully and purified using a Dynabeads mRNA direct kit (Invitrogen, 61012). To avoid RNA degradation, RNase-free tubes and RNase-free water were used. RNA samples were quantified, diluted and incubated at 95°C in a heat block for 3 min to disrupt secondary structures. The tubes were chilled on ice immediately after denaturation for 2 min. The RNA samples were dropped onto the membrane (Amersham Hybond-N +, GE) and allowed to air dry for 5 min. Then, RNA was crosslinked to the membrane with UV light (2 autocrosslink, 150 mJ/cm2 UV Stratalinker, STRATAGENE). The membrane was washed in TBST (1TBS, 0.1% Tween-20), dyed in methylene blue (Sigma–Aldrich) as a quantitative control, and incubated in blocking buffer containing 5% non-fat dry milk in TBST for 2 h at room temperature. Then, the membrane was incubated with m6A antibody (1:1000, Synaptic Systems, 202-003) at 4°C overnight. The membrane was washed 3 times for 10 min each in TBST and then incubated with HRP-linked secondary anti-rabbit IgG antibody (1:5000, CST, 7074) for 1 h at room temperature. Signals were detected with WesternBright ECL HRP substrate (Advansta). The dots were quantified using ImageJ.
Western Blot
Cells were lysed in ice-cold enhanced RIPA lysis buffer (Shanghai Wansheng Haotian Biological Technology) containing phosphatase inhibitor and protease inhibitor cocktail (Roche). Equivalent total protein extracts were separated by SDS–PAGE and transferred to nitrocellulose membranes (Merck Millipore). The membranes were blocked with 5% non-fat dry milk in TBST for 1 h at room temperature. The following antibodies were used in this study: anti-METTL3 (Proteintech, 15073-1-AP), anti-METTL14 (R&D, HPA038002), anti-MHC (R&D, MAB4470), and anti-GAPDH (CST, 2118). Immunoreactivities were determined using WesternBright ECL HRP substrate (Advansta).
RNA Isolation and Quantitative RT-PCR Assay
Total RNA was extracted from cells with TRIzol reagent (Invitrogen, 15596018). First-strand cDNA for PCR analyses was synthesized with HiScript III RT SuperMix for qPCR (+ gDNA wiper) (Vazyme, R323-01), and quantitative real-time PCR was performed using ChamQ Universal SYBR qPCR Master Mix (Vazyme, Q711-02). The GAPDH gene served as an endogenous control. The qRT-PCR results were analyzed and presented as relative RNA levels of the CT (cycle threshold) values, which were then converted as fold change. The results are presented as the means ± SD.
All primers for qPCR are listed in Supplementary Table 1.
RNA Library Construction and Sequencing
Total RNA was isolated and purified using TRIzol reagent (Invitrogen, 15596018) following the manufacturer’s procedure. The RNA amount and purity were quantified by a NanoDrop ND-1000 (NanoDrop). RNA integrity was assessed using the Agilent 2100 system with RIN number >7.0 and confirmed by electrophoresis with a denaturing agarose gel.
RNA-seq was performed by BGI Co., Ltd. Briefly, total RNA was used to purify the poly-A containing RNAs using poly-T oligo-attached magnetic beads. Following the purification, the remainder of the RNA was fragmented into small pieces using divalent cations under high temperature. Then, the cleaved RNA fragments were reverse transcribed to create the cDNA library according to the mRNA-Seq sample preparation kit protocol (Illumina). The average insert size for the paired-end libraries was 300 bp (±50 bp). Then, paired-end sequencing was performed on an Illumina NovaSeqTM 6000 (BGI Co., Ltd, Shenzhen, China) following the vendor’s recommended protocol.
MeRIP Sequencing
RNA extraction and quality control were performed as previously noted. Poly(A) RNA was purified from 50 μg of total RNA using Dynabeads Oligo (dT) 25-61005 (Thermo Fisher) and fragmented into 100-nucleotide-long oligonucleotides using Magnesium RNA Fragmentation Module (NEB). The cleaved RNA fragments were immunoprecipitated using an anti-m6A affinity purified antibody (Synaptic Systems, 202003). Then, the IP RNA was reverse transcribed to cDNA by SuperScriptTM II Reverse Transcriptase (Invitrogen, 1896649) followed by U-labeled second-stranded DNA synthesis. Then, after ligation with the adapter to the A-tailed fragmented DNA, the ligated products were amplified with PCR, and 2 × 150 bp paired-end sequencing (PE150) was performed on an Illumina NovaSeqTM 6000 (LC-Bio Technology CO., Ltd., Hangzhou, China) following the vendor’s recommended protocol.
m6A-IP-qPCR
m6A methylated RNA immunoprecipitation was performed as previously described (Diao et al., 2021a). Briefly, 200 μg of total RNA was incubated with 8 μg of m6A-specific antibody (Synaptic Systems, 202003) for 2 h at 4°C with gentle rotation. Then, 50 μl protein A/G magnetic beads were added and incubated for 2 h at 4°C with gentle rotation. Beads were pelleted at 2,500 rpm for 30 s. Then, the supernatant was removed, and the beads are resuspended in 500 μL RIP buffer. The process is repeated for a total of three RIP washes followed by one wash in PBS. Beads were incubated with DNase I for 30 min at 37°C and were then digested by Proteinase K for 2 h at 37°C with rotation. A MicroElute RNA Clean Up Kit (Omega) was used for RNA purification. Purified RNA was reverse transcribed and quantified by real-time RT-PCR.
RNA-Seq Analysis
The sequencing adapters were removed from raw fastq data using cutadapt (Kechin et al., 2017) software, and then clean reads were mapped to the mouse reference genome (mm10) using HISAT2 (Kim et al., 2019). According to reference gene annotation (Ensembl release 102), the raw counts of each gene were calculated using featureCounts (Liao et al., 2014). Raw counts were further normalized as reads per kilobase of genome per million mapped reads (RPKM) using the fpkm function in the DESeq2 (Love et al., 2014) package. Differentially expressed genes were identified using DESeq2 (Love et al., 2014) with adjusted P ≤ 0.05. LncRNAs were selected for downstream analysis based on gene type.
MeRIP-Seq Analysis
The sequencing adapters were removed from raw fastq data using cutadapt (Kechin et al., 2017) software, and then clean reads were mapped to the mouse reference genome (mm10) using HISAT2 (Kim et al., 2019). Unique mapped reads were selected using samtools (Li et al., 2009) with a mapping quality greater than 30. To visualize read coverage using IGV (Robinson et al., 2011), the bigwig format of mapped reads was generated using deeptools (Ramírez et al., 2014). Mapped reads of IP and input libraries were fed into the exomePeak (Meng et al., 2013) package for calling peaks and identifying distinct peaks, and statistical significance was defined as FDR ≤ 0.05. Called peaks were annotated by intersecting with gene architecture using bedtools (Quinlan and Hall, 2010) and custom Python script. Peaks located at lncRNAs were selected based on gene type for downstream analysis.
Motif Identification and Peak Distribution Among lncRNA Bodies
Sequence motifs enriched in m6A peaks were identified by HOMER (Heinz et al., 2010) with “−len 5,6,7,8 −rna” and other default settings as parameters. An m6A metagene plot was plotted using the Guitar (Cui et al., 2016) package. To investigate the peak distribution of lncRNA exon elements, an in-house Python script was developed to calculate the numbers of peaks located at the first exon, internal exons and last exon. Pie charts were plotted using ggplot2 package.
Relationship Analysis of lncRNA m6A Level and Expression Abundance
Significant differential m6A peaks were integrated with corresponding differential expression data regardless of significance. Based on change orientation, the lncRNAs with significant differential m6A peaks were classified into four groups: hypermethylated and upregulated, hypermethylated and downregulated, hypomethylated and upregulated, and hypomethylated and downregulated. A four-quadrant diagram was generated using the ggplot2 package. To estimate the relationship of lncRNA m6A level and expression abundance, Pearson correlation analysis was performed using R software (version 4.0.3).
Neighbor Gene Analysis of Differentially Expressed lncRNAs
Ten neighboring mRNAs (five upstream and five downstream) of differentially expressed lncRNAs were obtained using an in-house Python script, and mRNAs with significantly different expression (adjusted P ≤ 0.05) were fed into the clusterProfiler (Yu et al., 2012) package for GO and KEGG analysis. Statistical significance was defined as p ≤ 0.05. The top 10 most significant terms are shown using dot plots. For genes associated with muscle development, their neighboring lncRNAs were also required to have significant m6A changes. Finally, a set of lncRNAs with expression and m6A alterations related to muscle development were identified.
Visualization Analysis
A heatmap was drawn using the pheatmap package, and other charts that were not specified were drawn using the ggplot2 package.
Results
Dynamic Profile of lncRNAs in Myoblasts and Differentiated Myotubes
We used mouse C2C12 myoblast cells to mimic skeletal muscle differentiation. C2C12 cells constantly proliferate in the presence of serum and begin to differentiate in the absence of serum. After differentiation for 4 days, obvious morphological changes were observed when myocytes fused to form multinucleated myotubes (Supplementary Figure 1A), as described in our previous reports (Xie et al., 2018; Tan et al., 2021). Next, RNA dot blotting was performed to investigate the dynamics of m6A RNA modification during myogenesis, and decreased global m6A levels were observed in myotubes that had been differentiated for 4 days (D4) compared to myoblasts (cells on growth medium, GM) (Figure 1A). To test whether this change was due to altered expression of m6A methyltransferases or demethylases, we profiled the core components of m6A methyltransferases METTL3, METTL14, and m6A demethylases FTO and ALKBH5 during C2C12 differentiation. Consistent with the decline in m6A levels, METTL3 and METTL14 protein expression was decreased on D4 and negatively correlated with the myogenic marker MHC (Supplementary Figure 1B). However, the demethylases FTO and ALKBH5 had opposite changes during C2C12 differentiation (Supplementary Figure 1C), and the upregulation of FTO is consistent with previous report (Wang et al., 2017). These data revealed that m6A and its core methyltransferase decreased during C2C12 differentiation. Such a dynamic change in m6A may contribute to the regulation of muscle genes in myoblasts and myotubes.
Figure 1. Dynamic profile of lncRNAs in myoblasts and differentiated myotube. (A) RNA dot blot assay of m6A methylation levels of C2C12 myoblasts in GM and D4. (B) MA plot shows the relationship of expression abundance and fold changes of lncRNAs in myoblasts (GM) and myotube (D4). log2 (MeanRPKM) represents gene expression values, log2Fold Change represents the fold change of lncRNAs at D4 compared to GM. Red dots represent 582 significantly up-regulated lncRNAs at D4 in relation to GM, adjusted P ≤ 0.05; Green dots, represent 97 significantly down-regulated lncRNAs at D4 in relation to GM, adjusted P ≤ 0.05; Yellow dots represent lncRNA without significantly differential expression, adjusted P > 0.05. Top differentially expressed lncRNAs were marked in purple. (C) The volcano plot shows significantly differentially expressed lncRNAs at GM and D4. Red dots represent 582 significantly up-regulated lncRNAs at D4, adjusted P ≤ 0.05; Green dots, represent 97 significantly down-regulated lncRNAs at D4, adjusted P ≤ 0.05; Yellow dots represent lncRNA without significantly differential expression, adjusted P > 0.05. Top differentially expressed lncRNAs were marked in purple. (D) The Venn diagram shows 97 and 582 significantly differentially expressed lncRNAs in GM and D4, and 4904 lncRNAs without significant difference between GM and D4. (E) The heatmap shows significantly differentially expressed lncRNAs in GM and D4. Yellow color represents up-regulation in D4, while green color represents down-regulation. Row represents genes, column represents samples, and each cell represents expression value. (F) Pie chart show the gene types of lncRNAs up-regulated in GM. (G) Pie chart show the gene types of lncRNAs up-regulated in D4. (H) Quantitative real-time reverse transcription PCR (qRT-PCR) validated top differently expressed lncRNAs in developing muscle cells. RPKM: Reads Per Kilobase per Million mapped reads. adjusted P: adjusted p-value; Data are presented as Mean ± SD; p value: **p < 0.01, ***p < 0.001, ****p < 0.0001.
In addition to RNA modification, non-coding RNAs, especially long non-coding RNAs (lncRNAs), could represent robust gene expression regulators. To profile the dynamic changes in lncRNAs during C2C12 differentiation, RNA sequencing (RNA-seq) was performed to analyze the expression of lncRNAs in GM and D4. Using bioinformatic analysis, we detected 5,583 lncRNAs expressed in at least one sample. Of these, 679 differentially expressed lncRNAs were identified (adjusted P value ≤0.05, Supplementary Table 3), We found that these differentially expressed genes generally tended to exhibit increased expression (Figure 1B). Moreover, the number of upregulated genes was far greater than that of downregulated genes (Figures 1B,C,E). We identified 97 and 582 lncRNAs significantly differentially expressed in GM and D4, respectively (Figure 1D). The top differentially expressed lncRNAs are marked in purple and are listed in Table 1. Among these lncRNAs, the majority were lincRNAs with 36.08% in GM and 34.19% in D4. The second rank was antisense RNA in both GM and D4 samples (Figures 1F,G). To validate the RNA-seq results, we performed quantitative real-time reverse transcription PCR (qRT-PCR) for the top nine upregulated and top five downregulated lncRNAs. All tested lncRNAs showed significantly differential expression, which was consistent with the RNA-sequencing results (Figure 1H). Of note, the top highly expressed lncRNAs, including Gm14635, Gm28653, 2310043L19Rik and 5430431A17Rik, were also highly expressed in the limb. In particular, Gm28653, which is also known as linc-MD1, is a muscle-specific lncRNA that controls muscle differentiation (Cesana et al., 2011). In addition, lncRNA 1700025l06Rik has the same transcript locus as lincMyod, which is directly activated by MyoD and regulates skeletal muscle differentiation by blocking IMP2-mediated mRNA translation (Gong et al., 2015).
Table 1. Top differentially expressed lncRNAs in myoblasts and differentiated myotube.
Meanwhile, we analyzed the mRNAs expression pattern of RNA-seq data, and identified 8,817 out of 16,190 differentially expressed mRNAs in D4 versus GM (adjusted P value ≤0.05, Supplementary Table 4). The number of upregulated coding genes was a little less than that of downregulated genes (Supplementary Figures 2A,B,D), of which 4483 and 4334 mRNAs significantly differentially expressed in GM and D4, respectively, and 7737 mRNA showed no significantly different expression in both group (Supplementary Figure 2C). The distinct difference in the number of genes with significant change between mRNAs and lncRNAs may due to the fact that mRNAs account for the majority of cellular RNA contents. Taken together, our RNA-seq analysis revealed the temporal expression of lncRNAs and mRNAs during myoblast differentiation.
Features of lncRNA m6A Methylation in Undifferentiated and Differentiated Muscle
Given the great importance of m6A methylation and non-coding RNAs in skeletal muscle development (Martone et al., 2020) as well as the altered m6A modification and lncRNA expression profiles in myoblast differentiation, we conducted methylated RNA immunoprecipitation sequencing (MeRIP-seq) in GM and D4 samples to uncover the m6A methylation landscape of lncRNAs. We identified greater than 20 thousand unique peaks from each MeRIP-seq sequencing library (FDR ≤ 0.05, Supplementary Tables 5, 6). Among them, we found 1383 lncRNA m6A peaks in the GM sample and 1848 lncRNA m6A peaks in the D4 sample (FDR ≤ 0.05, Supplementary Tables 7, 8). And we also found 5122 significantly differently mRNA peaks in D4 versus GM, of which 2692 mRNA m6A peaks in the GM sample and 1442 mRNA m6A peaks in the D4 sample (FDR ≤ 0.05, Supplementary Table 9 and Supplementary Figure 3). To identify whether m6A peaks share common sequence elements, we analyzed m6A binding motifs using Homer. We detected a significantly enriched consensus motif DRACU within lncRNAs from myoblast GM and myotube D4 (Figure 2A). The motif matched the well-validated consensus m6A motif DRACH (where D = A, G or U; R = A or G; H = A, C or U) and was similar to a previous report in lincRNAs (Meyer and Jaffrey, 2017). We explored the distribution of peaks along the lncRNA gene body and found that its density gradually decreased from the transcription initiation site to the transcriptional termination site (Figure 2B), which differs from the coding genes (Liu et al., 2015).
Figure 2. Features of lncRNA m6A methylation in undifferentiated and differentiated muscle. (A) The enriched consistent motif of m6A peaks in lncRNAs in GM and D4. (B) Metagene profiles of enrichment of all m6A peaks across lncRNAs transcriptome. (C) The top: pie charts represent the proportion of m6A peaks in the three regions of lncRNAs at GM. The bottom: histogram represents the relative enrichment of m6A peaks in the three regions of lncRNAs at GM. (D) The top: pie charts represent the proportion of m6A peaks in the three regions of lncRNAs at D4. The bottom: histogram represents the relative enrichment of m6A peaks in the three regions of lncRNAs at D4. (E) The frequence of m6A Peak Numbers in lncRNAs in GM and D4. (F) Bar plot shows the Numbers of m6A methylated lncRNAs in GM and D4. Blue represents 74 hyper-methylated lncRNAs in GM, while yellow represents 83 hyper-methylated lncRNAs in D4.
We further analyzed the peak distribution of lncRNA exons. We found that m6A peaks were preferentially enriched in the last exon of lncRNAs expressed in GM and D4 samples. In total, 63.69% and 67.25% m6A peaks were identified in the last exon of lncRNAs expressed in GM and D4, respectively. We then analyzed the peak enrichment in each lncRNA. Interestingly, we found that m6A peaks are preferentially enriched in the first exon and internal exon, but not the last exons (Figures 2C,D). We then analyzed the numbers of m6A peaks within lncRNAs and identified a median value of 1.0 m6A peaks per lncRNA (Figure 2E), and no differences were noted between the GM and D4 data. Furthermore, we performed integrating analysis by coupling MeRIP-seq and RNA-seq data. We mapped m6A peaks to differentially expressed lncRNAs and found that 74 and 83 lncRNAs were significantly hypermethylated in myoblasts (GM) and myotubes (D4), respectively (Figure 2F). These data provide a fundamental reference for the m6A epitranscriptome for further study.
Differentially m6A-Modified lncRNAs in Undifferentiated and Differentiated Muscle
To explore the putative function of m6A on lncRNAs, we investigated the correlation between lncRNA transcript abundance and m6A methylation. In total, 123 significantly differentially m6A-methylated lncRNAs were expressed during myogenesis, as shown by RNA-seq. Among these lncRNAs, we identified 34 hypermodified lncRNAs (32 hyper-upregulated and 2 hyper-downregulated) and 14 hypomodified lncRNAs (6 hypo-upregulated and 8 hypo-downregulated) according to the criteria of adjusted P ≤ 0.05 and FDR ≤ 0.05 (Figure 3A). These results indicate a temporal difference in m6A methylation in differentially expressed lncRNAs. To verify the relationship between m6A level changes of lncRNAs and expression changes of lncRNAs, we randomly selected 15 significantly differentially expressed lncRNAs (Table 2) and performed qRT-PCR analysis. The results showed that 10 hyper-upregulated lncRNAs increased and 5 hypo-downregulated lncRNAs decreased in D4 (Figure 3B). These findings were consistent with sequencing data. It is worth mentioning that lncRNA Brip1os, which is significantly downregulated in D4 samples, as shown in Figure 1H, was accompanied by a decline in m6A modification.
Figure 3. Differentially m6A modified lncRNAs in undifferentiated and differentiated muscle. (A) Distribution of genes with a significant change in both the m6A methylation and RNA expression levels before (GM) and after differentiation (D4), different colors were used to identify representative genes. And 15 m6A methylated significantly differently expressed lncRNAs were marked. (B) qRT-PCR validated the 15 m6A methylated significantly differently expressed lncRNAs in developing muscle cells. (C) Real-time PCR detection of the expression of the 15 m6A methylated significantly differently expressed lncRNAs in immunoprecipitated RNAs. IgG Immunoprecipitation was used as negative control. Quantitative data was represented as Mean ± SD; p value: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001. ns, no significant difference. (D) Integrative Genomics Viewer (IGV) plots show the m6A peaks of lncRNA Xist, Ptgs2os2 and Brip1os were highly enriched in m6A-RIP data. (E) Real-time PCR detection of the expression of the 5 hyper-upregulated and 5 hypo-downregulated lncRNAs in immunoprecipitated RNAs of GM and D4. Data were normalized by IgG Immunoprecipitation. Quantitative data was represented as Mean ± SD; p value: *p < 0.05, ** p < 0.01, ***p < 0.001, ****p < 0.0001. ns, no significant difference.
Table 2. Randomly selected 15 altered m6A peaks in myoblasts and differentiated myotube.
To verify the significantly differentially m6A-modified lncRNAs, we used an antibody against m6A and performed RNA immunoprecipitation followed by real-time PCR (m6A-IP-qPCR). As shown in Figure 3C, compared to that in the IgG control, most of the lncRNAs in Figure 3B were significantly enriched in the m6A group, indicating that these transcripts were m6A enriched. For example, the enrichment of Xist, Ptgs2os2 and Brip1os was elevated up to hundreds of thousands of fold in the m6A group, which is consistent with the MeRIP-seq data that revealed clear m6A peaks around their RNAs (Figure 3D). Furthermore, we verified the m6A peaks enrichment of 5 hyper-upregulated and 5 hypo-downregulated lncRNAs in GM and D4. By normalization of each group of Normal IgG, the m6A-IP-qPCR data was consistent with Figure 3A (Figure 3E). In summary, these results demonstrated that m6A methylation in lncRNAs is involved in myogenesis.
m6A Methylation Levels Were Positively Correlated With the Abundance of lncRNAs
More recently, m6A modification of mRNA was established to influence RNA stability dynamics and translation efficiency, and rapidly accumulating evidence shows significant crosstalk between lncRNA methylation and m6A-mediated epigenetic mechanisms (Kan et al., 2021). We then examined the correlation of lncRNA expression abundance with m6A methylation levels. For lncRNAs with significant expression abundance changes, their m6A levels were positively correlated with their expression levels (R = 0.6, P = 6.8e-6) (Figure 4A). However, for lncRNAs without significant transcript abundance changes, no significant correlation was noted between their m6A levels and expression levels (R = 0.21, p = 0.071) (Figure 4B). Interestingly, in both GM and D4 samples, lncRNAs with m6A methylation showed higher expression levels than those without m6A methylation (Figures 4C,D). Our analysis reveals that m6A methylation levels exhibit a positive correlation with the expression levels of m6A-modified lncRNAs and highlights the importance of m6A methylation in myogenesis-related lncRNAs.
Figure 4. m6A methylation levels were positively correlated with the abundance of lncRNAs. (A) Scatter plot shows the positive correlation between m6A levels (significant changes) and expression values of lncRNAs with significantly differential expression between GM and D4, adjusted P ≤ 0.05. (B) Scatter plot shows no correlation between m6A levels (no significant changes) and expression values of lncRNAs without significantly differential expression between GM and D4. (C) Cumulative frequency of log2FC for lncRNAs containing m6A or without m6A methylation in GM. Kolmogorov-Smirnov test was used to estimated inter-group difference. (D) Cumulative frequency of log2FC for lncRNAs containing m6A or without m6A methylation in D4. Kolmogorov-Smirnov test was used to estimated inter-group difference.
Myogenesis-Related lncRNAs Are Regulated by the m6A Methyltransferase METTL3
To further investigate whether altered m6A modification levels could affect lncRNA expression, we overexpressed METTL3 in C2C12 cells (oe-M3), and GFP-overexpressing cells were used as a negative control (GFP). The real-time PCR and Western blot results validated that METTL3 was successfully overexpressed (Figures 5A,B). Then, we assayed the expression levels of six selected lncRNAs in GFP- and METTL3-overexpressing cells. As shown in Figure 5C, all 6 lncRNAs, including Brip1os, were significantly upregulated when METTL3 was overexpressed.
Figure 5. Myogenesis associated lncRNAs are regulated by m6A methyltransferase METTL3. (A) qRT-PCR shows the RNA expression level of METTL3 in METTL3-overexpressing C2C12 cells. GFP-overexpressing C2C12 cells as negative control. (B) Western blot detected the protein expression levels of METTL3 in METTL3-overexpressing C2C12 cells. (C) qRT-PCR shows the expression of myogenesis associated lncRNAs in METTL3-overexpressing C2C12 cells. (D) qRT-PCR shows the RNA expression level of METTL3 in METTL3 knockdown C2C12 stable cell lines. A nonsense sequence constructed to dCas9 repressor as negative control. (E) Western blot detected the protein expression levels of METTL3 in METTL3 knockdown C2C12 stable cell lines. (F) qRT-PCR shows the expression of myogenesis associated lncRNAs when METTL3 was knockdown. Data are presented as Mean ± SD; p value: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.
To validate the METTL3 overexpression results, we designed two gRNAs downstream of transcription start sites and used them to guide the fusion of inactive Cas9 (dCas9) to the Krüppel-associated box (KRAB) repressor to inhibit the transcription of METTL3 in C2C12 cells. METTL3 mRNA and protein levels were greatly reduced compared to those with control gRNA (Figures 5D,E). Then, we assayed the 6 lncRNA expression levels, and qPCR results showed that the expression levels of these lncRNAs decreased when METTL3 was knocked down (Figure 5F). Taken together, our results revealed that lncRNA expression is positively correlated with m6A modification levels during myogenesis, and it might be a universal regulation way that m6A modification levels affect lncRNAs abundance.
m6A-Methylated lncRNAs Regulate Nearby mRNAs and Contribute to Muscle Tissue Development
Previous studies have reported that lncRNAs function in various physiological and pathological processes by regulating their adjacent mRNAs, either positively or negatively (Engreitz et al., 2016). Thus, we analyzed the significantly differentially expressed lncRNAs (FDR ≤ 0.05) during myogenesis as well as their nearest 10 mRNAs (upstream and downstream 5, respectively, adjusted P ≤ 0.05). GO analysis of biological processes for these mRNAs showed that 94 mRNAs were related to muscle tissue development (Figure 6A, adjusted P ≤ 0.05). KEGG pathway enrichment analysis showed that the cell cycle and MAPK signaling pathways were significantly enriched (Figure 6B, adjusted P ≤ 0.05), and such results are consistent with our previous studies showing that the JNK/MAPK and P38/MAPK signaling pathways play essential roles in myogenesis (Xie et al., 2018; Liu et al., 2021).
Figure 6. Functional relevance between m6A methylated lncRNAs and their adjacent mRNAs. (A) The top ten GO terms of the adjacent mRNAs (adjusted P ≤ 0.05) that related to differentially expressed lncRNAs in muscle cells. (B) The top twelve KEGG pathways of the adjacent mRNAs (adjusted P ≤ 0.05) that related to differentially expressed lncRNAs in muscle cells. (C) Table shows seven pairs of significantly differently expressed and methylated lncRNAs and their adjacent mRNAs in muscle cells. All have a significant threshold of FDR-adjusted p value ≤0.05. FDR, False Discovery Rate. (D) qRT-PCR shows the adjacent mRNA expression in GM and D4. (E) The effects of si-lncRNAs on the RNA expression levels of the corresponding lncRNA and mRNA. Data are presented as Mean ± SD; p value: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.
To confirm that m6A-methylated lncRNAs could regulate their adjacent mRNAs, we performed a synthetic analysis of 94 muscle tissue development-related mRNAs and nearby m6A-methylated lncRNAs. Briefly, given that nearby lncRNAs of these 94 mRNAs have significantly different expression, we further estimated m6A difference of these lncRNAs between GM and D4. As a result, seven lncRNAs have a significant m6A difference, and these lncRNAs also correspond to seven mRNAs (Figure 6C). Furthermore, mRNA expression of these paired adjacent mRNAs was tested, and qPCR results showed that Pi16, Cdon and Col14a1 were upregulated in D4 samples, consistent with their adjacent lncRNAs Gm41556, 4930581F22Rik, and Has2os, respectively. In contrast, Hdac4, Usp2 and Tbx2 were upregulated, which is an opposite effect compared with that noted in their adjacent lncRNAs (Figure 6D). These data indicated two opposite regulatory mechanisms between m6A-methylated lncRNAs and their nearby mRNAs, including positive and retrograde regulation.
Next, we further confirmed the regulation between lncRNAs and their adjacent mRNAs by knocking down corresponding lncRNAs using siRNAs. As shown in Figure 6E, when lncRNAs Gm41556, 4930581F22Rik, and Has2os were knocked down, their adjacent mRNAs Pi16, Cdon and Col14a1 were downregulated correspondingly, verifying the positive regulation between these lncRNAs and their nearby mRNAs. In contrast, when lncRNA Brip1os was knocked down, its adjacent mRNA Tbx2 was significantly upregulated, suggesting negative regulation. Taken together, our results showed that knockdown of m6A-methylated lncRNAs impacted the expression of their adjacent mRNAs and suggested that the functional relevance of m6A-methylated lncRNAs by regulating their adjacent mRNAs.
The METTL3/m6A/Brip1os/Tbx2 Axis in Muscle Development
Given that the lncRNA Brip1os exhibits markedly decreased expression and m6A modification levels during myogenesis and a perfectly negative correlation is noted between Brip1os and its nearby gene Tbx2, we further clarified their relationship. Brip1os and Tbx2 are both located at chr11qC, and these two transcription units have the same orientation. Brip1os is greater than 10 kb downstream of Tbx2 (359,435 bp) in the genome (Figure 7A). We further examined the expression of Tbx2 in muscle development by using single-cell RNA-seq data from Tabula Muris1 and analyzed the expression levels in skeletal muscle satellite stem cell and smooth muscle cell groups, which could be considered generally representative of undifferentiated myoblasts and differentiated myotubes, respectively. Surprisingly, although there were 439 samples in the skeletal muscle satellite stem cell group and only 42 samples in the smooth muscle cell group, Tbx2 expression levels in the smooth muscle cell group were significantly greater than those in the skeletal muscle satellite stem cell group (Figure 7B). Such results implied that Tbx2 was upregulated during skeletal muscle development, which was consistent with our qPCR results.
Figure 7. METTL3/m6A/Brip1os/Tbx2 axis in muscle development. (A) The genomic location of lncRNA Brip1os and its adjacent mRNA Tbx2. (B) Tbx2 expression data from public record of single cell RNA-seq gene expression of Tabula Muris data. (C) Brip1os and Tbx2 expression levels in our RNA-seq data. Pearson correlation analysis was performed to estimate expression correlation between two genes. (D) qRT-PCR shows the expression of Tbx2 when METTL3 overexpressed. (E) qRT-PCR shows the expression of Tbx2 when METTL3 was knockdown. Data are presented as Mean ± SD; p value: **p < 0.01, ****p < 0.0001.
As shown in Figure 7C, Brip1os was highly expressed in two GM samples and decreased in two D4 samples, in which Tbx2 exhibited the opposite trend (Pearson R = −0.99, p = 0.009). These results confirm their negative correlation with RNA expression. Next, we assessed whether METTL3 affects Tbx2 expression. As shown in Figure 7D, the Tbx2 mRNA levels were greatly downregulated when METTL3 was overexpressed. Accordingly, Tbx2 mRNA levels were upregulated when METTL3 was knocked down (Figure 7E). These results suggest that Tbx2 could be regulated by METTL3. Taken together, we validated the retrograde regulatory relationship between Brip1os and Tbx2 as well as their responsive reaction to changes in the m6A methyltransferase METTL3.
Discussion
Skeletal muscle development is precisely regulated in a sophisticated spatiotemporal manner. Our previously identified changes in gene expression and epigenetic modifications during skeletal muscle development have greatly improved our understanding of the mechanism related to myogenesis, including coding gene and non-coding RNA modifications (Diao et al., 2021a, b), mRNA expression (Tan et al., 2021), miRNA regulation (Xie et al., 2013, 2018), and lncRNA function (Liu et al., 2021). Taking advantage of the sequencing approach and gene annotation, a significant number of lncRNAs have been shown to play crucial roles in skeletal muscle development (Luo et al., 2021). Given the robust function of m6A methylation, the functions of m6A-modified mRNAs in the process of skeletal muscle development have been well studied, and the role of m6A-modified non-coding RNAs has also been appreciated (Li et al., 2021). In the present study, we hypothesized that lncRNAs might also be modified by m6A and participate in skeletal muscle differentiation. In the current study, we provided the first evidence that both the lncRNA transcriptome and m6A epitranscriptome underwent highly dynamic changes throughout mouse skeletal muscle development. Such results are consistent with results from other studies investigating the role of m6A-lncRNA in tissue development, such as those for mouse embryonic stem cell differentiation (Wang et al., 2014). Specifically, we observed that m6A methylation levels of lncRNAs are positively correlated with the transcript abundance of lncRNAs. Moreover, our studies revealed that lncRNAs exhibit pairwise expression correlations with neighboring mRNAs. Our results highlight a potential role of m6A-modified lncRNAs during skeletal muscle development.
LncRNAs are more cell-specific than other RNAs, and their expression models are not completely understood. Given that m6A is a ubiquitous modification in RNAs and regulates gene expression, we systematically identified m6A-modified lncRNAs and uncovered the m6A marks affecting the expression of lncRNAs. Our data showed that for lncRNAs with significant expression abundance changes, their m6A levels were positively correlated with their expression levels. However, for lncRNAs without significant transcript abundance changes, no significant correlation was noted between their m6A levels and expression levels. These results were further verified by overexpression or knockdown of the m6A core methyltransferase METTL3.
Due to the poor conservation of lncRNAs, their function and regulatory mechanisms are not completely understood. It is known that lncRNAs can regulate the expression of neighboring genes by cis-acting mechanisms (Lee, 2012). Mancini-DiNardo et al. clarified that elongation of the Kcnq1ot1 transcript is required for genomic imprinting of neighboring genes (Mancini-DiNardo et al., 2006). Furthermore, Ponjavic et al. noted that spatiotemporal coexpression of ncRNAs and nearby protein-coding genes represents a general phenomenon and presented substantive and predictive criteria for prioritizing lncRNA and mRNA transcript pairs when investigating their biological functions (Ponjavic et al., 2009). This information led us to hypothesize that m6A-methylated lncRNAs regulate nearby mRNAs and contribute to muscle tissue development. This hypothesis was supported by our computational analysis and experimental results. In particular, Brip1os is the most significantly differentially expressed and m6A-modified lncRNA, and its nearby mRNA Tbx2 plays an important function in muscle tissue development. The change in Tbx2 mRNA expression was opposite to that of Brip1os in D4 compared to GM. These findings indicate that they exhibit a retrograde regulatory relationship, which was verified using a public dataset. In addition, it has been known for decades that Tbx (T-Box) genes play crucial roles in limb development (Zhu et al., 2014; Pflugfelder et al., 2017). Further studies on their regulation are warranted. Our analysis provided candidate lncRNAs and mRNAs for further examination of the gene regulation network in muscle development.
The transition from myoblast proliferation to differentiation is accompanied by drastic alterations in the transcriptome. Transcriptional changes influencing muscle status are affected by a number of processes involving DNA, RNA and proteins. This study uncovered the differential expression and m6A methylation status of lncRNAs with temporal-specific expression in developing muscle. Surprisingly, we found no differences in the methylation of lncRNAs that drive myoblast state changes, such as linc-MD1 (Cesana et al., 2011) or lncMyoD (Dong et al., 2020), suggesting that these lncRNAs are not direct targets of dynamic m6A modification. Kcnq1ot1, a differentially expressed and m6A-enriched lncRNA identified in our study, participates in the regulation of genes within the Kcnq1 imprinting domain (Zhang et al., 2014) and controls maternal p57 expression in muscle cells by promoting H3K27me3 accumulation in an intragenic MyoD-binding region (Andresini et al., 2019). Intriguingly, of the top m6A-enriched lncRNA transcripts in GM and in D4, many have not been previously linked to myoblast/myotube function but have been reported to be involved in other cell types. For example, Snhg14 (small nucleolar RNA host gene 14) is expressed at elevated levels among the top ten differentially m6A-enriched lncRNAs in D4 compared to GM. Snhg14 is highly expressed in Parkinson’s disease, and silencing Sngh14 mitigated dopaminergic neuron injury by downregulating a-syn by targeting miR-133b (Zhang et al., 2019). In addition, Sngh14 functions as a ceRNA in Ang II-induced cardiomyocytes to sponge both miR-322-5p and miR-384-5p to elevate PCDH17 levels (Long et al., 2020). The function of m6A-modified lncRNAs during skeletal muscle development needs to be further studied.
In summary, we described the expression and m6A methylation profiles of lncRNAs that display temporal expression in mouse myoblasts and differentiated myotubes. Our findings provide new insight into the pivotal regulatory role of m6A-modified lncRNAs in muscle development. Our data uncovered the novel posttranscriptional regulation underlying muscle differentiation and provide a molecular basis for further studies to determine the function and mechanism of m6A-lncRNAs in skeletal muscle development and muscle-related diseases.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Ethics Statement
The animal study was reviewed and approved by the Institutional Animal Care and Use Committee of Sun Yat-sen University.
Author Contributions
S-JX, Z-DX, QZ, and Y-WP conceived and designed this project. S-JX, ST, and L-TD performed the experiments with the help of Y-XH, Y-RH, HL and, W-YX. S-JX, W-JC, Z-GZ, P-LL, and W-CC analyzed the data. S-JX, Z-DX, ST, L-TD, Y-WP, and QZ wrote the manuscript. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by grants from the Natural Science Foundation of Guangdong Province (2021A1515010477), the National Natural Science Foundation of China (31701116, 81974436, 81900559, 81870449, 81770651, 82170674, and 82170216), the Program for Guangdong Introducing Innovative and Entrepreneurial Teams (2019ZT08Y485), Guangzhou Science and Technology Plan Projects (202102020104), the Key-Area Research and Development Program of Guangdong Province (2019B020236004 and 2019B020235002), and the Guangdong Province Key Laboratory of Malignant Tumor Epigenetics and Gene Regulation (2020B1212060018).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
The authors would like to thank Zelin Wang (Shuzhi Biotech, LLC, Guangzhou) who assisted our efforts in analyzing the sequencing data used in this study.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.762669/full#supplementary-material
Footnotes
References
Alarcón, C. R., Lee, H., Goodarzi, H., Halberg, N., and Tavazoie, S. F. (2015). N6-methyladenosine marks primary microRNAs for processing. Nature 519, 482–485. doi: 10.1038/nature14281
Andresini, O., Rossi, M. N., Matteini, F., Petrai, S., Santini, T., and Maione, R. (2019). The long non-coding RNA Kcnq1ot1 controls maternal p57 expression in muscle cells by promoting H3K27me3 accumulation to an intragenic MyoD-binding region. Epigenetics Chromatin 12:8. doi: 10.1186/s13072-019-0253-1
Bentzinger, C. F., Wang, Y. X., and Rudnicki, M. A. (2012). Building muscle: molecular regulation of myogenesis. Cold Spring Harb. Perspect. Biol. 4:a008342. doi: 10.1101/cshperspect.a008342
Borensztein, M., Monnier, P., Court, F., Louault, Y., Ripoche, M. A., Tiret, L., et al. (2013). Myod and H19-Igf2 locus interactions are required for diaphragm formation in the mouse. Development 140, 1231–1239. doi: 10.1242/dev.084665
Bryson-Richardson, R. J., and Currie, P. D. (2008). The genetics of vertebrate myogenesis. Nat. Rev. Genet. 9, 632–646. doi: 10.1038/nrg2369
Cesana, M., Cacchiarelli, D., Legnini, I., Santini, T., Sthandier, O., Chinappi, M., et al. (2011). A long noncoding RNA controls muscle differentiation by functioning as a competing endogenous RNA. Cell 147, 358–369. doi: 10.1016/j.cell.2011.09.028
Chen, L. L., and Carmichael, G. G. (2010). Decoding the function of nuclear long non-coding RNAs. Curr. Opin. Cell Biol. 22, 357–364. doi: 10.1016/j.ceb.2010.03.003
Chen, T., Hao, Y. J., Zhang, Y., Li, M. M., Wang, M., Han, W., et al. (2015). m(6)A RNA methylation is regulated by microRNAs and promotes reprogramming to pluripotency. Cell Stem Cell 16, 289–301. doi: 10.1016/j.stem.2015.01.016
Cong, X. X., Gao, X. K., Rao, X. S., Wen, J., Liu, X. C., Shi, Y. P., et al. (2020). Rab5a activates IRS1 to coordinate IGF-AKT-mTOR signaling and myoblast differentiation during muscle regeneration. Cell Death Differ. 27, 2344–2362. doi: 10.1038/s41418-020-0508-1
Cui, X., Wei, Z., Zhang, L., Liu, H., Sun, L., Zhang, S. W., et al. (2016). Guitar: an R/Bioconductor package for gene annotation guided transcriptomic analysis of RNA-related genomic features. Biomed. Res. Int. 2016:8367534. doi: 10.1155/2016/8367534
Diao, L. T., Xie, S. J., Lei, H., Qiu, X. S., Huang, M. C., Tao, S., et al. (2021a). METTL3 regulates skeletal muscle specific miRNAs at both transcriptional and post-transcriptional levels. Biochem. Biophys. Res. Commun. 552, 52–58. doi: 10.1016/j.bbrc.2021.03.035
Diao, L. T., Xie, S. J., Yu, P. J., Sun, Y. J., Yang, F., Tan, Y. Y., et al. (2021b). N-methyladenine demethylase ALKBH1 inhibits the differentiation of skeletal muscle. Exp. Cell Res. 400:112492. doi: 10.1016/j.yexcr.2021.112492
Dong, A., Preusch, C. B., So, W.-K., Lin, K., Luan, S., Yi, R., et al. (2020). A long noncoding RNA, LncMyoD, modulates chromatin accessibility to regulate muscle stem cell myogenic lineage progression. Proc. Natl. Acad. Sci. 117, 32464–32475. doi: 10.1073/pnas.2005868117
Engreitz, J. M., Haines, J. E., Perez, E. M., Munson, G., Chen, J., Kane, M., et al. (2016). Local regulation of gene expression by lncRNA promoters, transcription and splicing. Nature 539, 452–455. doi: 10.1038/nature20149
Fazi, F., and Fatica, A. (2019). Interplay between N (6)-methyladenosine (m(6)A) and non-coding RNAs in cell development and cancer. Front. Cell Dev. Biol. 7:116. doi: 10.3389/fcell.2019.00116
Fu, Y., Dominissini, D., Rechavi, G., and He, C. (2014). Gene expression regulation mediated through reversible m(6)A RNA methylation. Nat. Rev. Genet. 15, 293–306. doi: 10.1038/nrg3724
Gheller, B. J., Blum, J. E., Fong, E. H. H., Malysheva, O. V., Cosgrove, B. D., and Thalacker-Mercer, A. E. (2020). A defined N6-methyladenosine (m6A) profile conferred by METTL3 regulates muscle stem cell/myoblast state transitions. Cell Death Discov. 6:95. doi: 10.1038/s41420-020-00328-5
Gong, C., Li, Z., Ramanujan, K., Clay, I., Zhang, Y., Lemire-Brachat, S., et al. (2015). A long non-coding RNA, LncMyoD, regulates skeletal muscle differentiation by blocking IMP2-mediated mRNA translation. Dev. Cell 34, 181–191. doi: 10.1016/j.devcel.2015.05.009
He, R. Z., Jiang, J., and Luo, D. X. (2020). The functions of N6-methyladenosine modification in lncRNAs. Genes Dis. 7, 598–605. doi: 10.1016/j.gendis.2020.03.005
Heinz, S., Benner, C., Spann, N., Bertolino, E., Lin, Y. C., Laslo, P., et al. (2010). Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol. Cell 38, 576–589. doi: 10.1016/j.molcel.2010.05.004
Helm, M., and Motorin, Y. (2017). Detecting RNA modifications in the epitranscriptome: predict and validate. Nat. Rev. Genet. 18, 275–291. doi: 10.1038/nrg.2016.169
Hitachi, K., Nakatani, M., Takasaki, A., Ouchi, Y., Uezumi, A., Ageta, H., et al. (2019). Myogenin promoter-associated lncRNA Myoparr is essential for myogenic differentiation. EMBO Rep. 20:e47468. doi: 10.15252/embr.201847468
Kan, R. L., Chen, J., and Sallam, T. (2021). Crosstalk between epitranscriptomic and epigenetic mechanisms in gene regulation. Trends Genet. 19:S0168-9525(21)00170-0. doi: 10.1016/j.tig.2021.06.014
Kechin, A., Boyarskikh, U., Kel, A., and Filipenko, M. (2017). cutPrimers: a new tool for accurate cutting of primers from reads of targeted next generation sequencing. J. Comput. Biol. 24, 1138–1143. doi: 10.1089/cmb.2017.0096
Kim, D., Paggi, J. M., Park, C., Bennett, C., and Salzberg, S. L. (2019). Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat. Biotechnol. 37, 907–915. doi: 10.1038/s41587-019-0201-4
Kudou, K., Komatsu, T., Nogami, J., Maehara, K., Harada, A., Saeki, H., et al. (2017). The requirement of Mettl3-promoted MyoD mRNA maintenance in proliferative myoblasts for skeletal muscle differentiation. Open Biol. 7:170119. doi: 10.1098/rsob.170119
Lan, Y., Liu, B., and Guo, H. (2021). The role of m6A modification in the regulation of tumor-related lncRNAs. Mol. Ther. Nucleic Acids 24, 768–779. doi: 10.1016/j.omtn.2021.04.002
Lee, J. T. (2012). Epigenetic regulation by long noncoding RNAs. Science 338, 1435–1439. doi: 10.1126/science.1231776
Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., et al. (2009). The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079. doi: 10.1093/bioinformatics/btp352
Li, J., Pei, Y., Zhou, R., Tang, Z., and Yang, Y. (2021). Regulation of RNA N6-methyladenosine modification and its emerging roles in skeletal muscle development. Int. J. Biol. Sci. 17, 1682–1692. doi: 10.7150/ijbs.56251
Liang, Y., Han, H., Xiong, Q., Yang, C., Wang, L., Ma, J., et al. (2021). METTL3-mediated m6A methylation regulates muscle stem cells and muscle regeneration by notch signaling pathway. Stem Cells Int. 2021:9955691. doi: 10.1155/2021/9955691
Liao, Y., Smyth, G. K., and Shi, W. (2014). featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30, 923–930. doi: 10.1093/bioinformatics/btt656
Lim, Y. H., Ryu, J., Kook, H., and Kim, Y. K. (2020). Identification of long noncoding RNAs involved in differentiation and survival of vascular smooth muscle cells. Mol. Ther. Nucleic Acids 22, 209–221. doi: 10.1016/j.omtn.2020.08.032
Lin, J., Zhu, Q., Huang, J., Cai, R., and Kuang, Y. (2020). Hypoxia promotes vascular smooth muscle cell (VSMC) differentiation of adipose-derived stem cell (ADSC) by regulating Mettl3 and paracrine factors. Stem Cells Int. 2020:2830565. doi: 10.1155/2020/2830565
Liu, J., Li, K., Cai, J., Zhang, M., Zhang, X., Xiong, X., et al. (2020). Landscape and regulation of m(6)A and m(6)Am methylome across human and mouse tissues. Mol. Cell 77, 426–440.e6. doi: 10.1016/j.molcel.2019.09.032
Liu, N., Dai, Q., Zheng, G., He, C., Parisien, M., and Pan, T. (2015). N6-methyladenosine-dependent RNA structural switches regulate RNA–protein interactions. Nature 518, 560–564. doi: 10.1038/nature14234
Liu, S., Xie, S., Chen, H., Li, B., Chen, Z., Tan, Y., et al. (2021). The functional analysis of transiently upregulated miR-101 suggests a “braking” regulatory mechanism during myogenesis. Sci. China Life Sci. 64, 1–12. doi: 10.1007/s11427-020-1856-5
Long, Y., Wang, L., and Li, Z. (2020). SP1-induced SNHG14 aggravates hypertrophic response in in vitro model of cardiac hypertrophy via up-regulation of PCDH17. J. Cell. Mol. Med. 24, 7115–7126. doi: 10.1111/jcmm.15073
Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15:550. doi: 10.1186/s13059-014-0550-8
Lu, L., Sun, K., Chen, X., Zhao, Y., Wang, L., Zhou, L., et al. (2013). Genome-wide survey by ChIP-seq reveals YY1 regulation of lincRNAs in skeletal myogenesis. EMBO J. 32, 2575–2588. doi: 10.1038/emboj.2013.182
Luo, H., Lv, W., Tong, Q., Jin, J., Xu, Z., and Zuo, B. (2021). Functional non-coding RNA during embryonic myogenesis and postnatal muscle development and disease. Front. Cell Dev. Biol. 9: 628339. doi: 10.3389/fcell.2021.628339
Lv, W., Jin, J., Xu, Z., Luo, H., Guo, Y., Wang, X., et al. (2020). lncMGPF is a novel positive regulator of muscle growth and regeneration. J. Cachexia Sarcopenia Muscle 11, 1723–1746. doi: 10.1002/jcsm.12623
Ma, S., Chen, C., Ji, X., Liu, J., Zhou, Q., Wang, G., et al. (2019). The interplay between m6A RNA methylation and noncoding RNA in cancer. J. Hematol. Oncol. 12:121. doi: 10.1186/s13045-019-0805-7
Mancini-DiNardo, D., Steele, S. J. S., Levorse, J. M., Ingram, R. S., and Tilghman, S. M. (2006). Elongation of the Kcnq1ot1 transcript is required for genomic imprinting of neighboring genes. Gene Dev. 20, 1268–1282. doi: 10.1101/gad.1416906
Martinet, C., Monnier, P., Louault, Y., Benard, M., Gabory, A., and Dandolo, L. (2016). H19 controls reactivation of the imprinted gene network during muscle regeneration. Development 143, 962–971. doi: 10.1242/dev.131771
Martone, J., Mariani, D., Desideri, F., and Ballarino, M. (2020). Non-coding RNAs shaping muscle. Front. Cell Dev. Biol. 7:394. doi: 10.3389/fcell.2019.00394
Meng, J., Cui, X., Rao, M. K., Chen, Y., and Huang, Y. (2013). Exome-based analysis for RNA epigenome sequencing data. Bioinformatics 29, 1565–1567. doi: 10.1093/bioinformatics/btt171
Meyer, K. D., and Jaffrey, S. R. (2017). Rethinking m(6)A readers, writers, and erasers. Annu. Rev. Cell Dev. Biol. 33, 319–342. doi: 10.1146/annurev-cellbio-100616-060758
Milligan, L., Antoine, E., Bisbal, C., Weber, M., Brunel, C., Forne, T., et al. (2000). H19 gene expression is up-regulated exclusively by stabilization of the RNA during muscle cell di€erentiation. Oncogene 19, 5810–5816.
Pan, T. (2013). N6-methyl-adenosine modification in messenger and long non-coding RNA. Trends Biochem. Sci. 38, 204–209. doi: 10.1016/j.tibs.2012.12.006
Patil, D. P., Chen, C. K., Pickering, B. F., Chow, A., Jackson, C., Guttman, M., et al. (2016). m(6)A RNA methylation promotes XIST-mediated transcriptional repression. Nature 537, 369–373. doi: 10.1038/nature19342
Pflugfelder, G. O., Eichinger, F., and Shen, J. (2017). T-Box genes in Drosophila limb development. Curr. Top. Dev. Biol. 122, 313–354. doi: 10.1016/bs.ctdb.2016.08.003
Ping, X. L., Sun, B. F., Wang, L., Xiao, W., Yang, X., Wang, W. J., et al. (2014). Mammalian WTAP is a regulatory subunit of the RNA N6-methyladenosine methyltransferase. Cell Res. 24, 177–189. doi: 10.1038/cr.2014.3
Ponjavic, J., Oliver, P. L., Lunter, G., and Ponting, C. P. (2009). Genomic and transcriptional co-localization of protein-coding and long non-coding RNA pairs in the developing brain. PLoS Genet. 5:e1000617. doi: 10.1371/journal.pgen.1000617
Quinlan, A. R., and Hall, I. M. (2010). BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841–842. doi: 10.1093/bioinformatics/btq033
Ramírez, F., Dündar, F., Diehl, S., Grüning, B. A., and Manke, T. (2014). deepTools: a flexible platform for exploring deep-sequencing data. Nucleic Acids Res. 42, W187–W191. doi: 10.1093/nar/gku365
Ro, S., Lim, Y.-H., Kwon, D.-H., Kim, J., Park, W. J., Kook, H., et al. (2018). Identification of long noncoding RNAs involved in muscle differentiation. PLos One 13:e0193898. doi: 10.1371/journal.pone.0193898
Robinson, J. T., Thorvaldsdóttir, H., Winckler, W., Guttman, M., Lander, E. S., Getz, G., et al. (2011). Integrative genomics viewer. Nat. Biotechnol. 29, 24–26. doi: 10.1038/nbt.1754
Roundtree, I. A., Evans, M. E., Pan, T., and He, C. (2017). Dynamic RNA modifications in gene expression regulation. Cell 169, 1187–1200. doi: 10.1016/j.cell.2017.05.045
Shi, H., Wei, J., and He, C. (2019). Where, when, and how: context-dependent functions of RNA methylation writers, readers, and erasers. Mol. Cell 74, 640–650. doi: 10.1016/j.molcel.2019.04.025
Sui, Y., Han, Y., Zhao, X., Li, D., and Li, G. (2019). Long non-coding RNA Irm enhances myogenic differentiation by interacting with MEF2D. Cell Death Dis. 10:181. doi: 10.1038/s41419-019-1399-2
Sweta, S., Dudnakova, T., Sudheer, S., Baker, A. H., and Bhushan, R. (2019). Importance of long non-coding RNAs in the development and disease of skeletal muscle and cardiovascular lineages. Front. Cell Dev. Biol. 7:228. doi: 10.3389/fcell.2019.00228
Tan, Y.-Y., Zhang, Y., Li, B., Ou, Y.-W., Xie, S.-J., Chen, P.-P., et al. (2021). PERK signaling controls myoblast differentiation by regulating microRNA networks. Front. Cell Dev. Biol. 9:670435. doi: 10.3389/fcell.2021.670435
Ulitsky, I., and Bartel, D. P. (2013). lincRNAs: genomics, evolution, and mechanisms. Cell 154, 26–46. doi: 10.1016/j.cell.2013.06.020
Wang, L., Zhao, Y., Bao, X., Zhu, X., Kwok, Y.K.-y, Sun, K., et al. (2015). LncRNA dum interacts with dnmts to regulate Dppa2 expression during myogenic differentiation and muscle regeneration. Cell Res. 25, 335–350. doi: 10.1038/cr.2015.21
Wang, X., Huang, N., Yang, M., Wei, D., Tai, H., Han, X., et al. (2017). FTO is required for myogenesis by positively regulating mTOR-PGC-1alpha pathway-mediated mitochondria biogenesis. Cell Death Dis. 8:e2702. doi: 10.1038/cddis.2017.122
Wang, Y., Li, Y., Toth, J. I., Petroski, M. D., Zhang, Z., and Zhao, J. C. (2014). N6-methyladenosine modification destabilizes developmental regulators in embryonic stem cells. Nat. Cell Biol. 16, 191–198. doi: 10.1038/ncb2902
Xie, S., Zhang, Y., Qu, L., and Xu, H. (2013). A helm model for microRNA regulation in cell fate decision and conversion. Sci. China Life Sci. 56, 897–906. doi: 10.1007/s11427-013-4547-4
Xie, S. J., Diao, L. T., Cai, N., Zhang, L. T., Xiang, S., Jia, C. C., et al. (2021a). mascRNA and its parent lncRNA MALAT1 promote proliferation and metastasis of hepatocellular carcinoma cells by activating ERK/MAPK signaling pathway. Cell Death Discov. 7:110. doi: 10.1038/s41420-021-00497-x
Xie, S.-J., Lei, H., Yang, B., Diao, L.-T., Liao, J.-Y., He, J.-H., et al. (2021b). Dynamic m6A mRNA methylation reveals the role of METTL3/14-m6A-MNK2-ERK signaling axis in skeletal muscle differentiation and regeneration. Front. Cell Dev. Biol. 38, 4755–4772. doi: 10.3389/fcell.2021.744171
Xie, S.-J., Li, J.-H., Chen, H.-F., Tan, Y.-Y., Liu, S.-R., Zhang, Y., et al. (2018). Inhibition of the JNK/MAPK signaling pathway by myogenesis-associated miRNAs is required for skeletal muscle development. Cell Death Differ. 25, 1581–1597. doi: 10.1038/s41418-018-0063-1
Yang, D., Qiao, J., Wang, G., Lan, Y., Li, G., Guo, X., et al. (2018). N6-Methyladenosine modification of lincRNA 1281 is critically required for mESC differentiation potential. Nucleic Acids Res. 46, 3906–3920. doi: 10.1093/nar/gky130
Yang, Y., Hsu, P. J., Chen, Y. S., and Yang, Y. G. (2018). Dynamic transcriptomic m(6)A decoration: writers, erasers, readers and functions in RNA metabolism. Cell Res. 28, 616–624. doi: 10.1038/s41422-018-0040-8
Yong, H., Wu, G., Chen, J., Liu, X., Bai, Y., Tang, N., et al. (2020). lncRNA MALAT1 accelerates skeletal muscle cell apoptosis and inflammatory response in sepsis by decreasing BRCA1 expression by recruiting EZH2. Mol. Ther. Nucleic Acids 19, 97–108. doi: 10.1016/j.omtn.2019.10.028
Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. Omics 16, 284–287. doi: 10.1089/omi.2011.0118
Yu, X., Zhang, Y., Li, T., Ma, Z., Jia, H., Chen, Q., et al. (2017). Long non-coding RNA Linc-RAM enhances myogenic differentiation by interacting with MyoD. Nat. Commun. 10:181. doi: 10.1038/ncomms14016
Zhang, H., Shi, X., Huang, T., Zhao, X., Chen, W., Gu, N., et al. (2020). Dynamic landscape and evolution of m6A methylation in human. Nucleic Acids Res. 48, 6251–6264. doi: 10.1093/nar/gkaa347
Zhang, H., Zeitz, M. J., Wang, H., Niu, B., Ge, S., Li, W., et al. (2014). Long noncoding RNA-mediated intrachromosomal interactions promote imprinting at the Kcnq1 locus. J. Cell Biol. 204, 61–75. doi: 10.1083/jcb.201304152
Zhang, L. M., Wang, M. H., Yang, H. C., Tian, T., Sun, G. F., Ji, Y. F., et al. (2019). Dopaminergic neuron injury in Parkinson’s disease is mitigated by interfering lncRNA SNHG14 expression to regulate the miR-133b/α-synuclein pathway. Aging 11, 9264–9279. doi: 10.18632/aging.102330
Zhang, W., He, L., Liu, Z., Ren, X., Qi, L., Wan, L., et al. (2020). Multifaceted functions and novel insight into the regulatory role of RNA N(6)-Methyladenosine modification in musculoskeletal disorders. Front. Cell Dev. Biol. 8:870. doi: 10.3389/fcell.2020.00870
Zhang, Y., Li, Y., Hu, Q., Xi, Y., Xing, Z., Zhang, Z., et al. (2020). The lncRNA H19 alleviates muscular dystrophy by stabilizing dystrophin. Nat. Cell Biol. 22, 1332–1345. doi: 10.1038/s41556-020-00595-5
Zhao, Y., Cao, F., Yu, X., Chen, C., Meng, J., Zhong, R., et al. (2018). Linc-RAM is required for FGF2 function in regulating myogenic cell differentiation. RNA Biol. 15, 404–412. doi: 10.1080/15476286.2018.1431494
Zhou, L., Sun, K., Zhao, Y., Zhang, S., Wang, X., Li, Y., et al. (2015). Linc-YY1 promotes myogenic differentiation and muscle regeneration through an interaction with the transcription factor YY1. Nat. Commun. 6:10026. doi: 10.1038/ncomms10026
Keywords: m6A, lncRNAs, Brip1os, METTL3, skeletal muscle development
Citation: Xie S-J, Tao S, Diao L-T, Li P-L, Chen W-C, Zhou Z-G, Hu Y-X, Hou Y-R, Lei H, Xu W-Y, Chen W-J, Peng Y-W, Zhang Q and Xiao Z-D (2021) Characterization of Long Non-coding RNAs Modified by m6A RNA Methylation in Skeletal Myogenesis. Front. Cell Dev. Biol. 9:762669. doi: 10.3389/fcell.2021.762669
Received: 22 August 2021; Accepted: 14 September 2021;
Published: 13 October 2021.
Edited by:
Huilin Huang, Sun Yat-sen University Cancer Center (SYSUCC), ChinaReviewed by:
Zhanpeng Huang, The First Affiliated Hospital of Sun Yat-sen University, ChinaChao Shen, City of Hope National Medical Center, United States
Copyright © 2021 Xie, Tao, Diao, Li, Chen, Zhou, Hu, Hou, Lei, Xu, Chen, Peng, Zhang and Xiao. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Zhen-Dong Xiao, xiaozhd@mail2.sysu.edu.cn; Qi Zhang, zhangq27@mail.sysu.edu.cn; Yan-Wen Peng, pengyw@mail.sysu.edu.cn
†These authors have contributed equally to this work