SERBP1 affects the apoptotic level by regulating the expression and alternative splicing of cellular and metabolic process genes in HeLa cells

SERBP1 sequences cloning and plasmid construction

The primer pairs for Hot Fusion were designed with CE Design (V1.04). Fragments of the SERBP1-specific sequence and the pIRES-hrGFP-1a vector sequence were included in each primer.

Forward primer: agcccgggcggatccgaattcATGCCTGGGCACTTACAGGA

Reverse primer: gtcatccttgtagtcctcgagAGCCAGAGCTGGGAATGCC

The pIRES-hrGFP-1a vector was then digested for 2–3 h at 37 °C by XhoI (NEB) and EcoRI. Then, the enzyme-digested vector was leaked on 1.0% agarose gel and purified with Qiagen column kit. After isolating RNAs from HeLa and HEK293T cells with TRIzol, the oligo dT primer was utilized to reverse-transcribe purified RNA into cDNA. The digested vector and inserted sequences were prepared for ligation with a ClonExpress® II One Step Cloning Kit (Vazyme). Then, the plasmids were transformed into Escherichia coli strain that was plated onto LB agar plates containing 1 µL/ml ampicillin and incubated overnight at 37 °C. Universal primers (complementary to the backbone vector) were used to screen the colonies, which were finally verified by Sanger sequencing.

Assessment of SERBP1 expression and apoptosis experiment

We used RT-qPCR experiment to assess the expression levels of genes, which was performed on Bio-Rad S1000 (DBI Bioscience, Shanghai, China). Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as internal control to assess the effects of SERBP1 overexpression (SERBP1-OE). Information of primers is presented in Table S1. RNA isolation and cDNA synthesis procedures were performed following the published method (Li et al., 2018). The 2^−ΔΔCT method was used to normalize the mRNA expression levels (Livak & Schmittgen, 2001). Statistical significance was analyzed with paired Student’s t-test. Western blotting (WB) experiment was also performed to assess SERBP1 overexpression according to the published method (Li et al., 2018) with 25 μg protein as input. Anti-FLAG (1:2,000 dilution; polyclonal antibody; cat. no. 2368S; CST), Anti-CASP3 (1:2,000 dilution; polyclonal antibody; cat. no. 9664P; CST), Anti-GAPDH (1:5,000 dilution; polyclonal antibody; cat. no. 60004-1-Ig; Proteintech Group, Rosemont, IL, USA), and anti-ACTIN (1:2,000 dilution; polyclonal antibody; cat. no. AC026; ABClonal, Wuhan, China) were used as antibody. The raw data for RT-qPCR and WB were shown in Table S2 and Figs. S2A, S2B and S2E, respectively. For apoptosis experiment, flow cytometric analysis was performed by staining viable cells with 7-amino actinomycin D and FITC-conjugated Annexin V (4A Biotech Co., Ltd., Beijing, China). We then calculated the percentage of apoptotic cells by summing the right two quadrants.

RNA-seq library construction

Two biological replicates SERBP1-OE and control (Ctrl) cells were prepared for total RNA extraction. RNA was purified twice with phenol-chloroform treatment, and genomic DNA was removed with RQ1 DNase (Promega, Madison, WI, USA). The quality, quantity, and integrity of RNA were verified by absorbance at 260 nm/280 nm (A260/A280) (SmartSpec Plus; Bio-Rad, Hercules, CA, USA) and 1.5% agarose gel electrophoresis, respectively.

RNA-seq libraries were prepared by VAHTS Stranded mRNA-seq Library Prep Kit (Vazyme, Nanjing, China) with one 1 microgram RNA as input, which was conducted in accordance with the published method (Li et al., 2018). We used Illumina HiSeq X Ten system for sequencing with 150 nt paired-end as parameter.

RNA-seq raw data clean and alignment

FASTX Toolkit (version 0.0.13) was used to remove reads with N bases and to trim the adaptors and low-quality bases. Finally, short reads (<16 nt) were removed. Then, high-quality reads were aligned onto the GRCh38 genome by TopHat2 (Kim et al., 2013), allowing up to four mismatches. Fragments per kilobase of transcript per million fragments mapped (FPKM) value was calculated by uniquely mapped reads (Trapnell et al., 2010). The statistical power of this experimental design, calculated in RNASeqPower was 0.78.

Package edge R (Robinson, McCarthy & Smyth, 2010) from R Bioconductor was used to identify differentially expressed genes (DEGs) with false discovery rate <0.05 and fold change >2 or <0.5 as criteria.

Alternative splicing prediction and analysis

We used ABLas pipeline (Xia et al., 2017) to define and quantify alternative splicing events (ASEs) and SERBP1-regulated alternative splicing events (RASEs). Ten types of ASEs were detected by ABLas, including alternative 5′ splice site (A5SS), alternative 3′ splice site (A3SS), cassette exons, exon skipping (ES), intron retention (IR), mutually exclusive exons (MXE), mutually exclusive 5′UTRs (5pMXE), mutually exclusive 3′UTRs (3pMXE), A3SS&ES and A5SS&ES. With a false discovery rate cutoff of 5%, significance of the ratio alteration of AS events was calculated by Student’s t-test to assess SERBP1-regulated ASEs.

Validation of DEGs and AS events

Reverse transcription and quantitative polymerase chain reaction (RT-qPCR) experiment was conducted to confirm the validity of DEGs. Information of primers was presented in Table S1. Detailed procedures of RT-qPCR could be found in the published paper (Li et al., 2018). Moreover, RT-qPCR was also used to validate the RASEs. The primers for detecting ASEs are shown in Table S1. The detailed primer design method could be found in the published paper (Liu et al., 2020).

Other statistical analysis

The KOBAS 2.0 server was used to identify the Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of DEGs and RASGs to determine relevant functional categories (Xie et al., 2011). The hypergeometric test and Benjamini-Hochberg false discovery rate (FDR) controlling procedure were used to define the enrichment of each term. Expression levels of SERBP1 in multiple cancers were analyzed using GEPIA2 online tool (Tang et al., 2019). Kaplan-Meier (K-M) plotter tool (Szasz et al., 2016) was used to analyze the prognosis results of SERBP1 in multiple cancer types.

SERBP1 showed dysregulated expression pattern and affected prognosis in multiple cancer types

We first investigated the expression pattern of SERBP1 in all the cancer types from the cancer genome atlas (TCGA) database using GEPIA2 web server (Tang et al., 2019). The analysis results showed that the expression values of SERBP1 were significantly higher in several cancers, including cervical squamous cell carcinoma (CESC), cholangiocarcinoma (CHOL), colon adenocarcinoma (COAD), lymphoid neoplasm, diffuse Large B-cell lymphoma (DLBC), glioblastoma multiforme (GBM), Brain lower grade glioma (LGG), pancreatic adenocarcinoma (PAAD), rectum adenocarcinoma (READ), stomach adenocarcinoma (STAD), thymoma (THYM), and uterine carcinosarcoma (UCS). In no cancer did SERBP1 show significantly lower value, although it was highly expressed in acute myeloid leukemia (LAML) normal tissues (Fig. 1A). We then analyzed the prognosis effects according to SERBP1 expression levels in cancers. Using K-M plotter tool (Szasz et al., 2016), we found that higher expression levels of SERBP1 affected the prognosis results of several cancers (Table 1). Specifically, we presented the survival curves of kidney renal clear cell carcinoma (KIRC), liver hepatocellular carcinoma (LIHC), and sarcoma (SARC) with most significant prognostic results (Figs. 1B–1D). The survival plots for the other 18 cancer types were also shown in supplementary figure (Fig. S1). In Fig. 1A, we observed that expression levels of SERBP1 in these three cancers were not significantly different, which was thought to be attributed to the complex factors affecting prognosis and the distinct functions of SERBP1 in these cancers. To detect the expression levels of SERBP1 in other cancer cell lines, we performed RT-qPCR experiment and found that SERBP1 was extensively expressed in multiple colorectal cancer cells, HEK293T cells (Fig. 1E), and hepatocellular carcinoma cells (Fig. 1F), while their expression levels were not uniform among these cell lines. In summary, these results indicate that SERBP1 is dysregulated and has the ability to affect the progression of multiple cancers.

SERBP1 promotes the apoptosis of HeLa cells

To confirm the function of SERBP1, we constructed a functional cell model. We predicted that the overexpression of SERBP1 may modulate the apoptosis levels of cancer cells. We succeeded in establishing transfected cells overexpressing SERBP1 (SERBP1-OE) fused to GFP fluorescent tag and flag tag Western blotting and RT-qPCR demonstrated this effect in both HeLa and HEK293T cells after 48 h (Figs. 2A, S2B, Table S2). Then, a flow cytometric experiment was conducted to explore the cellular functions of SERBP1. According to the experimental results, SERBP1-OE HeLa cells were proven to induce cell apoptosis compared with the Ctrl cells (Figs. 2B and 2C), whereas SERBP1-OE significantly reduced cell apoptosis level in HEK293T cells (Figs. S2C and S2D), and we found the protein level of caspase 3 was decreased by SERBP1-OE in HEK293T cells (Fig. S2B). According to these experimental results, we focused on the functions of the SERBP1-OE-induced genes in these cells.

Profiling the expression of SERBP1-regulated genes

HeLa cells were used to examine the responses in SERBP1-OE and empty vector-transduced cells. According to the RT-PCR results, the efficiency of SERBP1-OE was high with about 150-fold change (Fig. 2A). Then, cDNA libraries of samples were generated from control and SERBP1-OE cells to systematically investigate SERBP1-mediated transcriptional and post-transcriptional regulation. The RNA-seq libraries of four samples were sequenced on the Illumina HiSeq2000 platform to produce high-quality transcriptome data.

According to the sequencing data, a total of 92.0 ± 4.4 M raw reads were generated. Then, after the removal of the adaptors of the autogenous and low-quality bases and reads, the final clean reads were 89.3 ± 4.4 M. We compared the available sequencing data to the reference genome and obtained a total of 79.0 ± 4.0 M reads. The uniquely mapped reads were 75.1 ± 3.7 M, and were used for our subsequent analysis to ensure the reliability of the results.

We used expected FPKM to assess the gene expression patterns of each individual gene. There were 60,498 genes in the genome, and at least 12,611 genes were detected with FPKM > 1 in all conditions. SERBP1 also showed a significantly higher FPKM value in SERBP1-OE samples (Fig. 3A). The Pearson’s correlation coefficients based on FPKM values of all genes were used to calculate the correlation matrix, which was presented in the heatmap. The heatmap demonstrated that SERBP1-OE and Ctrl samples were clearly separated (Fig. 3B), while the correlation coefficients among the four samples were consistently high (R > 0.99), indicating that SERBP1-OE did not globally alter the transcriptome profile and could just affect the expression levels of a set of genes.

SERBP1 regulates HeLa cell transcription

We next explored SERBP1-regulated transcription by predicting differentially expressed genes (DEGs). We used edge R, which could identify DEGs between two or more samples, to identify the DEGs between SERBP1-OE and control cells with the criteria 5% FDR and fold change (FC) ≥2 or ≤0.5 (Robinson, McCarthy & Smyth, 2010). Among the 28,561 genes, we identified 249 upregulated DEGs and 245 downregulated DEGs (Fig. 3C, Table S3), and we found that SERBP1 could certainly regulate gene transcription. Additionally, the RNA-seq samples showed a high consistency in SERBP1 mediation and transcription. Heatmap analysis revealed the differences in gene expression patterns under different experimental conditions (Fig. 3D).

SERBP1 regulates gene expression in multiple cancer-related functions

A total of 494 genes showing significant differences between SERBP1-OE and Ctrl samples were identified. Then, functional annotation analyses were conducted to determine the potential biological roles of these DEGs using GO and KEGG databases. Then, we used the cutoff criterion to enrich the terms and identified 17 GO terms from the upregulated DEGs and 30 GO terms from the downregulated DEGs. According to the biological process from GO analysis, the upregulated genes in SERBP1-OE cells were only enriched in multicellular organismal development (Fig. 3E). In contrast, the downregulated genes in SERBP1-OE cells were mainly enriched in synaptic transmission, platelet activation, positive regulation of cell proliferation, blood coagulation, signal transduction, negative regulation of cell proliferation, positive regulation of transcription from the RNA polymerase II promoter, and apoptotic process (Fig. 3F), which was consistent with the functions of SERBP1 in cancer and in the regulation of proliferation and apoptosis. Based on the KEGG analysis, we also obtained the significantly enriched KEGG pathways of DEGs. Among the upregulated genes, we identified enrichment in protein digestion and absorption, pancreatic secretion, glycine, serine, and threonine metabolism (Fig. S3A). And, among the downregulated genes, we identified enrichment in tight junction, staphylococcus aureus infection, and cholinergic synapse (Fig. S3B).

Identification of SERBP1-regulated ASEs

After we obtained the results of the functional terms and pathways of DEGs associated with SERBP1 overexpression, we next aimed to determine the role of SERBP1 in AS regulation. We identified SERBP1-dependent ASEs using transcriptome sequencing in HeLa cells. First, we mapped each sample from the RNA-seq data to a unique position on the reference genome for AS analysis. The results showed that the annotated exons contained 367,321 genes, representing over 60% of the samples. We used TopHat2 to analyze splice junctions and found that the novel splice junctions and the known junctions in all the RNA-seq samples were 165,725 and 375,094, respectively. Then, we used ABLas software to investigate the overall changes in AS, which we had previously analyzed in TopHat2. We counted the detected annotated variable splicing events (known ASEs) in the genome in each sample. A total of 21,373 known ASEs were detected among all known ASEs from all samples, and 60,724 novel ASEs among all novel ASEs. Altogether we detected all 10 types of AS.

To analyze the project that was subjected to experimental processing conditions of SERBP1-associated RASEs, we compared the variable splicing level of each splicing type of each gene between the two samples using Student’s t-test with P-value < 0.05. In total, 321 upregulated RASEs and 365 downregulated RASEs were identified. In addition, the overall RASEs contained 194 intron retention (IR) RASEs, and the rest included 136 A5SS, 118 A3SS, 66 cassette exon, 80 ES, 28 MXE, 31 5pMXE, 18 A3SS&ES, 9 A5SS&ES, and 6 3pMXE RASEs (Fig. 4A). The data above indicated that SERBP1 regulates ASEs in HeLa cells. We then performed hierarchical clustering heatmap analysis of all significant RASE ratios, and found a clear separation between SERBP1-OE and control samples and very consistent values within replicates (Fig. 4B).

To determine whether AS events could be attributed to a simple increase in transcription, we analyzed SERBP1 and its regulation of the expression of downstream genes and investigated whether SERBP1-regulated AS had the same effect. As a result, we found that only four genes overlapped between the DEGs and the regulated splicing genes (Fig. 4C).

Furthermore, we studied the potential biological roles of SERBP1-regulated AS genes. The functional analysis of these genes included GO categories of 428 genes and KEGG pathways of 529 genes. With the criterion of P-value < 0.05, we found that the SERBP1-regulated AS genes were enriched in 57 GO terms in cellular component, 45 GO terms in molecular function, and 109 GO terms in biological processes. The highly enriched biological processes were protein ubiquitination involved in ubiquitin-dependent protein catabolic process, pyruvate metabolic process, fatty acid metabolic process, and negative regulation of DAN-dependent transcription (Fig. 4C). These results indicated that the overexpression of SERBP1 prompted cell apoptosis, especially in cancer cells. There was no overlap between the DEGs and RAS gene (RASG) enrichment. The KEGG pathway analysis identified 167 pathways, and the highly enriched pathways were ubiquitin-mediated proteolysis, bacterial invasion of epithelial cells, focal adhesion and regulation of actin cytoskeleton (Fig. 4D). These pathways included pathways in cancer, EGFR tyrosine kinase inhibitor resistance, and the Ras signaling pathway. There was also no overlap between the DEGs and RASG enrichment. Genes involved in each GO term or KEGG pathway are provided.

Validation of SERBP1-regulated AS of metabolic process genes in HeLa cells

With the analysis of the GO biological process terms of SERBP1-regulated AS genes (RASGs), we used RT-qPCR to validate the genes associated with the metabolic process in RASGs. Among the GO terms, we found that pyruvate metabolic process and fatty metabolic process were highly enriched in RASGs. The pyruvate metabolic process was associated with genes such as PDP1, SLC16A3, ME3, and PC, and the fatty acid metabolic process was associated with LPIN3, CROT, ACSF3, ALKBH7, and SLC27A1. According to these results, the RASEs from ME3, LPIN3, and CROT genes were downregulated when SERBP1 was upregulated, and the RASEs from PDP1, SLC27A1, and ALKBH7 genes were upregulated when SERBP1 was upregulated. Two splicing events of PC genes were up- and downregulated when SERBP1 expression was changed, respectively. We then validated these AS events after SERBP1-OE in HeLa cells. After performing RT-qPCR experiment, we found that PDP1, ME3, SLC27A1, LPIN3, and ALKBH7 showed consistent AS ratio changes with RNA-seq results (Figs. 5A–5E), validating the AS regulation of SERBP1 in HeLa cells.