J-1 male mouse embryonic stem cells (ESCs) were purchased from the Cell Bank of Chinese Academy of Sciences. Wild-type and Mettl3-KO J-1 mouse ESCs were provided by Howard Y. Chang in Stanford University (Batista et al., 2014). Ythdc1 ^flox/flox (Control) and Ythdc1-cKO (Ythdc1-KO) mouse ESCs were provided by Yawei Gao in Tongji University (Chen et al., 2021). Mycoplasma testing was performed routinely and cells tested negative were used in further experiments. Tissue culture plates were initially treated with 0.2% gelatin for at least 30 min at 37°C. ESCs were subsequently grown in gelatinized plates containing 15% KnockOut Serum Replacement (Gibco), 1% non-essential amino acids (Gibco), 1% Glutamax (Gibco), 1% Penn-Strep (Life technologies), 0.1 mM β-mercaptoethanol (Invitrogen) and 1000 U/ml ESGRO Recombinant Mouse leukemia inhibitory factor (LIF) Protein (Millipore).

To generate embryoid bodies (EBs), suspension of each cell line at a concentration of 10 × 10⁴cells per milliliter in differentiating medium (Knock-out DMEM supplemented with 15% FBS (Gibco), 1% non-essential amino acids, 1% Glutamax, 1% Penn-Strep, 0.1 mm β-mercaptoethanol) was deposited in 20 μL hanging drops in Petri dishes. After 2 days, the suspension-cultured cells were plated onto 0.2% gelatin-coated dishes in cardiomyocyte differentiation media (CMD) containing DMEM (BI), 15% FBS (Gibco), 1% Penn-Strep, 1% non-essential amino acids (Gibco), 1% Glutamine (Gibco), 0.1 mm β-mercaptoethanol and 1 mm Ascorbic Acid (Sigma) for continued differentiation. Cell culture media was changed every other day and contracting patches of cells were collected at indicated time points.

The shRNA targeting sequences were cloned into pRSI9-U6-(sh)-UbiC-TagRFP-2A-Blasticidin. Lentiviral particles were packaged using Lenti-X 293T cells. Virus-containing supernatants were collected at 48 h after transfection and filtered to eliminate cell debris. Wild-type and Mettl3-KO mESCs were infected by the lentivirus for 48 h and stable mESC knockdown lines were generated using blasticidin (MCE) selection (10 µg/ml). shRNA targeting sequences are listed below.

Gm2379 (target sequence 1): 5′- GCA​CGG​CAT​TAC​AGA​ATT​TCA-3′; Gm2379 (target sequence 3): 5′- GAA​AGC​CTT​CAT​AGG​AAA​CAG-3′; Scramble Target Sequence: 5′- AAT​TCT​CCG​AAC​GTG​TCA​CGT-3′.

Total RNA from different cells was isolated from PBS-washed cells using Trizol reagent (Invitrogen). Isolated RNAs were used as templates for reverse transcription based on the manual for HiScript III 1st Strand cDNA Synthesis Kit plus gDNA wiper (Vazyme Biotech). Real-time PCR analysis was conducted with ChamQ SYBR qPCR Master Mix (Vazyme Biotech) and carried on a QuantStudio three Real-Time PCR System (Applied Biosystems). Primers for amplifying each target were listed in Supplementary Table S1.

ESCs were pelleted by centrifugation at 500 g for 3 min at 4°C. After washing with phosphate-buffered saline (PBS) buffer twice, cells were lysed in ice-cold cytoplasmic lysis buffer (10 mM HEPES pH 7.5, 0.05%TritonX-100, 5 mM MgCl₂, 10 U SUPERase-In and 1x protease inhibitor cocktail) on ice for 5 min. Samples were centrifuged at 2,000 g at 4°C for 10 min. The supernatant was collected as cytoplasmic fraction. The remaining nuclear pellets were washed with PBS twice and centrifuged at 2,000 g at 4°C for 2 min. After removing the supernatant, the nuclear pellet was resuspended in RIPA buffer (50 mM Tris-HCl, pH7.5, 150 mM NaCl, 1% NP-40, 0.5% sodium deoxycholate, 10 U SUPERase-In and 1 × protease inhibitor cocktail). This sample was collected as nuclear fraction. RNA from different fractions was isolated using TRIzol reagent (Invirtrogen), and the protein sample was prepared by adding SDS loading buffer.

Total RNA was isolated using TRIzol reagent (Invitrogen). Ribosomal RNAs depletion, strand-specific library preparation and sequencing were performed at the Annoroad Gene Technology (Beijing, China). Two biological replicates were prepared for each differentiation phases. The libraries were sequenced on an Illumina Novaseq 6,000 system using the 2 × 150-bp paired-end read mode.

The lncRNA-seq data has been deposited in NCBI’s Gene Expression Omnibus and is accessible through GEO Series accession number GSE197205. Trim Galore (v.0.6.6) was used to trim the low quality sequences (Martin, 2011). The paired-end reads were aligned to mouse genome (GRCm38) with HISAT2 aligner (v2.2.1) with --rna-strandness RF parameters and Ensemble genes (GRCm38.102) transcriptome annotation (Pertea et al., 2016). Then gene expression were quantified using featureCounts (v2.0.1) with -s 2 parameter (Liao et al., 2014), and differentially expressed genes were identifified using DESeq2 (v1.34.0) R package (Love et al., 2014). The FoldChange > 2 or < 0.5 and ajusted p-value < 0.05 genes were considered as differentially expressed genes. The TPM (Transcripts Per Kilobase of exonmodel per Million mapped reads) was used to produce heatmap.

LncRNAs with average TPM expressions in replicates from four differentiation stages were grouped into six different clusters using Mfuzz (v2.54.0) package in R with fuzzy c-means algorithm (Kumar and Futschik, 2007).

HOMER’s findMotifs.pl program was used to identify motifs in the promoters of lncRNA genes (encompassing regions 2000 bp upstream and 500 bp downstream of the transcription start site) from each cluster. If one lncRNA has multiple transcript isoforms, the longest transcript isoform was used for analysis and bedtools getfasta (v2.29.2) was used to extract sequnences. Enrichment analysis was performed on Gene Ontology website (Ashburner et al., 2000). Network plot was produced by Cytoscape (v3.8.2) (Paul Shannon et al., 1971).

MeRIP-seq and RIP-seq datasets were obtained from the Gene Expression Omnibus (GEO) public database under accession number GSE145315 and GSE157265, respectively (Chen et al., 2021; Liu et al., 2021). Quality control of MeRIP-seq and RIP-seq data were performed using FastQC (v.0.11.9). Low-quality bases and adaptor contamination were removed by Trim Galore. Remaining reads were aligned to mouse genome version GRCm38 with HISAT2 aligner with default parameters. For RIP-seq data SAMtools (v1.13) was used to remove the reads with mapping quality below 30 (Li et al., 2009). m⁶A peaks and differential m⁶A peaks were identified using exomePeak2 (v1.6.1) R package with p_cutoff = 0.00001 and log2FC_cutoff = 1 parameters (Meng et al., 2014). Log₂FoldChange < −1 and ajusted p value < 0.05 was considered as significantly downregulated peaks. For biological replicates, bedtools (v2.29.2) was used to intersect peaks with overlaped proportion > 0.5 (Meng et al., 2014). Peaks were annotated by HOMER’s annotatePeaks.pl program. The distribution of m6A peaks was plotted by R package Guitar (v.1.5.0) (Cui et al., 2016). The final bam files were transformed to bigwig format by using Deeptools bamCoverage (v3.5.1) with --normalizeUsing RPKM parameters (Ramírez et al., 2014). IGV (Integrative genomics viewer, v2.11.1) was used to visualize the aligned data (Thorvaldsdóttir et al., 2013). DESeq2 was used to identify Ythdc1 binding genes. The FoldChange > 2 and ajusted p-value < 0.05 genes were considered as Ythdc1 binding RNAs.

Ythdc1‐KO and Mettl3‐KO RNA-seq datasets were obtained from the Gene Expression Omnibus (GEO) public database under accession number GSE157262. The expression analysis was performed based on the available processed fpkm-matrix files in the database.

To investigate the global expression profiles of lncRNA transcriptome during ESC-derived cardiac differentiation, we collected the samples at four different phases including undifferentiated ESC, mesoderm (MES), cardiac progenitor (CP), and CM. Purified RNA after removal of ribosomal RNA was utilized for sequencing. The reproducible data revealed that 18,749 mRNAs and 6,061 lncRNAs were expressed during cardiac differentiation (Supplementary Figure S1A). As stage-specific marker genes, Nanog and Mesp1 were highly expressed at ESC and MES stages, while Mef2c and Tnnt2 were expressed at higher levels in the CP and CM phases (Figure 1A), indicating the high effectiveness of lncRNA-seq datasets. Comparative analysis of the expression patterns of mRNAs and lncRNAs showed that the expression of lncRNAs was lower than that of mRNAs. Over 50% of mRNAs had TPM values > 20 at four distinct stages, whereas only less than 10% of lncRNAs had similar TPM values (Supplementary Figure S1B). Mapping of these lncRNAs to genomic regions revealed that about half of the lncRNAs were derived from intergenic regions (Figure 1B).

To detect expression changes of lncRNAs during ESC-derived cardiac differentiation, we performed Fuzzy C-means Clustering based on the TPM values. This analysis led to the identification of six different gene clusters with differential lncRNA expression patterns across the differentiation phases (Figures 1C,D). Cluster one and six contained a batch of lncRNAs that were specifically expressed at MES and CP stages, respectively. Both cluster two and three comprised genes that were specifically expressed at the CM stage. However, lncRNAs in cluster two showed mildly increased expression from the ESC to CP stage but dramatically upregulated expression from the CP to CM stage. In contrast, the expressions of lncRNAs in cluster three were increased gradually with the progression of cardiac differentiation (Figures 1C,D). LncRNAs in cluster four showed decreased expressions from the ESC to CP stage and slightly recovered expression at the CM stage, while cluster five contained largest population of lncRNAs (n = 385) which expressed at highest levels in ESC and considerably downregulated afterwards. Altogether, our results suggested lncRNA transcriptome underwent dynamic regulation during ESC-derived cardiac differentiation and identified a batch of potential non-coding regulators involved in cardiogenesis.

It has been well documented that pluripotent genes such as Nanog and Sox2, are restricted to pluripotent cells and downregulated upon differentiation (Zhao and Jin, 2017). Similar as the regulatory pattern of pluripotent genes, we identified that 385 lncRNAs were expressed specifically in ESC and downregulated upon differentiation. Close examination of these lncRNAs revealed that the majority of them were less than 4,000 nucleotides in length (Figure 2A). Furthermore, 85.2% of ESC-specific lncRNAs harbored poly(A) tails, which might be associated with their stabilization in the cellular compartment. As canonical non-coding regulators of cell fate in the nucleus, Xist, Lncenc1, and Gas5, were included in the subgroup of polyadenylated lncRNAs specifically expressed in ESC (Figure 2B), further suggesting the accuracy and validity of our lncRNA-seq data. Intriguingly, this cluster also contained a large number of members from two lncRNA families, Platr and Snhg, some of which have been reported to regulate the stem cell fate transition (Bergmann et al., 2015; Lu et al., 2019; Li et al., 2020; Du et al., 2021; Yan et al., 2021).

To understand the key regulators that control the gene expressions of these ESC-specific lncRNAs, we next assessed the TF-binding cis-elements imbedded in the promoter region of their coding genes. We observed that about 72.73% of the lncRNAs could be bounded by Zfp281 (Figure 2C), which has been reported to be not required for maintenance of ESCs but indispensable for proper differentiation by targeting core pluripotent genes (Fidalgo et al., 2011). In addition, Smad2/3 and Zic2 have been implicated in the self-renewal and ESC specification (Sakaki-Yumoto et al., 2013; Luo et al., 2015), and our analysis suggested that approximately 49.35 and 73.25% of ESC-specific lncRNAs were enriched for the Smad2/3 and Zic2 motif, respectively (Figure 2C).

The clustering analysis also led to identification of ten distinct expression patterns of protein-coding genes across four differentiation phases (Supplementary Figure S2). We then performed gene ontology (GO) enrichment analysis of ESC-specific protein-coding genes in cluster two that were co-expressed with specifically expressed lncRNAs in ESC. Our analysis revealed that they were significantly enriched in biological processes including meiotic cell cycle, protein-DNA complex assembly, ribosome assembly and protein translation. Interestingly, these co-expressed lncRNAs were also involved in RNA methylation and histone modification (Figure 2D). Consistent with their potential biological processes, these ESC-specific lncRNAs were co-expressed with genes encoding proteins with translation regulator activity, methyltransferase activity and histone acetyltransferase activity (Figure 2E). These results suggest that the ESC-retained lncRNAs are forming regulatory networks with stage-specific TFs and protein-coding genes to modulate ESC fate determination. Supporting this notion, a human ESC-specific lncRNA, lncPRESS1, has been reported to be not only directly regulated by transcription factor P53, but also interact with SIRT6 and modulate histone H3K56 and H3K9 acetylation to effect ESC state (Jain et al., 2016).

LncRNA generally modulates target gene expression through recruiting specific functional proteins (Statello et al., 2021). We therefore searched for ESC-specific lncRNA-binding proteins by integrative analysis with two RNA-binding protein databases from both NPinter and RNAinter (Teng et al., 2020; Kang et al., 2022). Eventually we identified 1719 pairs of lncRNA-protein interactions. GO analysis of these binding proteins revealed that they harbored molecular activities of single- and double-stranded RNA binding, RNA helicase and chromatin binding. Surprisingly, m⁶A-containing RNA binding activity was also found in these lncRNA-interacting proteins (Figure 3A).

For individual lncRNAs, Gas5 and Sox1/2ot, which are linked to the regulation of pluripotency and specific linage differentiation (Knauss et al., 2018; Tu et al., 2018; Xi et al., 2022), were able to recruit many functional factors, such as Ago and Srsf family proteins (Figure 3B). Consistent with multi-faceted roles reported previously (Plath et al., 2003; Fang et al., 2004; Patil et al., 2016), Xist lncRNA could guide the regulation of gene expression through coupling with a variety of proteins including the key components of PRC and m⁶A machineries (Figure 3B). Among the identified proteins, the most representative ones are the core pluripotency component Oct4, histone modification activator PRC2 and m⁶A machinery proteins such as Igf2bp1/2/3 and Ythdc1/2, which modulate lncRNA functions potentially through direct interaction (Figure 3C). These results suggest that these ESC-specific lncRNAs mediate multi-layered regulation of gene expression through transcription, epigenetic and epitranscriptomic modifications.

Given that the ESC-specific lncRNAs are directly recognized by m⁶A machinery proteins, we therefore hypothesized that at least some of them could be modified by RNA methylation. Indeed, 97 out of 385 lncRNAs contained typical m⁶A peaks. Distinct from m⁶A distribution in the vicinity of the stop codons of mRNAs (Dominissini et al., 2012; Meyer et al., 2012), m⁶As in lncRNAs were highly enriched in close proximity to the middle of the transcript body (Figures 4A,C). In addition, over 85% of mRNA m⁶As were subject to the regulation of core methyltransferase Mettl3. While, only less than 13% of m⁶As in ESC-associated lncRNAs were regulated by Mettl3 (Figure 4B), suggesting the existence of additional methyltransferase catalyzing the m⁶A formation of these lncRNAs.

Since Ythdc1 has recently been identified to regulate the functionality of a number of non-coding RNAs (Wei and He, 2021), we next analyzed whether these m⁶A-containing lncRNAs could recruit Ythdc1 as the “reader” protein. As shown in Figure 4C, 28 out of 97 lncRNAs harbored m⁶A peaks that were targeted for regulation by Ythdc1 (Figure 4C). The remaining 69 lncRNAs may be regulated by other “reader” proteins. Surprisingly, although methylated adenosines of the 28 lncRNAs were directly bound by Ythdc1, almost all of them did not significantly change their expressions upon loss of Mettl3 or Ythdc1 (Figures 4D–F). A recent study has indicated that the RNA transcribed from long interspersed nuclear element-1 (LINE1) acts as a binding scaffold for Ythdc1 to regulate transcription activities in early embryo and mouse ESCs (Chen et al., 2021; Liu et al., 2021). It is possible that Ythdc1 may facilitate the formation of lncRNA-scaffold complex rather than regulate its abundance to modulate the nuclear gene expression in ESC. Intriguingly, one out of the 28 lncRNAs, Gm2379, was dually regulated by Mettl3 and Ythdc1, as loss of either m⁶A machinery protein resulted in obviously increased expression of Gm2379 (Figures 4E,F), suggesting that m⁶A may form different regulatory mode with lncRNAs to modulate their functionality.

We next focused on the m⁶A-containing transcript Gm2379, as it is the only lncRNA that responds to the loss of the m⁶A machinery proteins. We initially validated the expression pattern of Gm2379 during ESC-derived cardiac differentiation. As shown in Figure 5A, its expression was dramatically decreased at day 2 when embryonic bodies (EBs) were formed. Comparatively, the abundance of core pluripotency gene Nanog still maintained at high levels at this stage and declined after EBs generation. Nkx2-5 as one of the earliest markers of the cardiac differentiation, was found to be obviously upregulated at day 6. To decipher the function of Gm2379, we next probed its expression in separated nuclear and cytoplasmic fractions from cultured ESCs. Different from the nuclear 45S rRNA, Gm2379 was localized in both compartments (Figure 5B). Consistent with the global pattern of 97 ESC-specific lncRNAs, the two main m⁶A peaks of Gm2379 were distributed in proximity to the middle of the transcript body as well as subject to the regulation of Mettl3 and Ythdc1, as m⁶A signals of Gm2379 were overlapped with Ythdc1 RIP-seq signals and almost disappeared upon knockout of Mettl3 (Figure 5C). In addition, the transcript levels of Gm2379 were significantly increased in ESCs depleting Mettl3 or Ythdc1 (Figures 5D,E), which is in good agreement with previous data analysis. Altogether, our findings indicate that Gm2379 is specifically expressed in ESCs and dynamically regulated by nuclear “reader” protein Ythdc1 in a m⁶A-dependent manner.

m⁶A has been implicated in the regulation of ESC self-renewal and differentiation. Depletion of Mettl3 results in m⁶A erasure on target pluripotent mRNAs, facilitating ESC self-renewal and impairing directed differentiation (Batista et al., 2014; Geula et al., 2015). The nuclear “reader” Ythdc1 has been found to be required for the maintenance of ESC pluripotency in a m⁶A-dependent manner (Chen et al., 2021; Liu et al., 2021). We have revealed that the expression of Gm2379 lncRNA is regulated by both Mettl3 and Ythdc1, serving as a potential target of m⁶A in ESC. To determine whether Gm2379 is directly involved in the regulation of pluripotency, we generated stable ESC lines with shRNA-mediated gene knockdown of Gm2379 (Figure 6A). LIF is indispensable to maintain the pluripotent state of ESC (Martello and Smith, 2014). In the presence of LIF, depletion of Gm2379 selectively affected the expression of pluripotency genes. The transcript levels of Nanog and Dppa3 were significantly decreased in Gm2379-knock down (KD) cells, while Klf2 was upregulated in Gm2379-KD cells (Figure 6B). However, in the absence of LIF, the expression of Nanog and Klf2 rather than Dppa3 was comparable in control and Gm2379-KD ESCs, suggesting the selective regulation of pluripotency genes by Gm2379 lncRNA in a LIF-dependent manner.

In addition to the pluripotent genes, the crucial genes involved in the formation of three germ layers (ectoderm, mesoderm and endoderm) were dramatically dysregulated in Gm2379-KD ESCs (Figure 6C). Specifically, the ectoderm gene Fgf5 and endodermal gene Gata4 were significantly downregulated while the transcripts of the mesoderm gene T and ectoderm gene Otx2 were obviously increased in response to Gm2379 depletion. Mettl3, as the core m⁶A methyltransferase, has been demonstrated to regulate ESC pluripotency and differentiation (Batista et al., 2014; Geula et al., 2015). In contrast to Gm2379 KD, loss of Mettl3 resulted in increased expression of Fgf5 and decreased expression of Otx2 (Figure 6C). Considering the negative regulation of Gm2379 by Mettl3 shown above, it is very possible that Mettl3 directly targets Gm2379 lncRNA to modulate the expression of ectodermal genes.