Exon–Intron Differential Analysis Reveals the Role of Competing Endogenous RNAs in Post-Transcriptional Regulation of Translation

Stressful conditions induce the cell to save energy and activate a rescue program modulated by mammalian target of rapamycin (mTOR). Along with transcriptional and translational regulation, the cell relies also on post-transcriptional modulation to quickly adapt the translation of essential proteins. MicroRNAs play an important role in the regulation of protein translation, and their availability is tightly regulated by RNA competing mechanisms often mediated by long noncoding RNAs (lncRNAs). In our paper, we simulated the response to growth adverse condition by bimiralisib, a dual PI3K/mTOR inhibitor, in diffuse large B cell lymphoma cell lines, and we studied post-transcriptional regulation by the differential analysis of exonic and intronic RNA expression. In particular, we observed the upregulation of a lncRNA, lncTNK2-2:1, which correlated with the stabilization of transcripts involved in the regulation of translation and DNA damage after bimiralisib treatment. We identified miR-21-3p as miRNA likely sponged by lncTNK2-2:1, with consequent stabilization of the mRNA of p53, which is a master regulator of cell growth in response to DNA damage.


Introduction
The transcriptome is tightly regulated at different levels: along with the regulation of new transcription, RNA molecules can be modulated at the post-transcriptional level, and noncoding RNAs play a relevant role in this process [1,2]. The most famous mechanism of post-transcriptional regulation is mediated by microRNAs, which are small noncoding RNAs that inhibit translation of mRNAs or induce their degradation [3]. In addition, long noncoding RNAs (lncRNAs) take part in this mechanism. Some of them possess miRNA responsive elements (MREs) [4] and are known as competing endogenous RNAs (ceRNAs). Indeed, they can act as an endogenous miRNA decoy, and their expression modulates the amount of free miRNAs available to bind their targets [5]. Circular RNAs (circRNAs) constitute a subclass of exceptionally stable ceRNA molecules that contain a covalent circular structure formed by noncanonical 3 to 5 end-joining event called back-splicing. circRNAs are diffuse and sometimes conserved in eukaryotic organisms [6]. High throughput analysis of protein expression, as liquid chromatography mass spectrometry, detects changes in translation efficiency, but it is expensive and requires specific expertise. Here, we show that also transcriptome profiling can help us to identify post-transcriptional modulations. RNA-Seq can be used to investigate pre-mRNA dynamics, taking advantage of the intronic sequences that are also acquired although less abundantly than exonic sequences. Changes in the overall intronic read counts directly measure changes in transcriptional activity. We can discriminate whether RNA levels are modulated at the transcriptional and post-transcriptional level applying an exon-intron split analysis (EISA), which compares intronic and exonic changes across different experimental conditions [7]. Here, we applied this approach to determine whether post-transcriptional regulation mediated by lncRNAs might be an additional layer to quickly control protein expression in diffuse large B cell lymphoma (DLBCL) cells exposed to bimiralisib, a dual PI3K/mTOR inhibitor with proven preclinical and early clinical anti-lymphoma activity [8,9]. The mTOR pathway regulates cell growth and proliferation in response to mitogen, nutrient, and energy status and therefore controls the balance between anabolism and catabolism in response to environmental conditions [10][11][12]. In addition to various anabolic processes as protein, lipid and nucleotide synthesis, mTOR also promotes cell growth by suppressing protein catabolism. mTOR is a downstream mediator of several growth factor and mitogen-dependent signaling pathways but also reacts to intracellular stresses that are incompatible with growth such as low ATP levels, hypoxia, or DNA damage.
The exposure of DLBCL cell lines to bimiralisib induces important post-translational and transcriptional changes, affecting genes and proteins involved in the PI3K/AKT/mTOR pathway, BCR signaling, NF-kB pathway, mRNA processing, apoptosis, cell cycle, Myc pathway, MAPK/RAS signaling, and glycolysis [8]. Most of the pathways are largely modulated via both transcription regulation and protein phosphorylation changes [8]. We hypothesized that post-transcriptional modulation can play a role in the mechanism of action of bimiralisib also. A fast adaptation to lack of nutrients requires an optimization of the stability of transcripts that need to be translated into proteins. Here, we describe how bimiralisib induced a quick stabilization of transcripts needed to cope with amino acid deficiency and to modulate the translation. This event strongly correlated with the overexpression of a lncRNA, namely lncTNK2-2:1, which is associated to the increased stability of transcripts affected by mTOR inhibition and responsible for DNA damage response. In particular, we validated the stabilization of p53 transcript due to the sponge effect on miR21-3p mediated by lnc-TNK2-2:1.

Bimiralisib Reduces the Transcription of Genes Coding for Proteasome and Ribosome Components
We have previously reported that dual PIK3/mTOR pharmacological inhibition using bimiralisib has in vitro and in vivo anti-tumor activity and that it induces transcriptional changes in DLBCL cell lines [8]. Here, we have performed total RNA-Seq on RNA extracted from two DLBCL cell lines, U2932 and TMD8, exposed to DMSO or to bimiralisib for 4, 8, or 12 h. We noticed a general reduction for the most part of transcripts encoding for subunits of RNA pol I and III, which are responsible of rRNA transcription and of some subunits of RNA pol II, and of many transcripts encoding for the machinery processing rRNA. PIK3/mTOR inhibition induced also a downregulation of proteasome components, which was probably to balance the reduction of protein synthesis due to impaired ribogenesis (Supplementary Figure S1 and Table S1).

Post-Transcriptional Regulation of Many Transcripts Encoding for Riboproteins and Translation Regulators Is an Early Event after Dual PIK3/mTOR Inhibition
We applied EISA [7] to transcriptomic changes upon bimiralisib treatment. EISA measures changes in mature RNA and pre-mRNA reads across different conditions to quantify the transcriptional and post-transcriptional regulation of gene expression. After 4 h of PIK3/mTOR inhibition, we observed a marked post-transcriptional regulation, since many transcripts showed independent changes in exons and introns. After 8 h, changes in transcript levels between DMSO or bimiralisib-treated samples have been mainly due to alteration of transcription, as evident by the fact that changes in exons were well correlated with change in introns (Figure 1a). Focusing on transcripts early stabilized after PIK3/mTOR inhibition, we found they were mainly encoding ribosome components involved in the response to amino acids starvation. Among the most quickly degraded transcripts, we found mRNAs encoding for spliceosome components and pre-mRNA processing (Figure 1b and Table S2). many transcripts showed independent changes in exons and introns. After 8 h, chan in transcript levels between DMSO or bimiralisib-treated samples have been mainly to alteration of transcription, as evident by the fact that changes in exons were well co lated with change in introns (Figure 1a). Focusing on transcripts early stabilized a PIK3/mTOR inhibition, we found they were mainly encoding ribosome component volved in the response to amino acids starvation. Among the most quickly degraded t scripts, we found mRNAs encoding for spliceosome components and pre-mRNA cessing ( Figure 1b and Table S2).  exonic (∆exon) reads per each expressed transcript in DLBCL cell lines exposed to bimiralisib or DMSO for 4, 8, and 12 h. ∆exon/∆intron = 1 reflects transcriptional modulation, ∆exon/∆intron = 1 reflects post-transcriptional modulation. (b) Representative GSEA plots illustrating the transcriptional expression signature enrichment in transcripts ranked by their decreasing ∆exon/∆intron ratio in DLBCL cell lines exposed to bimiralisib (1 µM) or DMSO for 4 h. Green line, enrichment score; bars in the middle portion of the plots show where the members of the gene set appear in the ranked list of genes. Positive or negative ranking metrics indicate correlation or inverse correlation with the profile, respectively. FDR, false discovery rate; NES, normalized enrichment score.
2.3. The lincRNAs RP11-480A16.1 (lncTNK2-2:1) and GMDS-AS1 Are Differentially Expressed after Dual PIK3/mTOR Inhibition and Strongly Correlated to Significantly Stabilized Transcripts We hypothesized that a prompt alteration in the stability of several transcripts could be achieved by the expression of lncRNAs acting as miRNA sponges. We found 20 significant lncRNAs differentially expressed after pharmacological PIK3/mTOR inhibition: 15 upregulated and five downregulated. We also evaluated whether the 20 lncRNAs could be circRNAs, applying the algorithm CiriQuant [13] that accurately determines the junction of circRNAs from RNA-Seq paired samples. The output showed 28,521 back-spliced junctions (BSJ), most of which were cell specific (TMD8: 9519, U2932: 14,348). Six of the differentially expressed lincRNAs were among the reliably quantified circRNAs ( Figure 2a, Table 1). We calculated the Pearson correlation index of each selected lincRNAs with each stabilized or degraded transcript (∆ex/∆intr = 1). The expression of two of these lncRNAs, namely RP11-480A16.1 (lncTNK2-2:1) and GMDS-AS1, strongly correlated with transcripts differentially stabilized upon bimiralisib exposure ( Figure 2b and Table S3). The post-transcriptionally modified transcripts, ranked by their correlation index with lncTNK2-2:1 and GMDS-AS1, were enriched in genes involved in the regulation of translation and in the response to amino acid starvation ( Table 2). In particular, there was a highly significant enrichment of genes affected by mTOR inhibitor rapamycin, highlighting the prominent role of lncTNK2-2:1 and GMDS-AS1 in the post-transcriptional regulation due to bimiralisib exposure. We also noticed the enrichment of genes involved in DNA damage, and in particular, some of them were regulated by miRNAs as well. We focused on p53 and ATM due to their important biologic roles ( Figure 3 and Table S4).

lncTNK2-2:1 Induces Stabilization of p53 and ATM by Sequestering miR21-3p
We validated by qRT-PCR the correlation observed between the expression of lncTNK2-2:1 and GMDS-AS1 and the stability of ATM and p53 mRNAs. We measured the levels of the pre-mRNAs or the mature transcripts of ATM and p53 genes and the expression of lncTNK2-2:1 and GMDS-AS1, in U2932 and TMD8, at 4 h and 8 h after exposure to 1 µM of bimiralisib or DMSO. We confirmed the upregulation of both GMDS-AS1 and lncTNK2-2:1, although the GMDS-AS1 upregulation was significant in U2932 but not in TMD8 (Figure 4a). We measured the fold change of lncRNAs and of pre mRNAs or mature transcripts ( Figure S2a) and calculated their Pearson correlation index ( Figure  S2b). We clearly confirmed the predicted correlation of the expression of the lncRNA lncTNK2-2:1 with the mature transcripts for p53 and ATM (R = 0.581 and R = 0.596, respectively) (Figure 4b, bottom panels) but not with their pre-mRNAs (Figure 4b, top panels). GMDS-AS1 was significantly correlated neither to p53 nor to ATM stability ( Figure  4b). According to this, we searched for experimentally validated miRNA binding sites sheared by lncTNK2-2:1 and p53 in DIANA tools, LncBase (https://diana.e-  We validated by qRT-PCR the correlation observed between the expression of lncTNK2-2:1 and GMDS-AS1 and the stability of ATM and p53 mRNAs. We measured the levels of the pre-mRNAs or the mature transcripts of ATM and p53 genes and the expression of lncTNK2-2:1 and GMDS-AS1, in U2932 and TMD8, at 4 h and 8 h after exposure to 1 µM of bimiralisib or DMSO. We confirmed the upregulation of both GMDS-AS1 and lncTNK2-2:1, although the GMDS-AS1 upregulation was significant in U2932 but not in TMD8 (Figure 4a). We measured the fold change of lncRNAs and of pre mRNAs or mature transcripts ( Figure S2a) and calculated their Pearson correlation index ( Figure S2b). We clearly Non-coding RNA 2021, 7, 26 7 of 13 confirmed the predicted correlation of the expression of the lncRNA lncTNK2-2:1 with the mature transcripts for p53 and ATM (R = 0.581 and R = 0.596, respectively) (Figure 4b, bottom panels) but not with their pre-mRNAs (Figure 4b, top panels). GMDS-AS1 was significantly correlated neither to p53 nor to ATM stability (Figure 4b). According to this, we searched for experimentally validated miRNA binding sites sheared by lncTNK2-2:1 and p53 in DIANA tools, LncBase (https://diana.e-ce.uth.gr/lncbasev3 accessed 15 April 2021) [4] and TarBase (https://carolina.imis.athena-innovation.gr/diana_tools/web/index. php?r=tarbasev8%2Findex accessed 15 April 2021) [14] (Table S5). miR-21-3p and miR-22-3p are reported to target both lnc-TNK2-2:1 and p53 (Figure 4c), and miR-21-3p belongs to the miR-21-3p/TSC2/mTOR regulatory axis [15]. Thus, we focused on miR-21-3p, which was consistently reduced in both U2932 and TMD8 after bimiralisib treatment (Figure 4d) and the consequent upregulation of lncTNK2-2:1 (Figure 4a). The p53 stabilization was more robust in TMD8 than in U2932 ( Figure S2a), even if the lncRNA was upregulated in both cell lines (Figure 4a), which was an observation that was possibly justified by the lower basal expression of miR21-3p in U2932 compared to TMD8 ( Figure S2c).

lncTNK2-2:1 Degradation Reverts Stabilization of p53 and Releases miR21-3p
To validate the relationship between lncTNK2-2:1 expression and miR21-3p activity on its target p53, we electroporated 2 million TMD8 with 100 pmol of antisense oligonucleotides or of the negative control. After 72 h, we exposed the cells to 1 µM of bimiralisib or DMSO for 8 h. We measured the levels of the pre-mRNA or the mature transcript of the p53 gene and the expression of lncTNK2-2:1 by qRT-PCR. We confirmed that the antisense oligonucleotides degraded lncTNK2-2:1 efficiently both at basal condition and after its induction due to bimiralisib treatment (Figure 4e, left). We validated the relationship between the lncRNA lncTNK2-2:1 with the p53 transcript stability. Indeed, after lncTNK2:2-1 interference, we noticed the reduction of p53 stabilization (Figure 4e, middle), which was measured as the fold change between p53 mature mRNA with respect to the total pre-mRNA transcribed in each condition. According to this, we measured miR-21-3p after lncTNK2-2:1 degradation, and we showed that it increased in samples where lncTNK2-2:1 was knocked down with respect to the negative control (Figure 4e, right). This experimental observation confirmed the predicted miR-21-3p/lncTNK2-2:1/p53 regulatory axis.

Discussion
As a central controller of cell growth, mTOR regulates ribosome biogenesis. The latter is the most energy-demanding cellular process, and mTOR controls it by promoting the translation of riboproteins and by affecting ribosomal RNA (rRNA) synthesis. Ribosome synthesis requires all three nuclear RNA polymerases, Pol I for the synthesis of rRNA, Pol II for transcription of riboprotein genes, and Pol III for the synthesis of 5S RNA [16,17]. PIK3/mTOR inhibition by bimiralisib leads to the downregulation of all of them, a reduction of rRNA gene transcription, and, in addition, pre-rRNA processing impairment. Moreover, we describe here an additional mechanism that the cell activates when it must save energy: mTOR pathway inhibition leads to the stabilization of already present transcripts encoding for riboproteins. Since the new transcription of riboproteins genes is inhibited, the cell needs to save the already available riboprotein mRNA as long as possible in order to still translate essential proteins.
A growing number of miRNAs have been shown to control components or regulators of ribosome biogenesis [18]. In addition, lncRNAs have been increasingly found to play relevant roles in eukaryotic ribosome biogenesis that can be basally active or stress responsespecific [18]. These molecular mechanisms include protein binding, rDNA chromatin modifications, snoRNP formation, and transcript-specific translation modulations. Here, we report evidence of an additional example of lncRNAs involved in the regulation of ribogenesis. Two lncRNAs, lncTNK2-2:1 and GMDS-AS1, were modulated following PIK3/mTOR pathway inhibition, and their transcription changes are strongly correlated with the stabilization of transcripts encoding for many riboproteins. Other lncRNAs were significantly upregulated by bimiralisib treatment but were not correlated with transcripts early stabilized after mTOR inhibition. Thus, we postulated that these particular lncRNAs could be relevant players of mTOR-mediated modulation of translation in response to amino acid deficiency or other stressful events. One of the possible mechanisms that might mediate a quick transcript stabilization is the sequestration of miRNAs, and GMDS-AS1 is already known as an miRNA sponge in lung cancer [19]. As proof of principle, we looked for binding sites in the lincRNAs for miRNAs that could target stabilized transcripts correlated with the upregulation of lncTNK2-2:1 or GMDS-AS1. Both lncTNK2-2:1 and GMDS-AS1 were associated with the stabilization of genes involved in DNA damage response and regulated by miRNAs. Among them, p53 is also involved in ribosome biogenesis [20,21], and the p53 network is known to interact with several miRNAs [22]. We could validate the correlation between p53 mRNA stabilization and lncTNK2-2:1 level but not with GMDS-AS1.
We showed that miR-21-3p increased after lncTNK2-2:1 silencing, and concomitantly, p53 was not stabilized after bimiralisib treatment in the absence of lncTNK2-2:1. This evidence enforced our hypothesis that miR-21-3p was responsible for p53 mRNA degradation, which we formulated on the basis of in silico base pairing of lncTNK2-2:1 and miR-21-3p and on the high stability of predicted miR-21-3p binding in p53 3 UTR, which is compatible with mRNA degradation, instead of inhibition of translation [23]. Furthermore, in U2932, a cell line expressing low miR-21-3p levels, p53 was not strongly stabilized, despite the upregulation of lncTNK2-2:1. In support of our findings, miR21-3p is already known to modulate the mTOR pathway via TSC2 mRNA downregulation [15], and P53-dependent mTOR inhibition is mediated by TSC2 [24]. Here, we provide the evidence of an additional layer of regulation of p53-mTOR crosstalk through the rapid elimination of miR-21-3p and consequent stabilization of p53 and enhancement of TSC2 repressor activity of mTOR. The PI3K/mTOR pharmacological inhibition enforces the mTOR inhibition by a positive feedback loop mediated by the overexpression of lncTNK2-2:1.
We also selected ATM as potential interesting transcript, since it appeared related to lncTNK2-2:1 overexpression both in silico and in vitro, but we could not identify any miRNA that may mediate the connection between the lncRNAs and the mRNA. ATM mRNA stabilization after bimiralisib exposure might also be due to the downregula-tion of miRNAs directly regulated by PI3K/signaling mTOR and directly targeting ATM 3 UTR [25].
In conclusion, based on an alternative bioinformatic approach applied to RNA-Seq data, we selected candidate molecules that could be involved in a post-transcriptional mechanism of RNA competition, and we provided data suggesting a novel RNA network composed by lncRNAs, miRNAs, and mRNAs, which is affected by the dual PI3K/mTOR pharmacological inhibition in DLBCL cell lines.

Cell Culture and Bimiralisib Treatment
Established human DLBCL cell lines TMD8 and U2932 were grown as previously described [8]. Bimiralisib was kindly provided by PQR Therapeutics (Basel, Switzerland). TMD8 and U2932 were seeded 3 million cells/well in a non-tissue culture 6 well plate. Cells were treated for 4 and 8 h with 1 µM bimiralisib or 0.1% DMSO (Sigma, St. Louis, MO, USA), respectively. Treatment was stopped by washing the cell with RNAse-free PBS and followed by immediate RNA extraction.

RNA-Extraction
For each cell line and condition, cells were collected and resuspended in 1 mL of TRI Reagent (Sigma, St. Louis, MO, USA) for cell lysis, and extraction was performed. DNAse (Qiagen, Hilden, Germany) was added to the RNA samples and incubated for 15 min at room temperature. Total RNA was reprecipitated to remove salts and the enzyme.

Whole-Transcriptome Sequencing (RNA-Seq)
Two cell lines, U2932 and TMD8, were treated with 1 µM of bimiralisib or DMSO and RNA was extracted after 4, 8, or 12 h. Cells treatment and RNA extraction were described in [8]. Quality control for extract RNA was performed on the Agilent BioAnalyzer (Agilent Technologies, Santa Clara, CA, USA) using the RNA 6000 Nano kit (Agilent Technologies, Santa Clara, CA, USA), and concentration was determined by the Invitrogen Qubit (Thermo Fisher Scientific, Waltham, MA, USA) using the RNA BR reagents (Thermo Fisher Scientific). The TruSeq RNA Sample Prep Kit v2 for Illumina (Illumina, San Diego, CA, USA) was used for cDNA synthesis and the addition of barcode sequences. The sequencing of the libraries was performed via a paired end run on a NextSeq500 Illumina sequencer (Illumina). As an average, 25 million reads were collected per each sample.

Data Mining
We evaluated the RNA-seq reads quality with FastQC (v0.11.5), and we removed low-quality reads/bases and adaptor sequences using Trimmomatic (v0.35). The trimmedhigh-quality sequencing reads were aligned using STAR [26], which is a spliced read aligner that allows for sequencing reads to span multiple exons. On average, we were able to align 85% of the sequencing reads for each sample to the reference genome (HG38).
Then, the HTSeq-count software package [27] was used for the quantification of gene level expression. Differential expression analysis was performed on gene-level read count data using the 'limma' pipeline [28,29] We first subsetted the data to genes that had a counts-per-million value greater than one in 3 or more samples. The data were normalized per sample using the 'TMM' method from the edgeR package [30] and transformed to log2 counts-per-million using the edgeR function 'cpm'. Then, linear model analyses, with empirical-Bayes moderated estimates of standard error, were used to identify genes whose expression was most associated with phenotype of interest, and an FDR-adjusted p-value of <0.05 was set as a threshold for statistical significance.
Transcription rates were estimated based on the number of nascent unspliced transcripts using EISA [7]. For each gene, we used HT-Seq to compute the average number of reads mapping to the gene's introns (all exonic regions are excluded). Then, this number of intronic reads was divided by the total length of introns to yield a mean intronic coverage that was used as a proxy of the transcription rate.
Functional annotation was performed using Gene Set Enrichment Analysis (GSEA) [31] with all genes preranked by FC as determined by Limma test, or by delta exon/delta intron ratio as determined by EISA, or by Pearson correlation index between delta exon/delta intron ratio and lincRNAs expression. Gene sets were considered significantly enriched if p < 0.05 and FDR < 0.25.
The Pearson correlation was used to determine the relationship between exons/introns and gene expression. All statistical analysis was done with R (version 4.0.3) scripts. The significance of gene set overlap was determined by hypergeometric test.
MicroRNA binding prediction was performed searching in TarBase or LncBase and then intersection of miRNAs were computed by Venn diagram. MiRNA responsive elements (MRE) were calculated by the algorithm RNA22.
Prediction and quantification of circRNAs was carried out by the means of CIRIquant [13] with default parameters. In brief, CIRIquant uses HISAT2 [32] to align the RNA-seq reads to the reference genome and CIRI2 to identify putative circRNAs in the form of BSJ, which are then filtered to reduce the number of false-positives. Since the normalization of circRNA expression values is necessary for the differential expression analysis, TMM normalization factors were extracted from gene expression levels to remove the systematic technical effect of library size. Gene count matrix for the normalization was obtained using the script prepDE.py, from stringTie [33], on the aligned reads. Finally, the voomWithQualityWeights [34] function, from the limma package in Bioconductor [28], is applied to identify the statistical significance of circRNA expression change.

Reverse Transcription of Total RNA to cDNA
Total RNA (500 ng) was processed for each sample by the SuperScript III First-Strand Synthesis SuperMix for qRT-PCR (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. cDNA was stored at −20 • C.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/ncrna7020026/s1, Figure S1: Bimiralisib reduces transcription of genes encoding for proteasome and ribosome components, Figure S2: lncTNK2-2 induces stabilization of p53 by sequestering miR21-3p. Table S1: Details of GSEA of transcriptional changes in U2932 and TMD8, exposed to DMSO or to bimiralisib. Table S2: GSEA of post-transcriptional changes in U2932 and TMD8, exposed for 4 h to DMSO or to bimiralisib. Table S3: Pearson correlation indexes referred to each differentially expressed lincRNAs and differentially stabilized transcripts upon 4 h of bimiralisib exposure in TMD8 and U2932. Table S4: Details of genesets enriched after GSEA of post-transcriptionally modified transcripts ranked by decreasing Pearson correlation index between ∆exon/∆intron ratio and lncTNK2-2:1 and GMDS-AS1 expression Table S5: List of miR-21-3p targets according to LncBase and TarBase.
Author Contributions: N.M. designed and performed experiments, interpreted data, and co-wrote the manuscript; L.C. performed data mining, interpreted data, and co-wrote the manuscript; L.P. performed data mining; C.T. performed experiments; A.R. performed RNA sequencing; N.C. and V.C. designed bimiralisib, R.G. provided advice, F.B. interpreted data and co-wrote the manuscript; S.N. designed the study, performed data mining, interpreted data, and wrote the manuscript. All authors have read and agreed to the published version of the manuscript.