HeLa cells were purchased from the American Type Culture Collection (ATCC, http://www.atcc.org) and authenticated by short-tandem-repeat (STR) profiling. They were cultured under standard conditions with DMEM (Gibco, 11995065) containing 10% FBS (CLARK, FB25015), and 1% penicillin/streptomycin (Beyotime, C0222) at 37°C with 5% CO₂. HeLa cells were determined with a PCR-based method and DAPI staining to ensure no contamination of mycoplasma.

To induce ER stress, HeLa cells were cultured with the DMEM medium containing 300 nM TG (Sigma, T9033) for 6 h. For glucose deprivation, cells were cultured with DMEM without glucose (Gibco, 11966025) at 37°C for 6 h.

Total RNA was extracted by TRIzol reagent (Invitrogen, 15596026) according to the manufacturer’s instructions. The concentration and quality of extracted RNAs were verified by Nanodrop and gel electrophoresis, respectively. Libraries were constructed by the TruSeq Ribo Profile Library Prep Kit (Illumina, RPHMR12126) according to the manufacturer’s instructions. In brief, 10 μg total RNA was depleted rRNA with the Illumina Rio-Zero Gold Kit (Illumina, MRZE724) and next purified for end repair and 5′ adaptor ligation. Then, the reverse transcription was performed with random primers containing the 3′ adaptor and randomized hexamer sequences. Finally, complementary DNA (cDNA) was purified and amplified with a Thermal Cycler. The PCR products of 300–500 base pairs (bp) were purified, quantified and stored at −80°C before sequencing. The libraries were subjected to 150-nt paired-end sequencing generating a depth of 50–100 million read pairs with an Illumina Novaseq platform (Novogene).

For data processing, the adaptors were trimmed with Cutadapt to obtain clean reads. The data quality was then checked with FastQC and the low-quality (Q value ≤ 20) reads were removed. The remaining reads were subsequently aligned to the human reference genome (hg19) with Bowtie2 (-v 1). For linear RNAs including lncRNAs and mRNAs, the corresponding reads were counted with BEDtools and read per million (RPM) was used to calculate levels for lncRNAs and mRNAs. LncRNAs with an average RPM ≥0.1 were used for further analysis and the DE lncRNAs were determined by DEseq2 with a criterion of fold change ≥2 or ≤0.5 and p-value < 0.05. The DE mRNAs were determined by DEseq2 with the cutoff (the average RPM ≥10, fold change ≥2 or ≤0.5 and p-value < 0.05). For circRNA prediction, find_circ and CIRI2 were applied to identify high-confidence BSJs with default parameters. Only circRNA candidates predicted by both pipelines were used for further investigations and CIRI2-annotated BSJ reads (BSJ reads ≥2) were used to calculate circRNA levels. The DE circRNAs were determined by DEseq2 with the cutoff (fold change ≥2 or ≤0.5, and p-value < 0.05).

cDNA was synthesized from 500 ng total RNA with the GoScript Reverse Transcription System (Promega, A5000) according to the manufacturer’s protocol. For RT-PCR gels of XBP1 and FST mRNAs, amplification was performed with 30 cycles. Real-time quantitative PCR (RT-qPCR) was carried out with GoTaq SYBR Green qPCR Master Mix (Promega, A6001) on a QuantStudio Applied Biosystems (Thermo) according to standard procedures. All amplification curves reached the stationary stage before 35 cycles and the readings of the Ct value were obtained at the exponential stage. ACTB mRNA was used as an internal control. All PCR products were sequenced for confirmation and all primer sequences were included in Supplementary Table S3.

All plasmids for the luciferase reporter system were constructed with recombinant methods (Vazyme, c113-02). PRDM1 3′ UTR containing the binding sites of miR-9-5p and SOX12 3′ UTR containing the miR-744-5p binding sequences were PCR-amplified from the cDNA of HeLa cells and then inserted into the Firefly luciferase reporter vector pGL3-control (Promega, E1741) between XbaI (Thermo, FD0685) and FseI (NEB, R0588V) double-digested sites. Small interfering RNAs (siRNAs) targeting lncSLC25A1, TINCR and circBANP BSJ were synthesized by GenePharma (Shanghai, China). Transfection of plasmids and siRNAs was performed with Lipofectamine 2000 (Invitrogen, 11668019) according to the manufacturer’s protocol. Oligonucleotide sequences for primers used in plasmid construction and siRNAs are included in Supplementary Table S3.

The lncRNA/circRNA-mRNA co-expression network was constructed according to the expression levels in our dataset. Briefly, 51 lncRNAs, 39 circRNAs, and 279 mRNAs sensitive to TG and GD treatments were used to construct the network with a criterion of Spearman R ≥ 0.95 or ≤ -0.95 and p-value < 0.01. The constructed network consisted of 29 lncRNAs, 20 circRNAs, and 131 mRNAs, and was visualized with Cytoscape (https://cytoscape.org/).

The ceRNA network was constructed based on the lncRNA/circRNA-miRNA-mRNA axis. The lncRNA-miRNA, circRNA-miRNA, and mRNA-miRNA interactions were predicted with TargetScanHuman (https://www.targetscan.org/vert_72/). Briefly, mRNAs, lncRNAs, or circRNAs with at least two putative binding sites for the individual miRNA were used to construct the network.

Briefly, 1×10⁶ HeLa cells were co-transfected with 30 pmol siRNAs, 1 μg pGL3 Firefly luciferase plasmids and 100 ng pRL Renilla luciferase reporter vector (Promega, E2261). After transfection for 48 h, cells were lysed with passive lysis buffer on ice for 20 min and the luciferase activity was performed with Dual Luciferase Reporter Assay System Kit (Promega, E1910) according to the manufacturer’s protocol. The Firefly luciferase activities were measured and normalized to Renilla luciferase activities (F/R).

The RNA-RBP network was constructed based on the lncRNA-RBP and circRNA-RBP interactions predicted by RBPmap (http://rbpmap.technion.ac.il/). The lncRNA/circRNA-RBP interaction was determined by a criterion of the individual RNA sequence containing more than 2 motifs for RBP of interest (p-value < 10E-4). The RNA-RBP interactome map consisted of 46 lncRNAs, 17 circRNAs and 77 RBPs, and was visualized with Cytoscape.

FISH was carried out as previously described with minor modifications (Hu et al., 2016). RNA probe antisense to linc00612 was generated by the Transcript Aid T7 High Yield Transcription Kit (Thermo, K0441) with the corresponding insertion into the T vector (Promega, A3600) as a template. The probe was then labeled with Alexa Fluor546, by using the ULYSIS Nucleic Acid Labeling Kit (Thermo, U21652), which added a fluor on every G of the probe to amplify the fluorescence intensity. The primers for the antisense probe amplification were included in Supplementary Table S3. HeLa cells were fixed with 4% PFA for 10 min at room temperature after washing with PBS twice. RNA probes were denatured at 80°C for 10 min and placed on ice immediately. Fixed cells were incubated with RNA probes mixed with 20 ng/μL human Cot-1 DNA (Invitrogen, 15279011), 500 ng/μL yeast total RNA (Invitrogen, AM7118) and 10 units/mL RNase inhibitor (Promega, N2615) in 2 × hybridization buffer (4 × SSC, 40% dextran sulfate) at 37°C for 15–17 h, protected from light. After two 10-min washes in SSCT (2 × SSC and 0.4% Tween 20) buffer, nuclei were stained with DAPI (Sigma, F6057). Finally, images were captured using the LSM 980 confocal microscope (Zeiss).

Published eCLIP-seq data were obtained from Gene Expression Omnibus (GEO) database under the following accession numbers: GSE91952 (EIF4G2), GSE126263 (MSI1), GSE71096 (SRSF10), GSE69153 (RC3H1) and GSE107768 (FUBP3, DAZAP1, HNRNPA0 and PABPC4).

GO analysis of the DE mRNAs was performed using the GOrilla web-server with default parameters (http://cbl-gorilla.cs.technion.ac.il) (Eden et al., 2009). For data visualization, the plots were generated by the ggplot2 package in R software.

The physiological experiments were carried out in triplicates (N = 3), and statistical analysis of the data was performed with the two-tailed Student’s t-tests. Data were present as the mean from three independent experiments with SEM. The RNA-seq was performed with four replicates and statistical analysis for analysis was calculated by DEseq2.

To identify stress-related ncRNAs, we performed ribosomal RNA (rRNA) depleted RNA sequencing (RNA-seq) of human HeLa cells treated with TG (an ER stress inducer) or GD (a metabolic oxidative stress inducer). XBP1s is a well-characterized marker for ER stress and Follistatin (FST) is upregulated in response to GD (Yoshida et al., 2001; Gao et al., 2014). We first examined the levels of XBP1s and FST mRNAs, and found that XBP1s and FST detected by RT-PCR with specific primers were significantly increased upon TG and GD treatment, respectively (Supplementary Figure S1A). Then, we performed the bioinformatics to analyze lncRNAs and circRNAs in our dataset (Supplementary Figure S1B). To identify high-confidence back-splicing junctions (BSJs) of circRNAs, two published pipelines including CIRI2 (Gao et al., 2018) and find_circ (Memczak et al., 2013) were applied to annotate circRNAs (Supplementary Figure S1B). CircRNA candidates overlapped in these two pipelines were used for further investigations. Then, the differentially expressed (DE) circRNAs and lncRNAs were analyzed (Supplementary Figure S1B). Principal component analysis (PCA) plots with all identified transcripts also revealed that four biological replicates clustered together, and the controls, TG and GD groups were clearly separated (Supplementary Figure S1C), indicating reliability of RNA-seq and data analyses.

We characterized 2,406 lncRNAs in our dataset, among which 179 (97 upregulated, 82 downregulated) and 193 (82 upregulated, 111 downregulated) DE lncRNAs were identified (fold change ≥2 or ≤0.5, p-value < 0.05) upon TG and GD treatments, respectively (Figure 1A and Supplementary Table S1). Among them, 26 and 25 lncRNAs were significantly increased and decreased in response to both TG and GD treatments (Figure 1B). Simultaneously, we discovered ∼1‰ BSJ reads in our dataset and identified a total of 6531 circRNAs (BSJ reads ≥2) (Supplementary Figure S1D). We determined 161 (70 upregulated, 91 downregulated) and 147 (51 upregulated, 96 downregulated) DE circRNAs (fold change ≥2 or ≤0.5, p-value < 0.05) upon TG and GD treatments, respectively (Figure 1C and Supplementary Table S2). Among them, 14 and 25 circRNAs were markedly increased and decreased in response to both TG and GD treatments (Figure 1D). Then, we focused on 51 lncRNAs and 39 circRNAs both dysregulated in response to ER and metabolic stresses (Figures 1B, D). We randomly selected 6 lncRNAs and 6 circRNAs among them for experimental validation. Real-time quantitative PCR (RT-qPCR) analysis demonstrated that lincPINT, lncSARS1, lncZFAS1, circCRIM1, circHBS1L, and circPTBP2 were significantly increased after TG and GD treatments (Figures 2A, B). Conversely, BASP1-AS1, lncNT5C3B, TINCR, circABL2, circBANP, and circCEP72 were significantly decreased upon TG and GD treatments (Figures 2A, B). For these candidates examined, RT-qPCR-mediated verification was highly consistent with the RNA-seq data (Figure 2C). Taken together, the DE lncRNAs and circRNAs might play essential roles in coping with cellular stresses.

The functions of lncRNAs and circRNAs are tightly related to the roles of their co-expressed protein-coding genes (Quinn and Chang, 2016; Sheng et al., 2019; Chen et al., 2022). To investigate the roles of lncRNAs and circRNAs associated with ER and metabolic stresses, we constructed a co-expression network (Spearman R ≥ 0.95 or ≤ -0.95, p-value < 0.01) of lncRNAs, circRNAs and their co-expressed DE mRNAs (RPM ≥10, fold change ≥2 or ≤0.5, p-value < 0.05) for the 51 DE lncRNAs and 39 DE circRNAs identified in both stresses according to our RNA-seq data (Figure 3A). Our analysis revealed that 29 lncRNAs (17 downregulated, 12 upregulated) and 20 circRNAs (18 downregulated, 2 upregulated) interacted with 131 DE mRNAs (72 downregulated, 59 upregulated) (Figure 3A). Notably, we found that the genomic distances of all interacted nodes in the co-expressed network were more than 100 kilobases, indicating the trans roles of lncRNAs and circRNAs on the co-expressed mRNAs. RNA interference (RNAi) is a widely used approach to deplete lncRNA/circRNA of interest, although there are multiple other methods (Quinn and Chang, 2016; Zhang et al., 2021). We then performed loss of function studies of three candidates (two lncRNAs and one circRNA) using RNAi to validate the co-expressed correlation between DE lncRNAs/circRNAs and mRNAs (Figures 3B–D). Small interfering RNA (siRNA)-mediated TINCR silencing resulted in significantly decreased expressions of its co-expressed mRNAs (ERN1 and SOX12) in HeLa cells (Figure 3B). Consistently, lncSLC25A1 knockdown decreased its co-expressed targets (SOX12 and ATOH8), and circBANP depletion downregulated its co-expressed targets (ERN1 and ATF4) (Figures 3C, D). The heatmap further demonstrated that these 131 mRNAs were significantly dysregulated after TG and GD treatments (Supplementary Figure S2A). Gene Ontology (GO) analysis revealed that these mRNAs present in the co-expressed network significantly enriched in biological processes such as response to toxic substance, response to stress, response to extracellular stimulus, regulation of RNA metabolic process, programmed cell death, PERK-mediated unfolded protein response and cellular response to glucose starvation (Supplementary Figure S2B).

One of the molecular mechanisms for lncRNAs and circRNAs is to act as competing endogenous RNAs (ceRNAs), in which lncRNAs and circRNAs bind with miRNAs and decrease the corresponding binding of miRNAs on mRNA targets (Quinn and Chang, 2016; Chen et al., 2022). In order to explore the potential lncRNAs and circRNAs serving as ceRNAs responsive to cellular stresses via sequestering miRNAs and thus regulating mRNA targets, we constructed a ceRNA network among the DE lncRNAs, circRNAs, and mRNAs respond to both ER and metabolic stresses. 23 circRNAs, 2 lncRNAs, 48 miRNAs, and 32 mRNAs composed of the ceRNA regulatory network (Figure 4A). Among 32 mRNA targets, 17 were downregulated and 15 were upregulated in response to TG and GD treatments (Figure 4B). Then, we conducted GO analysis for the 32 mRNAs in the ceRNA network and found that they were associated with regulations of biological processes including transcription by RNA polymerase II, RNA metabolic process, macromolecule metabolic process, developmental process and cell differentiation (Figure 4C). To further verify the ceRNA roles of lncRNAs and circRNAs, we performed a dual-luciferase reporter assay in HeLa cells. In our network, circBANP interacted with miR-9-5p to derepress the PRDM1 level (Figure 4D). Knockdown of circBANP mediated by siRNA significantly reduced the luciferase activity of Firefly with the PRDM1 3’ UTR region containing the miR-9-5p binding site (Figure 4D). The lncSLC25A1-miR-744-5p-SOX12 axis exhibited a similar result according to our experiments (Figure 4D). Although, we have to point out that these were just results from overexpression and RNAi experiments, and further validations are required to examine the ceRNA regulatory network.

Given that lncRNAs and circRNAs could interact with RBPs to exert their crucial functions (Chen et al., 2022; Wang et al., 2022), we predicted the potential binding RBPs of both DE lncRNAs and DE circRNAs in ER or metabolic stresses via RBPmap, a tool mapping the interacting RBPs for RNAs of interest (Paz et al., 2014), and constructed the lncRNA/circRNA-RBP network based on the predicted interactions (Figure 5A). The network consisted of 17 circRNAs, 46 lncRNAs and 77 RBPs (Figure 5A). We also noticed that among the network, linc00612 and circSTAU2 interacted with the most RBPs among lncRNAs and circRNAs, respectively (Figure 5A). Linc00612 interacted with 16 different RBPs such as SRSF8, HNRNPDL, PUM2, and circSTAU2 interacted with 6 distinct RBPs such as PABPC1, SART3, SRSF10 (Figures 5B, C). Furthermore, the mRNA levels of 77 RBPs were dynamic after TG and GD treatments (Supplementary Figure S3A). GO analysis for these 77 RBPs indicated that they were related to biological processes including RNA 3’ end processing, mRNA stability involved in response to stress, regulation of translation, gene silencing by miRNA, nuclear export, mRNA splicing via spliceosome, etc. (Supplementary Figure S3B). RNA fluorescence in situ hybridization (FISH) displayed that linc00612 mostly localized to the cytoplasm, and a small portion resided in the nucleus in HeLa cells (Figure 5D), implying that linc00612 might possess both cytoplasmic and nuclear roles. In addition, integrative analysis of published eCLIP-seq data of RBPs revealed that EIF4G2, FUBP3, HNRNPA0, and MSI1 demonstrated binding signals on linc00612 (Figure 5E).