RNA-binding protein CELF6 modulates transcription and splicing levels of genes associated with tumorigenesis in lung cancer A549 cells

CELF6 expression in tumor samples

To have a global view of CELF6 expression pattern and its association with survival time in lung cancer, we used two web-based tools, including GEPIA2 (Tang et al., 2019) and KM-plotter (Gyorffy, 2021) to perform the analysis of gene expression and cancer survival using public datasets, including the cancer genome atlas (TCGA). The abbreviations of tumor types could be found at the website of GEPIA2.

The culture of A549 cells and plasmid transfections

The identified human NSCLC cell line A549 cells with STR analysis were obtained from Procell (CL-0016, Wuhan, Hubei, China). Before cell culture, the A549 cells were detected to prevent mycoplasma contamination by the provider and our laboratory.

Then the A549 cells were cultured under traditional cell culture conditions which have been described detailedly in a previous study (Ma et al., 2020). CELF6 sequences were embedded in pIRES-hrGFP-1a plasmid (VT1056; YouBio, Shanghai, China). CELF6-OE and empty plasmids were transfected into A549 cells using Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s protocol. After 48 h cultivation, transfected cells were harvested for following experiments.

Validation of CELF6 over-expression in A549 cells

We performed reverse transcription and quantitative polymerase chain reaction (RT-qPCR) and Western blot (WB) experiments to assess the expression level of CELF6 in OE and NC cells. Complementary DNA (cDNA) synthesis was done by standard procedures and RT-qPCR was performed on the Bio-Rad S1000 with Hieff™ qPCR SYBR® Green Master Mix (Low Rox Plus; YEASEN, Shanghai, China). Primer sequences for RT-qPCR were shown in Additional File 1. The relative level of each transcript was normalized to GAPDH mRNA level (internal control) using 2^−ΔΔCT method (Livak & Schmittgen, 2001). Raw and processed data of RT-qPCR results were shown in Additional File 2. Two-way ANOVA method was used to calculate statistical significance.

Western blot to assess CELF6 protein level

A549 cells were treated by a published protocol (Ma et al., 2020). The treated samples were then boiled for 10 min in water with 1× SDS buffer and separated on 10% SDS-PAGE. Then samples were treated with TBST buffer (5% non-fat milk power and 20 mM Tris-buffered saline and 0.1% Tween-20) for 1 h at room temperature. Membranes were then incubated with primary antibody: FLAG antibody (1:1,000 dilution; polyclonal antibody; cat. no. 2368S; CST), GAPDH (1:2,000; ABclonal; Wuhan, Hubei, China) and then incubated with HRP-conjugated secondary antibody. Bound secondary antibody (anti-mouse or anti-rabbit 1:10,000) (Abcam) was detected using the enhanced chemiluminescence (ECL) reagent (Bio-Rad, Hercules, CA, USA, 170506). The uncropped gels for WB were presented in Fig. S1.

RNA-seq preparation and data analysis

After cell culture, we extracted polyadenylated RNAs and constructed RNA-seq library according to published protocols (Song Ba et al., 2021). Briefly, we used TRIZOL (Ambion) method to extract total RNA and removed genomic DNA using RQ1 DNase (RNase free; Promega, Madison, WI, USA). After checking the quality and quantity of purified RNA, we used 10 μg total RNA to prepare a directional RNA-seq library by capturing polyadenylated mRNAs. Three biological replicates were prepared for both CELF6-OE and NC groups. The libraries were prepared according to the manufacturer’s instructions for high-throughput sequencing. Finally, we got 151-bp pair-end sequences for each sample. Low quality sequencing reads (containing N bases and quality score less than 20) and adaptor sequences from raw sequences were removed using FASTX-Toolkit (Version 0.0.13). Too short reads (<16 nt) were dropped. Next, quality-filtered reads were aligned to the GRCH38 genome by TopHat2 (Kim et al., 2013), with no more than four mismatches allowed. Reads with unique genomic locations were used for gene reads number counting and FPKM calculation (fragments per kilobase of transcript per million fragments mapped) (Trapnell et al., 2010). We calculated the replication and experimental power using R package RNASeqPower (Therneau, Hart & Kocher, 2022). The final power value 0.99 was used to detect a 2-fold expression. Then, the R Bioconductor package edgeR (Robinson, McCarthy & Smyth, 2010) was utilized to identify DEGs with the criteria p-value < 0.01 and fold change >1.5 or <2/3.

CELF6-regulated alternative splicing events (RASEs) analysis

The AS events (ASEs) and CELF6-regulated AS events (RASEs) between CELF6-OE and control were analyzed using the ABLas pipeline (Xia et al., 2017b). Ten types of ASEs were analyzed, including alternative 5′ splice site (A5SS), alternative 3′ splice site (A3SS), cassette exon, exon skipping (ES), A5SS&ES, A3SS&ES, intron retention (IR), mutually exclusive exons (MXE), mutually exclusive 5′UTRs (5pMXE), and mutually exclusive 3′UTRs (3pMXE). To assess CELF6-regulated ASEs (RASEs), Student’s t-test was performed to calculate the statistical significance of the ratio difference between CELF6-OE and control samples. p-value < 0.05 and ratio difference >0.2 were considered criteria for RASEs.

Functional enrichment analysis

To explore enriched functions of DEGs and RASGs, we used KOBAS 2.0 server (Xie et al., 2011) to identify enriched Gene Ontology (GO) terms and KEGG pathways. Hypergeometric test and Benjamini-Hochberg FDR were used to calculate the significant p-values and FDR values for each term or pathway.

Reverse transcription qPCR validation of DEGs and RASEs

We also used RT-qPCR experiment to elucidate the DEG and RASG validity of the RNA-seq data. The details of DEG and RASG primers were presented in Additional File 1. For RASEs, we used a method to design specific PCR primers covering the junction boundary of constitutive exon and alternative exons (Liu et al., 2020). Total RNA of the cultured CELF6-OE and control A549 cells was used for validation. The relative expression levels of all the DEGs and RASGs were normalized to the internal control GAPDH. Raw and processed data of RT-qPCR results were shown in Additional File 2. Two-way ANOVA method was used to calculate statistical significance.

Raw sequencing data

The raw sequencing data discussed in this manuscript are available under GEO series accession number GSE185865.

CELF6 expression pattern in tumor samples from TCGA

To investigate CELF6 expression pattern and its influence on prognosis in lung cancer, we used GEPIA2 (Tang et al., 2019) and KM-plotter (Gyorffy, 2021) to analyze gene expression level and survival time using public datasets. CELF6 had significantly higher expression levels in normal tissue compared with in tumor tissue in 16 out of the 31 cancer types, whereas significantly higher expression levels of CELF6 in tumor tissue were only observed in LAML and PCPG (Figs. 1A and 1B). We found that the expression of CELF6 was significantly down-regulated (FDR < 0.05) in lung tumor samples (LUAD and LUSC) (Fig. 1B). We then detected its expression changes in disease stages of LUAD and LUSC patients. The result showed significant differences among the four stages (Fig. 1C), although the overall expression levels of CELF6 were very low (median TPM < 1). Furthermore, KM-plotter server was used to comprehensively explore the association between CELF6 expression and survival rate of LUAD patients. The results of log-rank test demonstrated that higher CELF6 level is associated with the longer overall survival (OS) time in LUAD patients (Fig.1B).

CELF6 over-expression experiment was conducted in A549 cells

To explore the functions of CELF6 in lung cancer, a cell model was constructed by using cell culture and transfecting pIRES-hrGFP-1a-CELF6 into A549 cells to enhance CELF6 expression level. In order to assess the CELF6 over-expression (CELF6-OE) result in A549 cells, we detected CELF6 mRNA expression by RT‑qPCR and CELF6 protein expression by Western blot. The mRNA level of CELF6 was elevated significantly in transfected cells (p < 0.001, two-tailed t-test, Fig. 2A). The protein level of CELF6 was also elevated in transfected cells (Fig. 2B). These results suggest that CELF6-OE A549 cell models were successfully constructed. After collecting cells, cDNA libraries on CELF6-OE and control cells (three biological replicates per sample) were constructed for RNA‑seq experiment. The expression level of CELF6 was up-regulated after over-expression based on RNA-seq (p < 0.001, two-tailed t-test, Fig. 2C). Principal component analysis (PCA) based on expression values of all expressed genes showed that CELF6-OE was the major factor globally influencing the gene expression pattern (Fig. 2D).

CELF6 over-expression resulted in some transcriptional differences

To investigate CELF6-mediated transcriptional regulation, we performed a comprehensive analysis of the RNA‑seq dataset. The sequencing result yielded about 83 million paired‑end raw reads per sample. After initial quality control steps, about 81 million high‑quality reads were generated per sample. We then found that about 97% of the reads from each sample could be aligned to the GRCh38 genome using TopHat2, and about 94% of the aligned reads from each sample had a uniquely aligned genomic location in the GRCh38 genome. We used uniquely mapped reads to calculate the normalized transcription level as FPKM, and detected 28,121 expressed genes from the six RNA‑seq datasets.

Based on the expression level of all detected genes from the CELF6-OE and control samples, we used strict criteria, absolute fold change ≥1.5 or ≤2/3 and FDR < 0.05 together with edgeR package, to identify differentially expressed genes (DEGs) that were potentially regulated by CELF6 at transcriptional level (Robinson, McCarthy & Smyth, 2010). In total, 417 up-regulated and 1,351 down-regulated DEGs were identified (Additional File 3). Volcano plot of all expressed genes displayed that CELF6-OE had a preference to repress gene expression in A549 cells (Fig. 3A). Hierarchical clustering heat map analysis of the DEG expression levels showed high regulatory consistency among the three replicates in both CELF6-OE and control samples (Fig. 3B). These results demonstrated that CELF6 has the potential to extensively regulate gene expression pattern in A549 cells.

CELF6-OE preferentially down-regulated the expression of large number of inflammation genes

To further decipher the potential biological functions of these DEGs, we performed GO and KEGG enrichment analyses for both up- and down-regulated DEGs. The top 10 GO biological process (BP) terms enriched by up-regulated DEGs were shown with transport, axon guidance, signal transduction, cell adhesion, small molecule metabolic process, and cellular protein metabolic process. Up-regulated DEGs enriched in KEGG analysis were shown with MAPK signaling pathway, Estrogen signaling pathway, ECM-receptor interaction, protein processing in endoplasmic reticulum and so on. (Fig. S2, details in Additional Files 4 and 5). We then focused our attention on down-regulated DEGs. Strikingly, we found that the down-regulated DEGs were mostly associated with type I interferon-mediated signaling pathway, cytokine-mediated signaling pathway, immune response, interferon-gamma-mediated signaling pathway, inflammatory response, and innate immune response (Fig. 3C). According to the parallel KEGG analysis results (details in Additional File 6), down-regulated DEGs were significantly enriched in cytokine-cytokine receptor interaction, TNF signaling pathway, and cell adhesion molecules (CAMs) (Fig. 3D, details in Additional File 7), which were associated with immune response, cell migration and proliferation (Chen & Mellman, 2017; Makrilia et al., 2009). We then selected DEGs in A549 cells for further validation by RT‑qPCR. A total of twelve DEGs, including HSPA6, HSPA1B, TNFSF10, CCL5, JUNB, BIRC3, MLKL, PIK3R2, CCL20, STAT1, MYD88, and CFS1, were selected. The results showed high consistency with the RNA‑seq data (Fig. 3E). These results together suggest that CELF6 regulates DEGs associated with cancer immunity and progression.

CELF6 modulates AS events of apoptosis-related genes in A549 cells

One definite aim of this study is to gain insight into the alternative splicing regulatory role of CELF6. Transcriptome data were used to examine regulated alternative splicing events (RASEs) by CELF6-OE in A549 cells. Among the uniquely aligned reads from CELF6-OE and control cells, about 38.94–40.41% were splicing junction reads. Of the 367,321 annotated exons from the human genome, 249,163 were detected. Using TopHat2 software (Kim et al., 2013), we detected 163,082 annotated and 162,218 novel splicing junctions, accounting for 49.8% of the total junctions (details in Additional File 8). Using the self-developed ABLas software (Xia et al., 2017a), we conducted overall analysis of and classified splicing methods on the results of splice junctions detected by TopHat2, and conducted statistics on various variable splicing events. A total of 17,691 known RASEs (annotated by reference genome) were identified, accounting for 26.3% of the total. The results showed that 49,565 novel splicing events (excluding intron retention) were identified, accounting for 73.7% of the total (details in Additional File 8). These results suggested that splicing profile is more complex than previously hypothesized. For each RASE, we calculated the differential splicing patterns between CELF6-OE and control samples. We found that 462 significant RASEs were identified (AS ratio difference ≥0.2 and P‑value ≤0.05). In detail, these RASEs included 47 3pMXE, 58 5pMXE, 288 A3SS, 33 A3SS&ES, 34 A5SS&ES, 303 A5SS,197 ES,55 MXE, 132 cassette exon and 362 intron retention events (Fig. 4A, details in Additional File 9). Hierarchical clustering heat map of all significant RASEs based on PSI showed high consistency with the CELF6-mediated post-transcription in both data sets (Fig. 4B). GO functional enrichment analysis demonstrated that the alternative spliced genes were enriched in translation, cellular protein metabolic process, gene expression, mRNA metabolic process, RNA metabolic process, cellular lipid metabolic process, and translational initiation (Fig. 4C). Enriched KEGG pathways (p < 0.05) included those involved in apoptosis-multiple species, p53 signaling pathway, platinum drug resistance, and apoptosis (Fig. 4D). The data suggested that CELF6 globally regulated ASEs in A549 cells. Excluding the changes in ASEs attributed to transcriptional regulation. Genes whose expression levels and AS were both regulated by CELF6 were also examined, and 85 genes were shared between the DEGs and RASGs (Fig. 4E). In addition, the differential splicing genes are enriched in p53 signaling pathway and apoptosis, which are closely related to tumorigenesis of lung cancer. From these differential genes, we screened out several genes associated with the occurrence of lung cancer. We chose CD47, CD44, FN1, CASP8, TP53, BAX, and SEC31A to identify CELF6-regulated alternative splicing events (Fig. S3). These results suggested that CELF6-OE markedly regulated AS of apoptotic-related genes.

Validation of CELF6 RASEs associated with apoptosis in A549 Cells

To validate CELF6-regulated AS pattern of apoptotic genes, RT-qPCR experiments were conducted. The results showed that five AS events out of seven RASGs related to apoptosis and ECM receptor interaction showed high agreement with the transcriptome analysis results of RNA-seq data (Figs. 5A, 5B and S3). Two important ASEs were originated from TP53 and CD44 (Figs. 5A and 5B), which have been well studied as key genes in cell proliferation in cancers.