DDX17 modulates the expression and alternative splicing of genes involved in apoptosis and proliferation in lung adenocarcinoma cells

Cell culture and transfection

Human A549 lung adenocarcinoma cells were obtained from the China Center for Type Culture Collection (CCTCC; Wuhan, China). The identify of cell line was authenticated (STR, Isoenzymatic Assay, PCR). We cultured the A549 cells with 10% fetal bovine serum (FBS) in minimal essential medium (MEM). The cells were maintained at 37 °C under 5% CO2 and were treated with 100 µg/ml streptomycin and 100 U/ml penicillin. Then the siRNA was transfected into the A549 cells with Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA) per the manufacturer’s instructions. After 48 h of transfection, cells were harvested for reverse-transcription quantitative PCR (RT–qPCR) analysis.

Assessment of gene expression

Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as an internal control. We synthesized the cDNA using standard methods, and we performed the RT-qPCR analysis on the Bio-Rad S1000 using the Bestar SYBR Green RT-PCR Master Mix (DBI Bioscience, Shanghai, China). Detailed primer sequences are listed in Table 1. Using GAPDH as a normalization factor, transcript expression was normalized and calculated by the 2-ΔΔCT method (Livak & Schmittgen, 2001). Statistical comparisons were made using paired Student’s t test in GraphPad Prism 7.0 software (Graphpad Prism software; GraphPad, San Diego, CA, USA).

Cell proliferation assay

Cell proliferation was examined using a Cell Counting Kit-8 (CCK-8; Dojindo, Co., Shanghai, China). In brief, A549 cells were seeded at a density of 6000 cells/well in 96-well culture plates. Cell-free wells were used as blank controls (blank), and cells treated with an equal volume of phosphate-buffered saline (PBS) were used as controls (control). After 24, 48 and 72 h of siRNA transfection, CCK-8 solution (20 µl) was added to the culture medium and incubated for an additional 3 h. Finally, the OD450 values were measured with a PerkinElmer EnVision plate reader. The cell proliferation rate was calculated according to the following formula: proliferation rate = (experimental OD value − blank OD value)/(control OD value − blank OD value) × 100%.

Annexin V apoptosis assay

An Annexin V-fluorescein isothiocyanate (FITC) and propidium iodide (PI) detection kit (Shanghai Yeasen Corporation, Shanghai, China) was used to detect apoptosis according to the manufacturer’s instructions. Specifically, A549 cells were cultured in six-well plates for 24 h and were then transfected with siRNA for 24, 48, or 72 h. Cells were treated with an equal volume of PBS as a control. The treated and control cells were mixed gently, and 5 µl of Annexin V-FITC reagent and 10 µl of 20 µg/ml PI reagent were then added and incubated at room temperature for 10–15 min in the dark. Following incubation, 400 µl of PBS was added. Flow cytometric analysis was used to detect apoptosis within 1 h. Flow cytometric analysis was used to detect apoptosis within 1 h. Cells positive for Annexin V-FITC or negative for PI fluorescence were considered apoptotic cells.

Acquisition of RNA-seq data from TCGA

For RNA-seq, LUAD gene-level expression data (526 tumor samples and 59 normal samples) from TCGA were downloaded from the UCSC Xena hub (http://xena.ucsc.edu). Splice junction data in BED format for 518 LUAD tumor and 59 normal samples from TCGA were downloaded from the Genomic Data Commons (GDC) portal for identification of ASEs (Kahles et al., 2018).

RNA extraction and sequencing

A Total RNA Extractor TRIzol (Ambion, Life Technologies, Carlsbad, CA, USA) was used to extract total RNA. RQ1 DNase (Promega, Madison, WI, USA) was used to remove any remaining DNA from the total RNA. RNA quality and quantity were assessed by measuring the absorbance ratio at 260 nm/280 nm (A260/A280) using a Nanodrop One spectrophotometer (Thermo Scientific, San Jose, CA, USA). 1.5% Agarose gel electrophoresis was performed to ensure RNA integrity.

For each sample, 1 µg of total RNA was used as the input material for RNA-Seq library preparation using a Kapa Stranded mRNA-seq kit (KAPA Biosystems, Wilmington, MA, USA) according to manufacturer’s instructions. We purified the polyadenylated mRNAs using VAHTS mRNA capture Beads (N401-01) and digested in fragmentation buffer. The fragmented mRNA was converted into double-stranded cDNA which was then subjected to end repair, poly(A) tailing and ligation with a KAPA Unique Dual-Indexed Adapter Kit (KAPA Biosystems, Wilmington, MA, USA). After ligation, the cDNA was sheared to fragments (300–500 bp). Finally, the DNA fragments were amplified, purified, quantified, and stored at −80 °C until used for sequencing. As the second strand of cDNA is synthesized, dTTP is replaced with dUTP, which quenches the second strand during amplification. The libraries were prepared and were run on the Illumina Novaseq 6000 system to generate 150-bp paired-end reads following the manufacturer’s instructions.

Cleaning and alignment of raw RNA-seq data

To start, we discarded any raw read containing more than N bases. Then we used the ASTX-Toolkit (Version 0.0.13) to removed adaptors from raw sequencing reads and trimmed low-quality bases. Short reads containing fewer than 16 nt were discarded. Subsequently, clean reads were mapped to the GRCh38 genome by TopHat2, allowing no more than four mismatches (Kim et al., 2013). Based on uniquely mapped reads, read numbers and fragments per kilobase of transcript per million fragments mapped (FPKM) values were calculated for each gene (Trapnell et al., 2010).

Analysis of differentially expressed genes (DEGs)

The analysis of DEGs was performed with the R Bioconductor package edgeR (Robinson, McCarthy & Smyth, 2010). The cutoff criteria for identifying DEGs were set as a p value <0.01 and a fold change>1.5 or <2/3.

To increase the statistical power, experiments were performed in three biological replicates. The statistical power of this experimental design, calculated in RNASeqPower (https://doi.org/doi:10.18129/B9.bioc.RNASeqPower) is 1 (Table S1).

Analysis of AS

The ABLas pipeline was applied to define and quantify the ASEs and regulated ASEs (RASEs) between the samples according to previously described methods (Jin et al., 2017a; Xia et al., 2017). In brief, ten types of ASEs were identified by the ABLas pipeline on the basis of the splice junction reads: exon skipping (ES), alternative 5′ splice site selection (A5SS), alternative 3′ splice site selection (A3SS), intron retention (IR), mutually exclusive exons (MXE), mutually exclusive 5′UTRs (5pMXE), mutually exclusive 3′UTRs (3pMXE), cassette exon (CassetteExon), A3SS&ES and A5SS&ES.

For analysis of RBP-regulated ASEs, Student’s t test was applied to evaluate the significance of the alteration ratio of ASEs. A false discovery rate (FDR) of 0.05 was set as the significance cutoff.

Validation of DEGs and ASEs by RT–qPCR

To confirm the results of RNA-seq-based identification of DEGs in A549 cells, some DEGs were selected for verification by RT–qPCR. Information on the primers is shown in File S3. A549 cells transfected with DDX17 siRNA were harvested after 48 h for RT–qPCR analysis. cDNA was synthesized from RNA using M-MLV Reverse Transcriptase (Vazyme). RT–qPCR was carried out on a StepOne RealTime PCR System using HieffTM qPCR SYBR^® Green Master Mix (Low Rox Plus; Shanghai Yeasen Corporation, Shanghai, China). The thermal cycling conditions used for PCR were as follows: denaturation at 95 °C for 5 min, followed by 40 cycles of denaturation at 95 °C for 15 s and annealing and extension at 60 °C for 30 s. PCR amplification was performed in triplicate for each sample. The RNA expression levels of all genes were normalized to the expression level of GAPDH.

Functional enrichment analysis

The biological functions of DEGs were comprehensively detected by Gene Ontology (GO) term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis using KOBAS 2.0 server (Xie et al., 2011). Hypergeometric test was used to test enrichment for each term and false discovery rate (FDR) procedure was used to correct p values.

Knockdown of DDX17 inhibits the proliferation and promotes the apoptosis of tumor cells

A CCK-8 assay was performed to evaluate cell proliferation. Compared to transfection with the control vector, siRNA knockdown of DDX17 significantly inhibited the proliferation of A549 cells (Fig. 1A). In addition, apoptosis analysis showed that DDX17 knockdown led to a marked increase in the apoptosis rate in A549 cells (Fig. 1B).

DDX17 regulates the transcriptome in A549 cells

To examine DDX17-mediated transcriptional regulation in A549 cells, we constructed six cDNA libraries with three biological replicates each from DDX17 knockdown and control cells, which were subjected to RNA-seq analysis. DDX17 was knocked down by siRNA and quantified by RT–qPCR.

As shown in Fig. 2A, the expression of DDX17 was significantly reduced in A549 cells transfected with the targeting siRNA. Effective DDX17 knockdown was further confirmed by RNA-seq analysis and western blotting (Figs. 2B, 2C; Fig. S1). Principal component analysis (PCA) showed that the DDX17 siRNA group was obviously separated from the control group (Fig. 2D).

The DEGs, namely, 992 upregulated genes and 625 downregulated genes, are visualized in an MA plot (Fig. 2E). Heatmap analysis based on the expression patterns of the DEGs demonstrated clear separation between the DDX17 siRNA samples and control samples and high consistency among the three biological replicates (Fig. 2F). These results indicated that DDX17 markedly regulated the expression levels of many genes.

We performed GO analyses to further investigate the potential biological functions of these DEGs (Figs. 2G, 2H). Functional enrichment analyses indicated that genes upregulated in the DDX17 knockdown group were enriched mainly in biological processes, including regulation of transcription (DNA-dependent), transcription (DNA-dependent), cAMP catabolic process, extracellular matrix organization, and regulation of cell growth, whereas downregulated genes were related primarily to nucleosome assembly, cell cycle arrest, DNA repair, endoplasmic reticulum unfolded protein process, and cell–cell signaling.

DDX17 selectively regulates the expression of genes involved in the progression of cancer cells

C−X−C motif chemokine ligand 5 (CXCL5), C-C motif chemokine ligand 20 (CCL20), Growth differentiation factor 15 (GDF15), Stanniocalcin-2 (STC2) and RAD51 recombinase (RAD51) have been reported to play a role in promoting the proliferation, migration, and invasion of NSCLC cells (Wang et al., 2018b; Wang et al., 2016; Wang et al., 2018a; Na et al., 2015; Hu et al., 2019). As shown in Fig. 3 and Fig. S2, RNA-seq indicated that the mRNA expression of these genes was significantly reduced in the DDX17 knockdown group compared with the control group. To confirm the reliability of the RNA-seq data, RT–qPCR was performed. The qPCR and RNA-seq results exhibited a high degree of concordance.

Similarly, the expression of tumor suppressor genes (e.g., Ankyrin repeat domain 1 (ANKRD1); hexamethylene-bis-acetamide-inducible protein 1 (HEXIM1); insulin like growth factor binding protein 4 (IGFBP4); insulin like growth factor binding protein 3 (IGFBP3); ras homolog gene family, member B (RHOB); and P53 apoptosis effector related to PMP22 (PERP)) was significantly increased in the DDX17 knockdown group compared with the control group (Jiménez et al., 2017; Tan et al., 2016; Li et al., 2017b; Wang et al., 2017; Chen et al., 2016; Chen et al., 2011). These results suggested that DDX17 may regulate the apoptosis and proliferation of LUAD cells by altering the expression of these genes.

DDX17 selectively regulates AS of genes involved in cell proliferation, migration, invasion, and apoptosis

To investigate the impact of DDX17 knockdown on AS of genes in A549 cells, we evaluated splicing alterations in the RNA-seq data using the ABLas pipeline. A total of 2094 RASEs were detected (Fig. 4A). Among the nine types of ASEs, ES, A5SS, A3SS and CassetteExon were the main types detected. Hierarchical clustering of RASEs revealed obvious separation between the DDX17 siRNA samples and the control samples, with a high degree of consistency among the three biological replicates (Fig. 4B). To explore the possible association between the regulated alternatively spliced genes (RASGs) and DEGs, we examined the overlap between the two datasets (Fig. 4C). Both the expression level and AS of 71 genes were significantly affected by DDX17 knockdown.

GO functional enrichment analysis showed that RASGs were enriched mainly in biological processes related to the terms apoptotic process, DNA repair, cellular protein metabolic process, viral reproduction, positive regulation of translation, and positive regulation of apoptotic process (Fig. 4D). Among these RASGs, the following have previously been reported to mediate the proliferation, migration, invasion, and apoptosis of cancer cells: Ribosomal protein S27a (RPS27A); LSM1 homolog, mRNA degradation associated (LSM1); legumain (LGMN); inverted formin 2 (INF2); LDL receptor related protein 5 (LRP5); and filamin A (FLNA) (Yu et al., 2021; Watson et al., 2008; Roslan et al., 2019; Zhang & Lin, 2021; Qian et al., 2018b; Qian et al., 2018a). Significantly fewer ES events in LSM1, CassetteExon events in INF2, CassetteExon events in LGMN and A3SS events in LRP5 were detected by RNA-seq in the DDX17 knockdown group than in the control group (Figs. 4E, 4F; Fig. S3). Conversely, significantly more A5SS events in RPS27A and CassetteExon events in FLNA were detected in the DDX17 knockdown group than in the control group. These results were confirmed by RT–PCR. In particular, these findings indicate that DDX17 can regulate both the AS and mRNA expression of LGMN.

DDX17 regulates the transcriptome in LUAD in the TCGA database

To further confirm the transcriptional activity of DDX17, we analyzed gene expression data from TCGA datasets (including 526 LUAD samples and 59 normal lung samples) by using the UCSC Xena platform (https://xena.ucsc.edu/). We selected 20 samples with the highest DDX17 expression (all of which are cancer tissues) and 20 samples with the lowest DDX17 expression (including 19 cancer tissues and 1 normal lung tissue) to assess the transcriptional regulation of DDX17. A total of 5,243 genes were differentially regulated, among which 4,218 were upregulated and 1,025 genes were downregulated (Fig. 5A). The DEGs were visualized in a heatmap, which showed clear separation between the two groups (Fig. 5B). The Venn diagram of the DEGs showed that 67 genes were upregulated in LUAD samples with high DDX17 expression compared to those with low DDX17 expression, whereas these genes were downregulated in A549 cells after DDX17 knockdown (Fig. 5C). Similarly, 23 genes were downregulated in LUAD samples with high DDX17 expression compared to those with low DDX17 expression, whereas these genes were upregulated in A549 cells after DDX17 knockdown (Fig. 5D). Among these 23 genes, apoptosis-related genes such as IGFBP4, IGFBP3 and PERP were upregulated in LUAD samples with low DDX17 expression compared to those with high DDX17 expression in the TCGA database (Fig. 5E).

DDX17 regulates AS of genes in LUAD in the TCGA database

Finally, to investigate whether DDX17 influences AS of genes in LUAD patients, we analyzed BED format splicing data (including 518 LUAD samples and 59 normal lung samples) from the GDC portal (Kahles et al., 2018). We selected 20 samples with the highest DDX17 expression and 20 samples with the lowest DDX17 expression to assess ASEs regulated by DDX17. Similar to the observations in A549 cells, A5SS, A3SS and CassetteExon were the main types of ASEs regulated by DDX17 in LUAD patients (Fig. 6A). By analyzing the overlap between the RASEs in the TCGA database and the RASEs in siDDX17 cells, we found significant overlap (175 genes) between the RASGs in the LUAD data and the RNA-seq data (Fig. 6B). GO functional enrichment analysis indicated that these genes were enriched mainly in biological processes related to the terms RNA splicing; mRNA processing; phosphorylation; positive regulation of transcription (DNA-templated); mRNA splicing, via spliceosome; and proliferation (Fig. 6C). Then, these overlapping genes were subjected to KEGG pathway enrichment analysis and were found to be enriched mainly in biological processes related to bacterial invasion of epithelial cells, focal adhesion, Epstein −Barr virus infection, the mRNA surveillance pathway, and the proteasome (Fig. 6D). Among the overlapping genes, significantly fewer A5SS events in LSM1, CassetteExon events in INF2, A5SS events in FLNA, A5SS events in RPS27A, and A3SS events in LGMN were detected by RNA-seq in LUAD samples with low DDX17 expression than in those with high DDX17 expression, whereas A3SS events in LRP5 were more prevalent in LUAD samples with high DDX17 expression (Fig. 6E).