This article is mentioned in:

Introduction

RNA helicases are enzymes that unwind nucleic acids in an energy-dependent manner via hydrolysis of nucleotide triphosphate (1); these enzymes are implicated in processes involving RNA, including transcription, splicing, RNA export, ribosome biogenesis, mRNA translation and RNA decay (2–5). Moreover, novel functions of RNA helicases include their activity as transcriptional coregulators and regulators of post-translational modifications, and their ability to modulate cellular signaling pathways such as Wnt-ß-catenin, NF-κB and MAPK (6–9). Deregulation of these RNA helicase functions may contribute to tumor development and progression (10–12). Dysregulation of RNA helicase expression levels, its mutation status in tumors and its role in the regulation of molecules have been implicated in cancer (13). High-throughput sequencing technology, with its rapidly decreasing costs and increasing applications, allows researchers to investigate the biological function of RNA helicases in cancer, which may be the key to ‘unwind’ cancer in the future.

DEAD-box helicase 41 (DDX41) is a highly conserved protein considered essential for cell growth and viability (14), which is mutated in primary tumor and relapse samples from patients with acute myeloid leukemia (AML) (15). Polprasert et al (16) identified somatic DDX41 mutations in myeloid neoplasms that resulted in loss of tumor suppressor function due to altered pre-mRNA splicing and RNA processing. RNA-sequencing (seq) analysis on peripheral bold mononuclear cells of DDX41 mutation carriers with hematological malignancy (HM) revealed altered expression levels of genes involved in hemoglobin complex and innate immunity (17). The majority of germline mutations are frameshift mutations, suggesting loss of function, with DDX41 serving as a tumor suppressor, which may impact initiation, maintenance or progression of tumorigenesis (18). The identification of DDX41 mutations has improved understanding of the potential mechanisms underlying HM; however, the expression and mechanism by which DDX41 contributes to tumorigenesis is still unknown.

As an upstream DNA sensor of stimulator of interferon gene (STING), DDX41 may serve dual roles in various types of cancer. The preventive role of DDX41 in HM is associated with inhibition of proliferation and promotion of apoptosis and DDX41 negative regulates p21 at the translational level (19); in addition, DDX41 silencing promotes apoptosis of HCT116 cells in a p21-dependent manner (19). These studies suggest that the role of DDX41 may be tumor type-dependent. However, the precise function of DDX41, and how it influences the formation and development of tumors remains poorly understood.

DDX41 is also associated with the immune response. Recent studies revealed that DDX41 primarily senses the viral DNA/RNA hybrid and is required for the anti-retroviral innate immune response to murine leukemia and HIV in primary mouse macrophages and dendritic cells (DCs) (20,21). Zhang et al (22) identified DDX41 as an intracellular DNA sensor in myeloid DCs; silencing of DDX41 was shown to block the activation of STING-TANK-binding kinase 1 signaling and the transcription factors NF-κB and interferon regulatory factor 3, resulting in decreased production of type I interferon. As a pattern-recognition receptor, DDX41 also recognizes bacterial di-GMP and di-AMP to activate the immune response via the STING-dependent signaling pathway (23). Intravenous administration of c-di-GMP into mice efficiently activates DDX41-STING signaling, production of interferon and activation of natural killer (NK) cells, and also facilitates antigen-specific cytotoxic T cell activity, resulting in a significant antitumor effect in a lung metastasis mouse model using malignant melanoma cells (24,25). Thus, DDX41 may trigger the production of molecules involved in the immune response.

Given the involvement of helicases in RNA splicing (26), Polprasert et al (16) expressed epitope-tagged versions of wild-type (WT) and R525H-mutated DDX41 in 293 cells. Epitope pull-down coupled with peptide sequencing analysis revealed that spliceosome proteins constitute a DDX41-associated group, which includes pre-mRNA processing factor 8 and splicing factor 3B subunit 1 (16). Furthermore, DDX41 mutation in the DEAD domain (such as R525H) alters this interaction, especially for major components in the U2 and U5 spliceosomes (16). In addition, DDX41 has been reported to be involved in pre-RNA splicing and DDX41 silencing may be associated with increased exon skipping (16). DDX41 is specifically recruited to the catalytically active C complex, which performs the second step of splicing, in which the 5′ and 3′ exons are ligated and an intronic lariat is released (27,28). These findings indicated that DDX41 may lead to the impairment of RNA processing and pre-mRNA splicing, which may partially contribute to tumorigenesis.

The present study aimed to investigate the function of DDX41 in HeLa cells and cervical cancer. DDX41 was cloned and overexpressed in HeLa cells. RNA-sequencing (RNA-seq) analysis of the effect of DDX41 on the transcriptome and alternative splicing of the overexpressed and normal control cell samples was performed. These results may improve understanding of the biological role of DDX41 in cancer.

Materials and methods

DDX41 cloning and plasmid construction

Primer pairs used for Hot Fusion were designed using CE Design v1.04 (Vazyme Biotech Co., Ltd.) with gene-specific sequences and oligonucleotides encoding 3× FLAG tag were added before the TAG termination codon, along with a portion of the vector pIRES-hrGFP-1a sequence, each with a 17–30 bp overlap. The following primers were used: Forward, 5′-agcccgggcggatccgaattcATGGAGGAGTCGGAACCCG-3′; and reverse, 5′-gtcatccttgtagtcctcgagGAAGTCCATGGAGCTGTGGGC-3′. Non-capitalized letters indicate sequence for homologous recombination during cloning, whereas capitalized letters indicate sequence paired with DDX41 gene. Vector pIRES-hrGFP-1a was digested by EcoRI and XhoI (New England BioLabs, Inc.) at 37°C for 2–3 h. The enzyme-digested vector was then run on a 1.0% agarose gel and purified using a QIAquick Gel Extraction kit (cat. no. 28704; Qiagen, Inc.). Total RNA from HeLa cells was isolated using TRIzol® (Invitrogen; Thermo Fisher Scientific, Inc.). Then, 1 µg total RNA was transcribed into cDNA using a PrimeScript™ RT Reagent kit (Takara Bio, Inc.), according to the manufacturer's protocol. The thermocycling condition were as follows; 98°C for 30 sec; 28 cycles of 98°C for 10 sec and 60°C for 30 sec; final extension at 72°C for 1 min and holding at 4°C. Linearized vector digested by EcoRI and XhoI (New England BioLabs, Inc.) and the PCR insert (1,866 bp) were added to a PCR microtube and cloned using the ClonExpress® II One Step Cloning kit (Vazyme Biotech Co., Ltd.), according to the manufacturer's protocol. Plasmids were introduced into Escherichia coli by chemical transformation. Cells (5×106) were plated onto LB plates containing ampicillin and incubated overnight at 37°C. Colonies were screened using a TaKaRa Taq™ HS PCR kit (Takara Bio, Inc.) according to the manufacturer's protocol, with M13 primers (forward, 5′-TGTAAAACGACGGCCAGT-3′ and reverse, 5′-CAGGAAACAGCTATGACC-3′). The PCR insert was verified by Sanger sequencing.

Cell culture and transfection

The human HeLa cell line was purchased from the Institute of Biochemistry and Cell Biology, Chinese Academy of Sciences, Shanghai, China. Cells were seeded into a petri dish (100 mm) at a density of 1×105 cells/ml and cultured at 37°C with 5% CO2 in DMEM (Invitrogen; Thermo Fisher Scientific, Inc.) containing 10% fetal bovine serum (HyClone; GE Healthcare Life Sciences), penicillin (100 U/ml) and streptomycin (100 g/ml). Transfection of 1×105 HeLa cells with 0.5 µg DDX41 overexpression or empty plasmid was performed using Lipofectamine® 2000 (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. Transfected cells were harvested after 48 h for subsequent experimentation.

Assessment of DDX41 overexpression

GAPDH was used as an internal control. cDNA synthesis was performed according to standard procedure (29) and reverse transcription-quantitative (RT-q)PCR was performed on the Bio-Rad S1000 system with Bestar SYBR Green RT-PCR Master Mix kit according to the manufacturer's instructions. (DBI Bioscience; Shanghai Xinghan Biotechnology Co., Ltd.). The concentration of each transcript was then normalized to the GAPDH mRNA levels using the 2−ΔΔCq method (30). Comparisons were performed by one-way ANOVA using GraphPad Prism 8 software (GraphPad Software, Inc.).

Western blot analysis

In brief, for preparation of total cell lysates, normal and DDX41-overexpressing HeLa cells were lysed in RIPA buffer containing 50.0 mM Tris-HCl (pH, 7.4), 150.0 mM NaCl, 1.0% deoxycholate, 1.0% Triton X-100, 1.0 mM EDTA and 0.1% SDS. The samples were centrifuged (10,000 × g; 5 min; 4°C). The protein level was determined by Pierce™ BCA protein assay kit (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. Then, 50 µg protein was incubated for 10 min in boiling water with 1X SDS sample buffer, separated by SDS-PAGE on a 10% gel, then transferred to a PVDF membrane (EMD Millipore, Inc.). The membranes were blocked with 5% skimmed milk [in buffer containing 10 mM Tris (pH, 8.0), 150 mM NaCl and 0.05% Tween-20] for 1 h. DDX41 was detected using monoclonal Flag antibody (1:2,000; cat. no. F7425; Sigma-Aldrich; Merck KGaA) with actin as a control (1:2,000; cat. no. AC038; ABclonal Biotech Co., Ltd.) diluted in TBST (0.05% Tween-20) at 4°C overnight. The blots were incubated with horseradish peroxidase-conjugated secondary antibody (goat anti-rabbit IgG; 1:10,000; cat. no. ab7090; Abcam) for 1 h at room temperature. Bands were detected using enhanced chemiluminescence reagent (Thermo Fisher Scientific, Inc.).

MTT assay

MTT assay was used to measure cell proliferation. Briefly, cells (1×104) were seeded in 96-well culture plates containing 200 µl cell growth medium. Upon reaching 70% confluence, the HeLa cells were transfected with DDX41 overexpression vector or control vector using Lipofectamine 2000 according to the manufacturer's protocol. The cells were then incubated at 37°C for 48 h. Subsequently, 0.025 ml MTT solution (5 mg/ml) was added to each well and cells were incubated for 4 h. Following centrifugation (5,000 × g; 5 min; room temperature), the supernatant was removed from each well. The colored formazan crystals produced from MTT in each well were dissolved in 0.15 ml DMSO and the optical density values were measured at 490 nm.

Flow cytometric analysis of cell apoptosis

HeLa cells (5×104) were seeded in 24-well culture plates. At 70% confluence, the cells were transfected with the DDX41 overexpression or control vector using Lipofectamine 2000 (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol. The cells were then incubated at 37°C for 48 h and viable cells were harvested and washed twice with PBS. The cell pellets were resuspended in 400 µl ice-cold 1X binding buffer at ~1×106 cells/ml, then incubated with 10 µl FITC-conjugated Annexin V/7-AAD (Beijing 4A Biotech Co., Ltd.) for 10 min in the dark at room temperature. Samples were analyzed within 1 h of staining using a A00-1-1102 CytoFLEX flow cytometer (Beckman Coulter, Inc.) and Flow-Jo software (v8.8.6; Tree Star, Inc.).

Library preparation and sequencing

Total RNA was extracted from HeLa cells using TRIzol. The RNA was further purified with two phenol-chloroform treatments and then treated with RQ1 DNase (Promega Corporation) at 37°C for 10 min to remove DNA. The quality and quantity of the purified RNA were determined by measuring the absorbance at 260/280 nm using Smartspec Plus (Bio-Rad Laboratories, Inc.). The integrity of RNA was further verified by 1.5% agarose gel electrophoresis. For each sample, 1 µg total RNA was used for RNA-seq library preparation using the VAHTS Stranded mRNA-seq Library Prep kit (cat. no. NR612-01; Vazyme Biotech Co., Ltd.). Polyadenylated mRNA was purified and fragmented, and then converted into double stranded cDNA. Following end repair and A tailing, DNA was ligated to VAHTS RNA Adapters (Vazyme Biotech Co., Ltd.). Purified ligation products corresponding to 200–500 bp were digested with heat-labile uracil-DNA glycosylase and the single stranded cDNA was amplified, purified, quantified and stored at −80°C before sequencing. Next, the 3M library was sequenced on an Illumina HiSeq X Ten platform in a 2 × 150 bp paired-end mode using HiSeq X Ten reagent kit v2.5 (cat. no. FC-501-2501; Illumina, Inc.). The fastq files were produced using the CASAVA pipeline v1.8.2 (Illumina, Inc.). Preliminary quality control analysis of the fastq files was performed using FASTQC software v 1.11.4 (fastqc.software.informer.com/).

Data processing and alignment

Raw reads containing >2-N bases were first discarded. Adaptors and low-quality bases were then trimmed from the raw sequencing reads using the FASTX-Toolkit (Version 0.0.13; hannonlab.cshl.edu/fastx_toolkit/). Short reads (<16 nucleotides in length) were also excluded. Clean reads were aligned to the GRCh38 genome by Tophat2 2.0.0 (ccb.jhu.edu/software/tophat/index.shtml) allowing four mismatches. Uniquely mapped reads were used to calculate the read number and fragments per kilobase of transcript per million fragments mapped for each gene.

Differentially expressed genes (DEGs) analysis

The R Bioconductor package edgeR 3.12 (31) (bioconductor.org/packages/release/bioc/html/edgeR.html) was utilized to screen out DEGs. A false discovery rate (FDR) <0.05 and fold change >1.5 or <0.7 were set as the cut-off criteria for identifying DEGs. In order to predict the gene function and calculate the functional category distribution frequency, Gene Ontology (GO) and enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were identified using the KOBAS 2.0 server (kobas.cbi.pku.edu.cn/) (32). Hypergeometric tests and the Benjamini-Hochberg FDR controlling procedure were used to define the enrichment of each pathway (corrected P-value <0.05).

Alternative splicing analysis

The alternative splicing events (ASEs) and regulated ASEs (RASEs) in the samples were defined and quantified using the ABLas pipeline, as described previously (33). In brief, detection of eight types of ASE was based on the splice junction reads. The eight types of ASE included cassette exon (CassetteExon), exon skipping (ES), mutually exclusive ES (MXE), alternative 5′ splice site (A5SS), alternative 3′ splice site (A3SS), MXE combined with an alternative 5′ promoter (5pMXE) and MXE combined with an alternative polyadenylation site (3pMXE). After detecting the ASEs in each RNA-seq sample, the paired t-test was chosen to calculate the P-value, with the alternative reads and model reads of samples as input data. The changed ratio of alternatively spliced reads and constitutively spliced reads between compared samples was defined as the RASE ratio. P<0.05 and RASE ratio >0.2 were set as the threshold for RASE detection.

RT-qPCR validation of DEGs and ASEs

In order to confirm the validity of the RNA-seq data, RT-qPCR was performed for selected DEGs and normalized to the reference gene GAPDH. The primer sequences are presented in Table SI. The same RNA samples for RNA-seq were used for qPCR. The PCR conditions consisted of denaturation at 95°C for 10 min, followed by 40 cycles of denaturing at 95°C for 15 sec and annealing and extension at 60°C for 1 min. A total of three independent repeats was performed for each sample. RT-qPCR assay was also used to analyze ASEs in HeLa cells. PCR primer pairs (Table SI) were designed to amplify the long and short splicing isoforms in the same reaction. In order to detect one of the alternative isoforms, one primer was designed for the alternative exon, and an opposing primer was designed for the constitutive exon; to detect the other alternative isoforms, a boundary-spanning primer was designed for the sequence encompassing the exon-exon junction, with the opposing primer in a constitutive exon. The changed ratio of AS events in RNA-seq was calculated using the formula: AS1 junction reads/AS1 junction reads + AS2 junction reads. The altered ratio of AS events in RT-qPCR was calculated using the formula: AS1 transcripts level/AS2 transcripts level.

Downloading data

The RNA-seq data of cervical cancer samples were downloaded from The Cancer Genome Atlas (TCGA) database (xenabrowser.net/datapages/) to analyze the expression levels of DDX41-regulated genes. Tumor suppressor genes were downloaded from TSGene 2.0 database (bioinfo.mc.vanderbilt.edu/TSGene/). The oncogenes were downloaded from an oncogene database (ongene.bioinfominzhao.org/index.html).

Gene expression profiling interactive analysis (GEPIA)

GEPIA was used to analyze RNA-seq expression data of 9,736 tumors and 8,587 normal samples from TCGA and Genotype-Tissue Expression (GTEx) projects, using a standard processing pipeline (gepia.cancer-pku.cn/). GEPIA provides customizable functions, such as tumor or normal differential expression analysis, profiling according to cancer types or pathological stages, patient survival analysis, similar gene detection, correlation analysis and dimensionality reduction analysis.

Immunofluorescence assay

HeLa cells were placed on glass slides using a cytospin (Shandong Junteng Medical Technology Co., Ltd.). Cells were fixed in 4% paraformaldehyde and permeabilized in 1% Triton X-100/PBS at room temperature for 15 min. Cells were incubated with primary antibody for DDX41 (1:100; cat. no. A6576; ABclonal) for 1 h at room temperature in a humidified chamber, washed with 0.1% Tween-20/PBS and incubated with species-appropriate FITV-conjugated secondary antibody (1:500; cat. no. AS011; ABclonal) for 1 h at room temperature. Slides were mounted with 100 ng/ml 4,6 diamidino-2-phenylindole Vectashield (Vector Laboratories, Inc.; Maravai Life Sciences) to stain nuclei at room temperature for 15 min and cells were visualized using an Olympus AX70 epifluorescence microscope (Olympus Corporation) with triple color filters and a 20× air objective. Numeric aperture varied automatically according to the light emitted. Images were captured with a Zeiss Axiocam camera (Zeiss GmbH).

Statistical analysis

Data are presented as the mean ± SD of ≥3 independent repeats. The statistical difference between groups was analyzed using one-way ANOVA with SPSS software (version 17.0; SPSS, Inc.). Patient survival analysis was conducted using the Kaplan-Meier method and analyzed using the log-rank test. The heatmap analysis compared the data obtained from DEG of RNA-seq (34) and sashimi plot representation of RNA-seq data was used in the Integrative Genomics Viewer (35). P<0.05 was considered to indicate a statistically significant difference.

Results

DDX41 inhibits proliferation and promotes apoptosis of HeLa cells

Previous studies have reported that in AML cell line xenografts in a mouse model, DDX41 knockdown accelerates tumor growth compared with control cells (16,48). In order to investigate the potential biological function of DDX41 in cervical cancer, DDX41-3×Flag plasmids was stably overexpressed in HeLa cells and compared with empty vector-transfected cells. The results showed significantly increased expression levels of DDX41, as assessed by western blotting and RT-qPCR (Fig. 1A). MTT assay was performed to determine cell proliferation; overexpression of DDX41 suppressed HeLa cell proliferation (Fig. 1B). In order to investigate the function of DDX41 in cell apoptosis, flow cytometry was performed. The results revealed that DDX41 overexpression significantly induced apoptosis of HeLa cells (Fig. 1C). By immunostaining with DDX41 antibody, DDX41 protein was detected both in nucleus and cytoplasm (Fig. S2B). Therefore, DDX41 may serve as a tumor suppressor in HeLa cells.

Figure 1. DDX41 inhibits proliferation and promotes apoptosis. (A) HeLa cells were transiently transfected with DDX41 OE or Con vectors followed by determination of Flag-tagged DDX41 protein levels by western blot analysis. The mRNA expression levels were measured by reverse transcription-quantitative PCR. (B) Cell proliferation was measured by MTT assay. Results are presented as the mean ± SD. (C) Effects of DDX41 on HeLa cell apoptosis were determined by flow cytometric analysis (n=3). *P<0.05, **P<0.01 and ***P<0.001. DDX41, DEAD-box helicase 41; Con, control; OE, overexpression; 7-AAD, 7-amino actinomycin D; OD, optical density; PE, phycoerythrin.

RNA-seq analysis of the impact of DDX41 overexpression on the HeLa cell transcriptome

In order to analyze the gene expression profile of DDX41, RNA-seq was performed for 83±4 million total raw reads per sample (Table SII). Filtered reads were aligned to the human GRCH38 genome by Tophat2 (36) and expression patterns of 28,914 known annotated genes were detected and characterized (Table SIII). edgeR (31) and (cutoff, fold change ≥1.5 or ≤0.7; P<0.01) was used to identify the DEGs between the DDX41-overexpressing cells and controls. An M-versus-A plot was constructed to display the DEGs that were associated with DDX41 overexpression and identified 504 up- and 455 downregulated genes (Fig. 2A; Table SIV). Heatmap analysis of the expression patterns of all DEGs in each sample was performed; the results showed high consistency for DDX41-mediated transcription in all datasets (Fig. 2B).

Figure 2. Differential gene expression and ASEs in response to DDX41 overexpression. (A) Detection of DDX41-regulated genes on the M-versus-A plot. Red, upregulated (FC≥1.5, P-value<0.01); blue, downregulated (FC≤0.7, P-value<0.01). (B) Hierarchical clustering of the expression levels of 959 DEGs regulated by DDX41 overexpression in HeLa cells expressing either the control or DDX41 plasmid. (C) Classification of different alternative spicing types regulated by DDX41 protein. FC, fold change; RASE, regulated alternative splicing events. 3pMXE, mutually exclusive 3′ untranslated regions; 5pMXE, mutually exclusive 5′ untranslated regions; A3SS, alternative 3′ splice site; A5SS, alternative 5′ splice site; ES, exon skipping; IntronR, intron retention.

The ABLas software tool (under permission) was used to analyze global changes in DDX41-mediated ASEs in HeLa cells. A total of 20,532 annotated ASEs were detected when comparing these uniquely mapped reads to the reference genome, as well as 60,654 novel ASEs (Table SV). A stringent cutoff of P≤0.05 was set to identify high-confidence RASEs and 707 RASEs were identified (Tables SVI and SVII). DDX41-regulated ASEs included 24 known intron retention events, 63 CassetteExons, 83 ES, 89 A3SS and 131 A5SS. The other event types included 11 3pMXE, 26 5pMXE, 13 A3SS + ES, 9 A5SS + ES and 18 MXE (Fig. 2C). These results suggested that DDX41 primarily regulated A5SS, A3SS and ES events in HeLa cells.

DDX41 regulates transcription and alternative splicing of certain tumorigenesis-associated genes

In order to test the hypothesis that DDX41 may regulate tumorigenesis-associated gene transcription, a gene dataset was constructed by downloading tumor suppressor genes from TSGene 2.0 (bioinfo.mc.vanderbilt.edu/TSGene/) (37) and oncogenes from an oncogene database (ongene.bioinfominzhao.org/index.html). The tumorigenesis-associated gene dataset was overlapped with DDX41-regulated genes (RNA-seq-identified DEGs). Heatmap analysis of the expression patterns of the overlapped genes was performed and showed a high consistency for DDX41-regulated gene expression in each sample (Fig. 3A). Certain oncogenes were downregulated by DDX41, such as LCN2, TP63, KIT, MYB, CDK5R2 and TNFRSF1B, whereas others were upregulated, including PANO1, GLI1, HOMER2 and RND3 (Fig. 3A). In order to assess the association between DDX41-regulated expression levels of oncogenes and tumor suppressors, expression data from CESC and control samples were downloaded (xenabrowser.net/datapages/) and analyzed. Expression levels of five oncogenes, TP63, LCN2, MYB, CSF2 and CYP24A1, were upregulated (Fig. 3B), whereas the expression levels of tumor suppressors GLI1 and PANO1 were downregulated in cancer samples (Fig. 3C).

Figure 3. Analysis of tumorigenesis-associated gene expression levels regulated by DDX41. (A) Hierarchical clustering of the expression of tumorigenesis-associated differentially expressed genes regulated by DDX41 overexpression in HeLa cells compared with controls. (B) Gene expression profiling of tumorigenesis-associated genes in the CESC dataset. *P<0.05. DDX41, DEAD-box helicase 41; CESC, cervical and endocervical squamous cancer.

In order to identify the cellular function of DDX41-regulated alternatively spliced genes, high-confidence regulated alternatively spliced genes were analyzed by performing GO and KEGG enrichment analysis. DDX41-regulated genes with alternative splicing were enriched in a number of tumorigenesis-associated pathways, including ‘EGFR signaling’, ‘FGFR signaling’, ‘MAPK signaling’ and ‘prostate cancer’ (Fig. 4A). DDX41-regulated tumorigenesis-associated genes included FGFR1 and FGFR4 (Fig. 4B).

Figure 4. Functional clustering of alternatively spliced genes regulated by DDX41. (A) Top ten GO biological processes and KEGG functional pathways of DDX41-regulated alternatively spliced genes. (B) Results for FGFR1 ES and FGFR4 A5SS. Integrative Genomics Viewer-sashimi plots showing alternative splicing changes in Con or DDX41-OE HeLa cells. Schematic diagram depicts the structure of model splicing variant AS1 (purple) and alternative splicing variant AS2 (green). Boxes indicate exon sequences; horizontal line indicates intron sequences. RNA-seq quantification of alternative splicing regulation. *P<0.05. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes ES, exon skipping; A5SS, alternative 5′ splice site; RNA-seq, RNA-sequencing; Con, Control; OE-DDX41, overexpression DEAD-box helicase 41.

DDX41 selectively regulates transcription of immunity-associated genes in HeLa cells

In order to assess the functional role of DDX41, all 959 DEGs were subjected to GO and KEGG annotation to identify their potential biological roles. Based on the cut-off criteria, the results showed that upregulated DEGs were enriched in six GO terms and downregulated DEGs were enriched in 30 GO terms (Tables SVIII and SIX). In the GO analysis, the upregulated genes were enriched in ‘regulation of transcription’ (Fig. S1A), while the downregulated genes were enriched in ‘cell adhesion’, ‘extracellular matrix organization’ and ‘angiogenesis’ (Fig. 5A). KEGG analysis showed that upregulated genes were associated with ‘antigen processing and presentation’, ‘endocytosis’ and ‘MAPK signaling pathway’ (Fig. 5B; Table SX). Moreover, the downregulated genes were enriched in ‘cell adhesion molecules’, ‘cytokine-cytokine receptor interaction’ and ‘B cell receptor signaling pathway’ (Fig. S1B; Table SXI. DDX41-downregulated immune-associated genes included CD22, CD177, CSF2, C4A, SELPLG, TNFRSF1B, HLA-DOB, HLA-DMA, GH1 and PDGFC.

Figure 5. Analysis of immune-associated gene expression and functional clustering of differentially expressed genes regulated by DDX41. (A) Top ten GO biological processes of DDX41-downregulated genes. (B) Top ten KEGG functional pathways of DDX41-upregulated genes. (C) Transcription of genes was assessed by RNA-seq and qPCR. (D) HPV gene expression levels measured by qPCR. *P<0.05, **P<0.01 and ***P<0.001. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; RNA-seq, RNA-sequencing; qPCR, quantitative PCR; Con, Control; OE-DDX41, overexpression DEAD-box helicase 41.

In order to investigate the DDX41-regulated transcription of immune-associated genes, heatmap analysis of these genes was performed, demonstrating a high consistency for DDX41-regulated transcription in both replicates (Fig. S1C). Expression levels of these genes were confirmed via qPCR analysis, showing a significant change in expression levels of CD22, CD177, SCF2, IL12A and TNFRSF1B, in agreement with results of RNA-seq analysis (Fig. 5C). HeLa is a well-characterized HPV+ cervical cancer cell line, so activated expression of immune response genes by DDX41 overexpression may have resulted from an effort to reactivate expression of the HPV genome (38,39). In order to assess this possibility, expression levels of HPV genes E6, E7 and L1 were evaluated via qPCR with the primer pairs shown in Table I, in ref. 40. None of these HPV genes showed significantly increased expression upon DDX41 overexpression (Fig. 5D). From these results, it was concluded that DDX41 regulated immune-associated gene expression.

Antigen processing genes upregulated by DDX41 in HeLa cells are associated with survival rate in patients with cervical cancer

Previous studies have reported that antigen processing presents antigen-derived peptides into MHC-I and MHC-II molecules to be recognized by T cell receptors expressed by CD8+ and CD4+ T cells, respectively (41–43). The top KEGG pathway in analysis of DDX41 overexpression was ‘antigen processing and presentation’, with significant upregulation of HSPA1A, HSPA1B, HSPA6, HLA-DMB and HLA-G detected (Fig. 6A). Expression levels of these genes were confirmed via RT-qPCR, and the results were consistent with the RNA-seq analysis (Fig. 6A).

Figure 6. Transcription of antigen processing and presentation genes regulated by DDX41. (A) Transcription of several genes were measured by RNA-seq and qPCR. (B) Survival analysis of two genes in the Gene Expression Profiling Interactive Analysis databases of CESC (cutoff point based on the quartile method). (C) Association between DDX41 expression and immunity in CESC, reproduced from TIMER2.0 database. ***P<0.001. Con, control; OE-DDX41, overexpression DEAD-box helicase 41; RNA-seq, RNA-sequencing; qPCR, quantitative PCR; CESC, cervical and endocervical squamous carcinoma.

DDX41 regulation of the expression of immune-associated genes in HeLa cells may reflect an important role of DDX41 in cervical cancer. To test the hypothesis, the association between the expression levels of DDX41-regulated genes and the survival rates of patients with CESC was analyzed in the GEPIA databases. Expression levels of HLA-DMB were positively associated with survival rate of patients with CESC in the whole disease course of >200 months (Fig. 6B). The expression levels of HLA-G, HLA-1A and HLA-1B were positively (though not significantly) associated with survival rate in either the early or late disease stage, while those of HSPA6 were negatively associated with survival rate (Figs. 6B and S2A).

Infiltration of tumor immune cell populations in solid tumors can be estimated from the transcriptomic changes in the bulk tumor samples (44). TIMER2.0 (timer.cistrome.org/) is a user-friendly web interface for dynamic exploring and visualizing the interaction between tumor-infiltrating immune cells and important tumoral transcriptome changes using deconvolution-based approaches (45). The TIMER2.0 database was next searched to estimate the correlation between DDX41 mRNA expression and immune cell infiltration in CESC. As illustrated in the scatter plot (Fig. 6C), expression of DDX41 was positively correlated with the population of B, CD4+ Th1 and NK cells, macrophage M2 and myeloid DCs. The correlation between DDX41 expression and CD8+ T cells was not significant. Collectively, the association between DDX41 expression levels and tumor immune cell populations supported the hypothesis that DDX41-regulated expression of immune-associated and antigen presentation genes evident in HeLa may regulate immune cell infiltration and cancer immunity in patients with CESC.

DDX41 regulates the alternative splicing of certain immune-response genes

Genes regulated by DDX41-mediated alternative splicing were enriched in the ‘innate immune response’ GO pathway (MAP4K2, PGLYRP1, NFATC3, TNRC6B, RAF1, MAPKAP1, TANK, NLRX1, PLA2G6, NOD1, ZC3HAV1, PML, ATG9A, ERBB3, NR4A1, FGFR4, MLST8, FGFR1, LRRFIP1, CFD, UBA52, GRB2, PRKACB, CHID1, AGO3 and ATF2; Table SXII). Enriched KEGG pathways included those involved in ‘B cell receptor signaling pathway’ (GRB2, IFITM1, PPP3CC, RAF1 and NFATC3; Table SXIII). In order to validate the ASEs identified in the RNA-seq data, 10 potential ASEs were analyzed by RT-qPCR. Out of 11 tested events, four ASEs validated by RT-qPCR showed a significant change, which was in agreement with the RNA-seq results. The four validated splicing events were located in the genes IFITM1, CFD, NLRX1 and NFATC3 (Figs. 7 and S3).

Figure 7. Validation of DDX41-regulated ASEs in the innate immune response. Integrative Genomics Viewer-sashimi plots showing alternative splicing changes that occurred in the control or DDX41 overexpression groups of HeLa cells. The results for IFITM1 (A), TANK (B) and ZC3HAV1 (C) are shown. Transcript isoforms and the density map of the distribution all RNA-seq reads around the gene regions involved in the ASEs from the four samples are shown. Alternative splicing isoforms of each ASE are depicted with the number of reads supporting each isoform shown. Schematic diagram depicts the structure of model splicing variant AS1 (purple) and alternative splicing variant AS2 (green). Boxes indicate exon sequences; horizontal line indicates intron sequences. RNA-seq quantification of alternative splicing regulation is shown. The changed ratio of AS events in RNA-seq was calculated using the formula: AS1 junction reads/AS1 junction reads + AS2 junction reads. The altered ratio of AS events in RT-PCR was calculated using the formula: AS1 transcripts level/AS2 transcripts level. Transcripts for genes are shown. *P<0.05 and **P<0.01. Con, control; OE-DDX41, overexpression DEAD-box helicase 41; RNA-seq, RNA-sequencing; qPCR, quantitative PCR; IntronR, intron retention; A5SS, alternative 5′ splice site.

Discussion

The present study demonstrated that DDX41 acted as a tumor suppressor by inhibiting cell proliferation and promoting apoptosis in HeLa cells. To the best of our knowledge, the present study is the first to profile changes in the HeLa cell transcriptome following DDX41 overexpression. A total of 959 DEGs were identified; GO/KEGG analysis of these DEGs revealed significant enrichment in ‘cell adhesion’, ‘angiogenesis’, ‘MAPK signaling pathway’ and ‘immune response’. Genes regulated by DDX41-mediated alternative splicing were highly enriched in ‘EGFR signaling pathway’ and ‘FGFR signaling pathway’. Identification and characterization of DDX41-regulated genes in cancer cells provided mechanistic insights into the link between tumorigenesis and cancer immunity, which may aid the development of novel therapeutic strategies against cervical and endocervical cancer (46,47).

A number of previously published works have indicated the diverse tumor cell type-specific roles of DDX41 in regulating tumorigenesis. For example, knockdown (KD) of DDX41 by short hairpin RNA promoted the proliferation of K562 cells, whereas DDX41 overexpression in U937 cells inhibited growth (16,48). Compared with that of WT DDX41 or DDX41-knockdown 293 cells, overexpression of DDX41R525H mutant increased cell soft agar colony formation (48). Similarly, DDX41 knockdown in CD34+ hematopoietic progenitor cells significantly enhanced colony formation compared with controls (16). An in vivo xenograft experiment of K562 cells in which DDX41 was knocked down showed accelerated tumor growth (16). Taken together, these studies support DDX41 as a tumor suppressor gene in multiple cell lines. The present results revealed that DDX41 inhibited proliferation and promoted apoptosis of HeLa cells and expression of DDX41 was associated with immune infiltration in CESC; this extends knowledge of the tumor suppressor function of DDX41 in cancer.

In the present study, DDX41 downregulated the expression levels of oncogenes LCN2, TP63, KIT, MYB and CDK5R2. Moreover, certain tumor suppressor genes were upregulated by DDX41, such as PANO1, GLI1 and HOMER2; expression of these genes was also deregulated and associated with survival rate in patients with CESC. p63 protein is a homologue of transcription factor p53, a tumor suppressor which serves key functions in tumorigenesis and oncogenesis (49,50). p63 protein is a useful immunohistochemical marker of differentiation of squamous neoplasms within the cervix and is a marker of premalignant lesions of the cervix used to predict malignant potential (51–54). The oncogenic roles of transcription factor MYB have been widely studied in different types of cancer including leukemia, colon and breast cancer and adenoid cystic carcinoma (55–57). Overexpression of MYB has been reported to be associated with poor prognosis in leukemogenesis (58), colorectal (59) and breast cancer (60) and Burkitt lymphoma (61). Therefore, DDX41 may function in tumorigenesis by regulating the expression levels of certain oncogenes and tumor suppressors.

In the present study, DDX41 regulated alternative splicing pre-mRNAs from genes enriched in cancer progression pathways, including ‘EGFR signaling’, ‘FGFR signaling’, ‘MAPK signaling’ and ‘prostate cancer’ pathways. Tumorigenesis-associated genes were FGFR1, FGFR4, NR4A1, CSK, MLST8, IGF1R and RRAF1. Alternative inclusion of one of two unique exons results in three versions of Ig-like domain III in FGFR 1–3 (62–64), but FGFR4 is not alternatively spliced in this region (65). Two FGFR1 IIIc splice variants (FGFR1a and FGFR1β) have been shown to be expressed at similar levels in normal urothelial cells, but FGFR1β is more highly expressed in tumor cells and results in increased proliferation (66). The FGFR IIIb isoform is primarily expressed in epithelial tissue, whereas the IIIc isoform is primarily expressed in mesenchymal tissue (67). Therefore, the FGF/FGFR axis may serve an important role in cancer progression (68,69). The present finding of DDX41-regulated alternative splicing of FGFR1 and FGFR4 suggested that DDX41 may regulate tumor progression via the FGF/FGFR axis.

GO/KEGG analysis revealed that DDX41 regulated the expression of genes enriched in immune-associated pathways including ‘leukocyte migration’, ‘cytokine-cytokine receptor interaction’, ‘antigen processing and presentation’ and ‘autoimmune disease’. This finding suggested a potential function of DDX41 in regulating cancer immunity. Among these pathways, the mRNA expression levels of CD22, CD177 and CSF2 were significantly downregulated and IL-12A was upregulated by DDX41. It has been reported that administration of IL-12 promotes the development of T helper type 1 CD4+ T cells, inducing rapid onset of insulin-dependent diabetes mellitus in non-obese diabetic mice (70). The interaction of CD22 with α2,6-linked sialic acid ligands has been proposed to regulate B lymphocyte function and migration (71). DDX41 selectively upregulated antigen processing and presentation genes including HSPA1A, HSPA1B and HSPA6. When these HSPs are released into the extracellular space, they act as a source of antigens due to their ability to chaperone peptides and induce DCs to cross-present antigens to T cells, and also stimulate the innate immune system independently of peptides (72). Castelli et al (73) reported that HSP70 purified from human melanoma activated T cells recognizing melanoma differentiation antigens by gaining access to the class I HLA presentation pathway. Noessner et al (74) reported that tumor-derived HSP70-peptide complexes have the immunogenic potential to instruct cross-presentation of antigenic peptide for T cell recognition by DCs. It has been suggested that combinatorial approaches encompassing DCs and tumor-associated antigens, such as HSP70, may be useful options for cancer therapies (75,76). The unique immune properties of HSP70 have enabled development of innovative prophylactic and therapeutic vaccines, particularly those for cervical and endocervical cancer (77).

DDX41 is an RNA helicase, which also acts as a dsDNA sensor in the cytoplasm, thus defending against pathogen invasion and activating the type 1 IFN response (22). Bruton's tyrosine kinase has been shown to phosphorylate Tyr414 of DDX41, which is critical for DDX41 recognition of DNA and binding to STING (78). In a zebrafish model, DDX41 contributed to STING-STAT6-mediated chemokine production via its DEADc domain (79). STING signaling is crucial for activation of transcription of type 1 IFN production genes, and therefore protects cells against pathogen invasion and the development of cancer by promoting antitumor immune response (80–82). The present results revealed that DDX41 regulated the expression of alternative splicing of hundreds of genes. Cytoplasmic and nuclear localization of DDX41 was also demonstrated. DDX41 may mediate gene regulation indirectly by affecting the function of its protein partners and/or directly by its RNA and DNA binding activity.

In summary, the present findings suggested that DDX41 selectively regulated the alternative splicing of numerous important genes involved in tumorigenesis, as well as expression of genes involved in the immune response in HeLa cells. The expression levels of DDX41-regulated cancer-associated genes were deregulated in cervical cancer samples and the expression of DDX41-regulated antigen processing and presentation genes was associated with survival rate of patients with CESC. These findings identified a DDX41-regulatory network connecting transcription and alternative splicing that predicted the biological functions of DDX41 in suppressing tumor cell growth, as well as regulating cross-presentation of antigens and other antitumor immune responses. Further study of DDX41 may lead to novel vaccine strategies in the treatment of cancer.

Supplementary Material

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Acknowledgements

The authors would like to thank Professor. Yi Zhang (ABLife, Inc., Wuhan, China) for helpful discussions.

Funding

The present study was supported by the Health Commission of Hubei Province Scientific Research Project (grant no. WJ2019M118) and ABLife (grant no. ABL-7702157).

Availability of data and materials

The datasets generated and/or analyzed during the current study are available in the National Center for Biotechnology Information Gene Expression Omnibus repository (series accession no. GSE122986).

Authors' contributions

KQ, YZ and XLY designed and supervised the experiments. KQ, XLY, DNJ, YC, PZ, JZ, and HHX performed the experiments. YQX, YXW and YZ analyzed the data. KQ, YQX and YZ drafted the manuscript. All authors reviewed, read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References