Identification of the whole genome of alternative splicing and RNA-binding proteins involved in nintedanib-induced apoptosis in gastric cancer cells

Cell culture

The human gastric cancer cell lines MKN-28 (Shanghai Fuheng Biotechnology Co., Ltd.) and AGS (Procell Life Science & Technology Co., Ltd) were maintained in 1640 medium supplemented with 10% fetal bovine serum and 1% penicillin-streptomycin. The cells were cultured at 37 °C in a 5% CO₂ environment.

Cell proliferation assay

Cell proliferation assay was performed using the CCK-8 kit (APExBIO, Houston, TX, USA). According to the manufacturer’s instructions, cells were seeded at 6000 cells per well in a 96-well plate. After incubating at 37 °C for 24 h, different concentrations of nintedanib (0, 0.01, 0.1, 1, 10, and 100 μM) were added to gastric cancer cells. Then, 10 μL of the CCK-8 reagent was added to each well at 24, 48, and 72 h after treatment, followed by incubation at 37 °C for 2 h. The optical density (OD) value was measured at 450 nm wavelength using a microplate reader.

EdU assay

The EdU cell proliferation assay was performed using the EdU kit (MeilunBio; EdU-555, China). MKN-28 and AGS cells were seeded at a density of 5 ×10⁴ cells per well in a 24-well plate. After treating with different concentrations of nintedanib (0, 1, 5, 10 μM) for 48 h, the cells were incubated with 50 μM EdU for 2 h and then fixed with 4% paraformaldehyde for 30 min. After permeabilization with 0.3% Triton X-100 for 15 min and washing with PBS, the cells were incubated with 100 μL of Click reaction mixture in the dark for 30 min. Finally, the cells were incubated for 10 min with Hoechst 33342 and observed under a fluorescence microscope to measure the fluorescence intensity.

Flow cytometry

AGS human gastric cancer cell lines were cultured in RPMI-1640 medium (with 10% FBS) at 37 °C and 5% CO₂ until they reached 50–60% confluence. Based on the experimental design, the cells were divided into a control group and treatment groups (nintedanib: 1 µM, 5 µM, 10 µM) and treated for 48 h. The Annexin V-FITC/PI dual staining kit was employed for analysis. A total of 1 × 10⁶ cells were resuspended in 300 µL of PBS, and 5 µL of Annexin V-FITC and 5 µL of PI were added, followed by a 30-minute incubation in the dark. Subsequently, 200 µL of Binding Buffer was added to mix, and flow cytometry analysis was performed immediately. Using the BD FACSCanto II flow cytometer, 30,000 events were recorded, and data were analysed to generate scatter plots distinguishing live cells, early apoptotic cells, late apoptotic cells, and necrotic cells. Results were expressed as percentages, and statistical methods were used to validate the significance of the findings.

RNA-seq

AGS cells were selected for subsequent sequencing analysis. The minimum sample number required for effective data analysis was determined using RNASeqPower (https://rodrigo-arcoverde.shinyapps.io/rnaseq_power_calc/). Estimating that the average depth of coverage from RNA sequencing would be 20 reads per gene, the biological coefficient of variation as 0.2 and alpha was set to 0.05. We determined that each group requires a minimum of 3 replicate samples to detect differential expression, with an effect size of 2 and power of 0.8. Log-phase AGS cells were treated with nintedanib for 48 h, and cellular RNA was extracted. RNA sequencing was performed with the assistance of Genewiz (Azenta Life Sciences, Suzhou, China). A total of 1 µg of total RNA was used for library preparation. Poly(A) mRNA was isolated using Oligo (dT) beads. The mRNA fragments were then subjected to fragmentation under divalent cation and high-temperature conditions. The first-strand cDNA and second-strand cDNA were synthesized using random primers. The purified double-stranded cDNA was then treated to repair the ends, followed by the addition of a dA tail in a single reaction and T-A ligation to add adapters at both ends. Size selection of the adapter-ligated DNA was performed using DNA purification beads. Each sample was PCR amplified using P5 and P7 primers, and the PCR products were validated. Subsequently, libraries with different indices were multiplexed and loaded onto the Illumina HiSeq/Illumina Novaseq/MGI2000 instrument for sequencing, using a 2x150 paired-end (PE) configuration as per the manufacturer’s instructions. The raw data were quality controlled using Cutadapt (V1.9.1) software. Image analysis was performed using HiSeq Control Software version 1.8.4. The DESeq2 Bioconductor package was used for differential expression analysis.

Alternative splicing analysis

To define and quantify alternative splicing (AS) and regulated alternative splicing (RAS) events between the nintedanib-treated group and the control group based on the RNA-seq data, we utilized the SUVA software (Cheng et al., 2021). Statistical analysis was performed on the read counts of all samples, and all detected alternative splicing events (ASEs) were included in the differential analysis (splicing ratio difference is more than 0.15 and P value isless than 0.05). The proportion of events was obtained as the pSAR (percentage of events) using SUVA.

Differentially expressed RBP analysis and RBP–RAS network construction

We screened for all differentially expressed genes between the nintedanib-treated group and the control group (FDR ≤ 0.01, FC ≥ 1.5 or ≤ 1/1.5) and identified an overlap with the existing 2,141 RBP genes. Co-expression analysis was performed for the Differentially expressed RBPs (DE-RBPs) and RAS events (pSAR ≥ 50%) related to cell proliferation and angiogenesis (|Pearson’s correlation| ≥ 0.99, p-value ≤ 0.01), and a regulatory network of RBP–RAS was constructed.

Functional enrichment analysis

To distinguish the functional categories of differentially expressed genes (DEGs), the KEGG Orthology-Based Annotation System (KOBAS) 2.0 (Xie et al., 2011) server was utilized to identify GO terms and KEGG pathways. The enrichment of each term was determined using the hypergeometric test and Benjamini–Hochberg FDR controlling procedure.

Other statistical analyses

Principal component analysis (PCA) analysis was performed by R package factoextra (https://cloud.r-project.org/package=factoextra) to show the clustering of samples with the first two components. The quantified experimental data of CCK-8 and EdU were presented as mean ± standard deviation (SD), and statistical analysis was performed using GraphPad Prism 9.0 software (GraphPad, La Jolla, CA, USA). Student’s t-test was used for comparison between two groups, and one-way ANOVA was used for comparison among multiple groups. A p-value <0.05 was considered statistically significant.

Nintedanib inhibits the proliferation of gastric cancer cells

To determine whether nintedanib has an inhibitory effect on the growth of gastric cancer cells, AGS and MKN-28 gastric cancer cells were treated with different concentrations of nintedanib for 48 h, and cell viability was assessed using the CCK-8 assay. The results demonstrated that nintedanib inhibited the viability of both gastric cancer cell lines in a dose-dependent manner (Figs. 1A–1B). The results revealed that after 48 h of treatment with nintedanib, the IC50 value of AGS cells was 6.51 ± 0.49 μM, and the IC50 value of MKN-28 cells was 8.07 ± 0.13 μM. To visualize the effects of nintedanib on gastric cancer cell proliferation, we examined the changes in EdU-positive cells after treatment with different concentrations of nintedanib (1, 5, 10 μM) for 24 h. The EdU labeling showed a dose-dependent decrease in the percentage of EdU-positive cells compared to the control group (Figs. 1C–1F). Based on the observed strong inhibitory effect of nintedanib on the proliferation of gastric cancer cells, we further investigated its impact on apoptosis in these cells. The results indicated that after 48 h of nintedanib treatment, gastric cancer cells exhibited concentration-dependent apoptosis (Figs. 1G–1H). Taken together, these data indicate that nintedanib significantly inhibits the proliferation of gastric cancer cells in vitro and induces apoptosis.

Identification of highly conserved regulated AS events responded to nintedanib therapy in gastric cancer cells

Log-phase AGS gastric cancer cells were incubated with 5 μM nintedanib for 48 h, and an untreated control group was also set up. RNA was extracted and subjected to RNA-Seq analysis. SUVA (Cheng et al., 2021) software was used to analyze the RNA-Seq data from the nintedanib-treated and control groups. We identified a significant number of splicing events in both datasets and determined the differential events among them (Fig. S1A). SUVA has reclassified the types of alternative splicing events (Fig. 2A). The number of five types of AS events identified by SUVA in RNA-Seq data and the number of RAS events analyzed for differences were counted. The main types of RAS events that occurred were alt5p and alt3p events (Figs. 2B and S1B). After mapping the differential splicing events identified by SUVA to classic splicing events, we found that A5SS, cassette exon, ES, and A3SS were the most frequently occurring alternative splicing events (Figs. 2C and S1C). Furthermore, approximately 80% of the AS events among these RAS events were complex splicing events (Figs. 2D and S1D), indicating the complexity of splicing events during nintedanib action in gastric cancer cells.

Although there were relatively many splicing events, not all of them can function at the transcriptional level. The splicing rate pSAR represents the proportion of ASEs in the transcriptome of a gene relative to the entire transcriptome. We selected dominant splicing events (pSAR >= 50%) from the RAS events (915 RAS events in total) for further analysis (Fig. 2E). Using their splicing ratio values, we performed principal component analysis (PCA), which showed that the nintedanib treatment and control groups were separated along the first principal component (Fig. 2F). Heatmap visualization was used to show the selected RAS events with pSAR >= 50% (Fig. 2G). GO and KEGG pathway analyses were performed on the genes containing differentially spliced events between the two groups, and the enriched entries were visualized (Figs. 2H, S1E). Finally, we found that the DEGs were mainly related to the biological pathway of cell apoptosis. This aligns with our findings regarding the apoptosis induced by nintedanib, as presented in Fig. 1.

Alternative splicing genes involved in the apoptosis pathway

In the previous section of our results, we focused on the differential ASEs that were closely related to cell apoptosis between the two groups. Further studies were conducted on the RASEs in the cell apoptosis pathway, and it was found that there were seven variable splicing genes shared by the experimental and control groups. Among them, MAP3K10, BBC3, BAG6, NR4A1, and IFI27 were significantly upregulated after nintedanib treatment. The heatmap showed the ratio values of the relevant RASGs in the pathway (Fig. 3A). The splicing event of the gastric cancer drug resistance-related gene IFI27 significantly increased the selection proportion of the 3′end exon in the nintedanib treatment group (Fig. 3B). The selection proportion of the 5′end exon of the cell apoptosis-related gene BBC3 was significantly higher in the nintedanib treatment group than in the control group (Fig. 3C). Figures S2A–S2B show the differential expression levels and read distribution of the genes MAP3K10 and BAG6 in the two groups.

Construction of co-expression network between RNA binding proteins and RAS of genes involved in apoptotic process pathway

RBPs are typically identified as upstream regulatory factors in AS events during tumor development (Hentze et al., 2018). Differences in RBPs between the nintedanib-treated group and the control group may have led to the differential AS events currently being presented. We used DESeq2 to analyze all differentially expressed genes (FDR ≤ 0.01, FC ≥ 1.5 or ≤ 1/1.5) between the nintedanib-treated group and the control group (Fig. 4A). When we overlapped all differentially expressed genes with the currently known 2141 RBPs, we found that 134 differentially expressed RBPs (DERBPs) were involved in the process of nintedanib drug treatment (Fig. 4B). The heatmap showed that most of the differentially expressed RBPs were downregulated in the nintedanib group (Fig. 4C). GO functional enrichment analysis of all upregulated DERBPs between the two groups showed that these RBPs were significantly enriched in biological pathways related to type I interferon signaling pathway, innate immune response, and immune system process (Fig. S3A). Downregulated DERBPs were significantly enriched in biological pathways related to rRNA processing, ribosome biogenesis, tRNA processing, and tRNA modification (Fig. S3B). KEGG functional enrichment analysis showed that upregulated DERBPs were mainly enriched in biological pathways such as cell adhesion molecules (CAMs) and Inflammatory mediator regulation of TRP channels (Fig. S3C). The significantly enriched items in the KEGG analysis of downregulated DERBPs were spliceosome, ribosome biogenesis in eukaryotes, and the hippo signaling pathway, among others (Fig. S3D).

To further explore the interaction relationship between DERBPs and differential AS events, we conducted co-expression analysis using the pSAR values of deRBP genes and relevant differential AS events. We identified typical DERBP genes that were co-expressed with differential AS events, extracted the genes where the co-expressed differential splicing events occurred, and conducted GO functional analysis. It was found that co-expressed genes were mainly enriched in pathways related to positive regulation of transcription by RNA (Fig. 4D). Co-expression analysis was also performed on differential RBPs and differential events of cell proliferation and angiogenesis-related genes that we focus on, and a regulatory network of RBP-RAS was constructed (Fig. 4E). It was found that the relevant genes included NOC2L, TAGLN, DHX9, and other abnormally expressed DERBPs, which could regulate the activity of downstream cell proliferation and angiogenesis-related genes such as NR4A1, BBC3, and IFI27. Box plots showed the expression levels of these DERBPs (Fig. S3E). We focused on the DHX9-NR4A1 relationship and presented the distribution plot and ratio value of DHX9 expression and clualt5p12044:NR4A1 reads (Fig. 4F).