We downloaded and obtained the gene expression (RNA-seq) and survival information datasets of 24 cancer types from TCGA project, including tumor and normal samples, from the UCSC XENA database (https://xenabrowser.net/datapages/). Then we analyzed the expression levels of RBM6. The survival package in R (version 4.2.3) was used to perform survival modeling and Kaplan-Meier (KM) analyses.

The short hairpin RNA (shRNA) duplexes of RBM6 were purchased from GenePharma (Suzhou, China). The negative control shRNA (shNC) sequence was: 5^′-TTCTCCGAACGTGTCACGT-3^′ (sense), whereas the shRNA targeting RBM6 (shRBM6) sequence was: 5^′-TTCAGGTTCCTTATCAGCTGG-3^′ (sense).

The HeLa cell line (CL-0101; Procell, Wuhan, China) was validated by STR profiling and tested negative for mycoplasma, and then cultured at 37 °C with 5% CO₂ in RPMI-1640 medium (CM-0063; Procell) supplemented with 10% fetal bovine serum (10091148; Gibco, Beijing, China), 100 µg/mL streptomycin, and 100 U/mL penicillin (SV30010; Hyclone, Logan, UT, USA). Then shRBM6 and shNC were transfected into HeLa cells using Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocol. After 48 h, the HeLa cells were collected for subsequent experiments and analysis.

The cell proliferation assay was performed as previously described [25] using the Cell Counting Kit-8 assay (40203ES76; Yeasen, Shanghai, China). Briefly, shRBM6- and shNC-transfected HeLa cells were seeded in culture plates at a density of 10,000 cells/well. After 24, 48, and 72 h of incubation, 10 µL CCK-8 solution was added to the cells followed by incubation for an additional 3 h at 37 °C. Then the optical density (OD) was measured at an absorbance of 490 nm using a microplate reader (ELX800; Biotek, Winooski, VT, USA). The OD values were used to calculate the proliferation rate.

Apoptotic cells were stained using the Annexin V-PE/7-AAD Apoptosis Detection Kit (KGA-1017; KeyGEN, Nangjing, China) after 72 h of transfection. Then the shRBM6- and shNC-transfected cells (2 $\times$ 10⁵) were incubated with 5 µL 7-AAD suspended in 50 µL binding buffer per sample for 15 min, and then incubated with 1 µL Annexin V-PE (BD Pharmingen, San Diego, CA) and 450 µL binding buffer per sample for 15 min. Flow cytometry was used to detect cell apoptosis levels for each sample.

DNA was removed by RQ1 DNase (Promega Corp., Madison, WI, USA) to extract total RNA. After quality and quantity assessment of RNA by the SmartSpec Plus spectrophotometer (BioRad Laboratories, Inc., Hercules, CA, USA), we used 1 µg integral RNA for RNA-seq library construction with the NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (#E7420; New England Biolabs, Ipswich, MA, USA) by purifying polyadenylated RNAs. Then the RNAs were fragmented and converted into cDNAs. After end repair and A tailing, the DNAs were ligated with the Diluted NEBNext® Adaptor (New England Biolabs, Ipswich, MA, USA). The ligation products corresponding to 300–500 base pairs and digested with heat-labile uracil-DNA glycosylase, were processed and stored at –80 °C for sequencing. For next-generation sequencing, the libraries were applied to the Illumina NovaSeq 6000 sequencing system (Novogene, Beijing, China) for 150 nucleotide (nt) paired-end sequencing.

We initially discarded raw reads with N bases, removed adaptors and low-quality bases using FASTX-Toolkit (version 0.0.13, https://github.com/Debian/fastx-toolkit), and finally kept reads longer than 16 nt. The quality filtered reads were aligned to the GRCh38 genome by HISAT2 [26]. Mapped reads with only one genomic location were used to calculate the read count of genes and the normalized fragments per kilobase of transcript per million fragments mapped (FPKM).

We used the splicing site usage variation analysis (SUVA, version 1.0.0, https://github.com/ablifedev/SUVA) pipeline [27] to calculate and quantify the AS events (ASEs) and RASEs (p-value 0.05) between shRBM6 and shNC samples, including the alternative 5^′ splice site (alt5p), alternative 3^′ splice site (alt3p), overlapped (olp), intron retention (ir), and contained splice junction (contain), which are fully explained in the SUVA paper. The frequency and proportion of each SUVA ASE read (pSAR) were used to analyze statistical significance. Meanwhile, we also calculated the canonical AS types, including intron retention (IntronR), alternative 5^′ splice site (A5SS), alternative 3^′ splice site (A3SS), exon skipping (ES), mutual exclusive exon (MXE), cassette exon, A5SS&ES, A3SS&ES, 5^′ MXE (5pMXE), and 3^′ MXE (3pMXE).

The differentially expressed genes (DEGs) were predicted using the DEseq2 software package (version 1.42.0, https://github.com/thelovelab/DESeq2) [28], which analyzes the differential expression of genes. A two-fold change (FC) and false discovery rate (FDR) 0.05 were selected as the cutoffs.

The iRIP-seq technique was conducted as previously described [29]. Briefly, after irradiation at 400 mJ/cm², the HeLa cells were lysed in cold wash buffer supplemented with RNase inhibitor (2313U; Takara Bio Inc., Shiga, Japan) and protease inhibitor cocktail (B14001; Bimake, Shanghai, China) on ice for 30 min. Then the RNA/protein mixture was vibrated vigorously and centrifuged to remove cell debris. Then RNAs were digested by MNase and RNaseT1 (Thermo Fisher Scientific, Waltham, MA, USA). Ethylenediamine tetra-acetic acid (EDTA, Sangon Biotech., Shanghai, China) was added to stop the digestion. The supernatant was incubated with 13 µg RBM6 antibody (14360-1-AP; Proteintech, Hongkong, China) and control IgG-antibody (AC005; ABclonal, Wuhan, China) overnight at 4 °C. After applying the mixture to the magnet and removing the supernatants, the beads were sequentially washed twice with lysis buffer. Next, the suspension was incubated in a heat block to release the immunoprecipitated RBPs with crosslinked RNAs and vortexed for 30 min. Proteinase K (Sangon Biotech., Shanghai, China) was added to the 10% input (without being immunoprecipitated) and immunoprecipitated RBM6 with crosslinked RNAs. The RNAs were purified with phenol: chloroform: isopentyl alcohol (25:24:1 pH 5; Solarbio, Beijing, China).

The cDNA libraries were prepared with the KAPA RNA Hyper Prep Kit (KK8541, Roche, Basel, Switzerland) according to the manufacturer’s instructions. Then the libraries were prepared following the manufacturer’s instructions and applied to the Illumina NovaSeq sequencing system for 150 nt paired-end sequencing.

We aligned the quality filtered reads onto the genome by HISAT2 (version 2.2.1, https://github.com/DaehwanKimLab/hisat2) [26], kept the uniquely aligned reads, and removed the PCR duplicates with the same aligning position. Then, Piranha was used to perform peak calling between RBM6-IP and input samples [30]. The target genes of RBM6 were finally determined by the peaks. The binding motifs of RBM6-bound peaks were called by HOMER software (version 4.10, https://github.com/javrodriguez/HOMER) [31].

Quantitative PCR (qPCR) was conducted as previously described [32]. Briefly, cDNAs were synthesized on a thermocycler (T100; Bio-Rad, Hercules, CA, USA) using a Reverse Transcription kit (R323-01; Vazyme Biotech, Nanjing, China), and qPCR was performed on the ABI QuantStudio 5 Real-Time PCR System (Thermo Fisher Scientific) following the standard operational process. The concentration of each transcript or splicing event was normalized to glyceraldehyde-3-phosphate dehydrogenase (GAPDH), and their levels were calculated using the 2^{-Δ⁢Δ⁢CT} method [33]. The unpaired and two-tailed Student’s t-test were used to compare the significant differences of each transcript or splicing event, the primer sequences are presented in Supplementary Table 1. Each sample had three technical replicates.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genome (KEGG) databases were used to calculate functions and enriched score and significance for each functional pathway using the KOBAS 2.0 server [34]. The statistical significance of each pathway was defined using hypergeometric tests and the Benjamini-Hochberg FDR controlling procedure (corrected p 0.05).

The two-tailed unpaired Student’s t-test was used for other comparisons that were not defined above between the shRBM6 and control groups. The threshold for a statistically significant difference was set as p 0.01 or FDR 0.05.

Previous studies have demonstrated that RBM6 can serve as a tumor suppressor and repress tumor growth and progression [22, 23, 35]. We compared RBM6 expression in tumor vs. adjacent normal samples in 24 cancer types from TCGA. Among the 24 cancer types, RBM6 expression was upregulated (p 0.05) in 7 cancer types and downregulated in 5 cancer types (Fig. 1A). Prognosis analysis revealed that higher RBM6 expression was significantly and positively associated with longer overall survival in 12 of the 18 cancer types (p 0.05), including head and neck squamous cell carcinoma (HNSCC; Fig. 1B), lung adenocarcinoma (LUAD; Fig. 1C), pancreatic adenocarcinoma (PAAD; Fig. 1D), and other cancer types (Supplementary Fig. 1). Meanwhile, higher expression of RBM6 was associated with a worse prognosis in the other six cancer types (Supplementary Fig. 1). These results demonstrate that RBM6 expression is related to the prognosis of multiple tumors and may play an important regulatory role in the pathogenesis or development of many tumors.

Fig. 1.

Higher RBM6 expression was associated with better prognosis in several tumor types by pan-cancer analysis. (A) The boxplot shows RBM6 expression pattern in tumor vs. normal tissues in the 24 cancer types from TCGA. RBM6 expression was up-regulated in seven cancer types, and down-regulated in five cancer types. * p 0.05; ** p 0.01; *** p 0.001; Student’s t-test. (B) The survival curve displays the prognosis of HNSCC from TCGA. (C) The same as (B) but for LUAD. (D) The same as (B) but for pancreatic adenocarcinoma. RBM6, RNA-binding motif protein 6; TCGA, The Cancer Genome Atlas; HNSCC, head and neck squamous cell carcinoma; LUAD, lung adenocarcinoma.

To further investigate the functions of RBM6 in tumor cells, we constructed an RBM6 knockdown (shRBM6) cell model by transfecting an RBM6-shRNA sequence into HeLa cells; HeLa cells transfected with an empty vector were served as the NC. We assessed the knockdown efficiency of shRBM6 by the real-time polymerase chain reaction (RT-PCR) and found that RBM6 was significantly down regulated in the shRBM6 samples (p 0.0001; Fig. 2A). We performed flow cytometry analysis and the CCK8 assay to determine whether RBM6 could modulate cell apoptosis (Fig. 2B) and proliferation of HeLa cells, respectively. The results demonstrated that the apoptotic cell percentage was significantly decreased in shRBM6 samples (p 0.01; Fig. 2B and Supplementary Fig. 2), whereas the cell proliferation rate reflected by the optical density was significantly increased in the shRBM6 samples (p 0.01; Fig. 2C). These results indicate that RBM6 has a suppressive role in the tumor phenotype or function of HeLa cells.

Fig. 2.

RBM6 knockdown promoted proliferation and inhibited apoptosis. (A) The bar plot shows reverse transcription and quantitative polymerase chain reaction (RT-qPCR) validation of RBM6 knockdown efficiency. (B) The bar plot shows the decreased cellular apoptosis rate in the shRBM6 group. (C) The bar plot indicates the increased cellular proliferation level in the shRBM6 group. The values are presented as the mean $\pm$ standard deviation. ** p 0.01; **** p 0.0001; Student’s t-test.

To further elucidate how RBM6 modulates the synthesis and processing of RNA in HeLa cells, we performed unbiased RNA-seq analysis for shRBM6 and normal controls (NC) HeLa cells by capturing the total polyadenylated RNAs, with three biological replicates for each group. After quality filter and genome alignment of the raw RNA-seq data, more than 95% of the aligned reads had a uniquely aligned location in the genome and were used for further analyses (Supplementary Table 2). By calculating the read count and FPKM value for each gene, we detected a total of 28,218 expressed genes from the six RNA-seq samples (Supplementary Table 3). We presented the FPKM values of RBM6 and confirmed its successful knockdown from the RNA-seq analysis results (Fig. 3A).

Fig. 3.

RBM6 modulated the AS pattern of genes enriched in DNA damage response- and cell survival-related pathways. (A) Boxplot shows the Fragments per kilobase per million mapped (FPKM) of RBM6. (B) Barplot shows the number of RASEs detected by SUVA between shRBM6 and NC. (C) Barplot presents the number of RASEs by dividing them into complex and simple ASEs using SUVA. (D) Barplot shows the corresponding classical ASE types for the RASEs detected by SUVA. (E) Barplot shows the number of RASEs with increased pSAR values. RASEs with pSAR $\ge$50% are labeled as brown. (F) PCA clustering of RASEs with pSAR $\ge$50%. (G) The bubble plot displays the top 10 GO pathways for RASGs (pSAR $\ge$50%) in the shRBM6 vs. NC groups. (H) The table shows proliferation and apoptosis pathways in the shRBM6 group. (I) Read distribution diagram shows clualt3p78233 MARCHF7. The boxplot in the right panle shows the splicing ratio of clualt3p78233 MARCHF7. AS, alternative splicing; RAS, regulated AS; RASEs, RBM6-regulated AS events; SUVA, site usage variation analysis; NC, negative control; ASEs, AS events; pSAR, proportion of each SUVA ASE read; PCA, principal component analysis; GO, Gene Ontology; RASGs, The genes involve RASE.

As an AS factor, RBM6 is supposed to modulate the tumorigenesis or development probably by playing an essential role in regulating the AS profile. Since cancers undergo complex splicing events [36], we detected thousands of ASEs in RNA-seq using SUVA [27] (Supplementary Fig. 3A). Among these ASEs, 1897 high-confidence RASEs were identified by applying stringent cutoffs of p $\le$ 0.05 and splicing ratio difference $\ge$0.15 (Supplementary Table 4, Fig. 3B). There were 5 RASE types: 852 alt3p, 818 alt5p, 95 contain, 19 IR, and 113 olp. Strikingly, 11.4% of SUVA RASEs were simple and 88.6% were complex. Complex events contained splicing junctions that were involved in two or more splicing events (Fig. 3C), indicating that RBM6 prefers to regulate the complex ASEs in HeLa cells. We then corresponded the SUVA RASEs to classical splicing events, with A5SS, A3SS, ES, and cassette exon events accounting for a large proportion (Fig. 3D and Supplementary Fig. 3B). One splicing event was involved in two or more transcripts, whereas only one transcript was probably highly expressed and dominantly spliced for this gene. Thus, we decided to predict the dominant transcript undergoing AS that could be quantified by the “pSAR” value using SUVA. We calculated the number of RASEs with different pSAR values, and found that 598 dominant splicing events had pSAR $\ge$50% (Fig. 3E). These 598 RASEs were selected for further analyses. Meanwhile, we performed principal component analysis (PCA) for RASEs and found that the shRBM6 group and the control group were clearly separated by the splicing ratio of these 598 RASEs (Fig. 3F and Supplementary Fig. 3C).

These 598 RBM6-RASEs were distributed in 463 genes, which were enriched in GTPase and transcription pathways (Fig. 3G). We found that RBM6 regulates the AS of many genes involved in proliferative and apoptotic processes (Fig. 3H). The reads distribution diagram showed an RASE in the membrane associated ring-CH-type finger 7 (MARCHF7) gene, and the splicing model showed the significantly dysregulated splicing pattern of this event (Fig. 3I). In such a complex AS region, SUVA identified the RASE that was significantly altered between the shRBM6 and control samples. As shown in Fig. 3I, the exons shown as black boxes were skipped after knocking down RBM6. Another RASE, located in PRKCD, showed a tendency for intermediate exons (Supplementary Fig. 3D). These results indicate that RBM6 may affect the survival and development of cancer cells by regulating variable splicing of genes involved in cell proliferation or apoptosis.

To investigate RBM6-mediated transcriptional regulation, the DESeq2 package was used to identify the DEGs between shRBM6 and NC samples using FC $\ge$2 and FDR 0.05 as criteria [28] (Supplementary Fig. 4A). We finally obtained 202 up-regulated and 1348 down-regulated DEGs (Fig. 4A, Supplementary Table 5), indicating that RBM6 knockdown tended to result in repression of gene expression as there were more downregulated than upregulated DEGs.

Fig. 4.

RBM6 significantly altered the expression of genes associated with tumor development, which may be mediated by TFs from RBM6-RASGs. (A) Volcano plot presents the upregulated- and downregulated-DEGs in the shRBM6 vs. NC groups using DESeq2. (B) The top 10 enriched GO pathways of up-regulated DEGs in the shRBM6 group. (C) The top 10 enriched GO pathways of down-regulated DEGs in the shRBM6 group. (D) Venn diagram shows the upregulated- and downregulated-genes and the known oncogenes or tumor suppressor genes. (E) Venn diagram shows the overlap of TFs and RASGs. Venn diagram showing the overlap of TF-targets and DEGs (above). Venn diagram shows TF-targets and TF-related DEGs overlay (below). (F) Network showing TF-target pairs were screen out and the co-disturbed. (G) Bar plot showing the RT-qPCR results. TF, ranscription factors; DEGs, differentially expressed genes (** p 0.01; *** p 0.001; **** p 0.0001.)

To further identify the potential functions of these DEGs, we performed GO (Supplementary Table 6) and KEGG (Supplementary Table 7) enrichment analyses. The top 10 upregulated- and downregulated GO terms and the top 10 upregulated- and downregulated KEGG pathways were shown in Fig. 4B,C and Supplementary Fig. 4B,C, respectively. The up-regulated genes in the shRBM6 cells are primarily enriched in type I interferon (IFN) signaling pathways, defense response to virus, and the innate immune response (Fig. 4B). Type I IFN is tightly associated with viral infection and has double-edged effects on tumor development and resistance to treatment [37, 38]. Interestingly, the down-regulated DEGs were mainly associated with cell cycle and DNA damage response (DDR) pathways, including cellular responses to DNA damage stimulus, cell cycle, and DNA repair. Positive regulation of GTPase activity, RNA splicing, and cilium assembly were also enriched in the down-regulated DEGs (Fig. 4C). Cellular responses to DNA damage stimulus are closely related to apoptosis in multiple cancers [39, 40]. These results suggest that shRBM6 can activate the type I IFN signaling pathway and suppress the DDR, ultimately inhibiting the apoptosis of cancer cells. Meanwhile, we also found that the expression of many oncogenes or tumor suppressor genes was detected in the DEGs caused by shRBM6 (Fig. 4D and Supplementary Fig. 4D).

RBM6 is an important AS regulator. The genes involve RBM6-RASE (RASGs) were enriched in transcriptional regulatory pathways (Fig. 3G). We speculate that the AS of many transcription factors (TFs) is regulated by RBM6, which in turn has an impact on the transcriptome. In total, 63 of 464 RASGs were identified as TFs (Fig. 4E). We conducted a correlation analysis of the expressions of DEGs and the splicing ratios of RASEs from these 63 TFs (Pearson’s correlation coefficient $\ge$0.8, p $\le$ 0.01). Meanwhile, TF-target pairs of these 63 TFs with RASEs were selected from three databases (TRRUST, Encode, and TFBSDB). Finally, 202 TF-target pairs were identified and the co-disturbed network was constructed (Fig. 4F). To confirm the dysregulation by shRBM6, we performed RT-qPCR and validated the decreased expression of interleukin 24 and ETS1, and increased expression of interferon regulatory factor 7 (IRF7) (Fig. 4G). Based on these results, we propose that RBM6 may indirectly play a regulatory role at the transcriptional level through regulating the splicing of these TFs (Supplementary Fig. 4E).

To study the RBM6-RNA interactions and the mechanism underlying the function of RBM6, we obtained the global RBM6-RNA interaction profile in HeLa cells by iRIP-seq [29, 41]. We used ultraviolet cross-linking and RNA digestion to capture the exact RNA-protein interaction locations. The iRIP cDNA libraries were sequenced on the Illumina platform, generating abundant iRIP-seq reads for RBM6 IP and input control samples with two biological replicates for each (Supplementary Table 8). We processed and aligned the cDNA sequence reads to the human GRCh38 genome using HISAT2. The fractions of the aligned reads were 78.14–86.99%, from which we selected the uniquely aligned reads for following analyses (Supplementary Table 8). By classifying the uniquely mapped reads to the genomic location, we found that the RBM6-associated reads showed a higher percentage in the 5^′UTR and intergenic regions but lower percentage in the intron and coding sequence (CDS) regions (Fig. 5A). The normalized reads density along the uniformed gene region showed that RBM6 had a stronger tendency for the 5^′UTR regions compared with the CDS and 3^′UTR regions (Fig. 5B).

Fig. 5.

RBM6-mRNA interaction map in HeLa cells. (A) Barplot shows the distribution of sequencing reads on different genome regions. (B) Coverage intensity analysis of iRIP-seq reads with the length of the transcription unit was carried out. The 5^′UTR, CDS, and 3^′UTR of the gene were divided into 100 evenly, each of which was called a bin. The sum of the average reads in each bin was calculated to obtain the overall reads coverage. (C) Venn diagram shows the common and unique peaks in two IP replicate treatment groups by using the Piranha method. The peaks with overlapping positions are clustered in the two replicates, and the overlapping peaks are classified under the same cluster. (D) The proportion of peak reads in different regions of the genome. (E) The top five enriched motifs in the overlapped RBM6 peaks. (F) RASGs co-expressing with RBM6 were screened. Venn diagram shows the peaks bound by RBM6-RASGs. (G) Reads plot shows one RASE in the ZNF276 and RBM6 binding profiles. iRIP-seq, RNA immunoprecipitation coupled with sequencing; CDS, coding sequence; IP, immunoprecipitation; ZNF276, zinc finger protein 276.

To further elucidate the characteristics of the RBM6-RNA interaction, we called the RBM6-binding peaks from iRIP-seq reads using Piranha. In total, 584 shared peaks between the two replicates were identified (Fig. 5C). The genomic locations of these peaks showed that most of these peaks came from the antisense, CDS and intron regions (Fig. 5D). Further analysis showed that RBM6 binding reads were slightly and specifically enriched in the transcription start site, whereas peak reads were significantly enriched in the intron region (Supplementary Fig. 5A,B). Next, we applied the HOMER algorithm to identify the enriched motif sequences over-represented in the RBM6 peaks. Consistent with a previous report [21], the CUCU motif was enriched. However, the top motif was GGCGAUG, indicating that the motif selection could be cellular heterogeneity (Fig. 5E).

Functional clustering by GO was performed to annotate the functions of the RBM6-bound genes containing at least one RBM6 peak; the enrichment fell into the cell cycle, cell differentiation, and cell migration (Supplementary Fig. 5C). As shown in Fig. 5F, 39 RASGs overlapped with RBM6-bound genes, indicating that RBM6 might directly bind to the RNAs of these targets and regulate their AS pattern. In addition, several genes were found to be involved in proliferation or apoptosis, including cyclin-dependent kinase 13 (CDK13), laminin subunit alpha 5 (LAMA5), obscurin (OBSCN), SMAD Family Member 6 (SMAD6), and transcription factor AP2 (TFAP2A). Four TFs, AT-Hook DNA Binding Motif Containing 1 (AHDC1), SMAD6, TFAP2A, and zinc finger protein 276 (ZNF276), were also directly bound by RBM6. Visualization of the reads profile showed an RASE in ZNF276 and RBM6-binding profiles (Fig. 5G). RBM6 specifically binds to the first 5^′ exon at the proximal end, and this exon is more selected in the control group, whereas the first 5^′ exon at the distal end is more selected after RBM6 knockdown.