Skip to main content

Identification and functional analysis of novel oncogene DDX60L in pancreatic ductal adenocarcinoma

Abstract

Background

Pancreatic ductal adenocarcinoma (PDAC) is a lethal cancer. Approximately 80% of patients initially diagnosed with locally advanced or metastatic disease survive only 4–11 months after diagnosis. Tremendous efforts have been made toward understanding the biology of PDAC.

Results

In this study, we first utilized next-generation sequencing technique and existing microarray datasets to identify significant differentially expressed genes between PDAC and non-tumor adjacent tissue. By comparing top significant survival genes in PDAC Gene Expression Profiling Interactive Analysis database and PDAC transcriptome data from patients, our integrated analysis discovered five potential central genes (i.e., MYEOV, KCNN4, FAM83A, S100A16, and DDX60L). Subsequently, we analyzed the cellular functions of the potential novel oncogenes MYEOV and DDX60L, which are highly expressed in PDAC cells. Notably, the knockdown of MYEOV and DDX60L significantly inhibited the metastasis of cancer cells and induced apoptosis. Further RNA sequencing analyses showed that massive signaling pathways, particularly the TNF signaling pathway and nuclear factor-kappa B (NF-κB) signaling pathway, were affected in siRNA-treated cancer cells. The siDDX60L and siMYEOV significantly inhibited the expression of chemokine CXCL2, which may potentially affect the tumor microenvironment in PDAC tissues.

Conclusions

The present findings identified the novel oncogene DDX60L, which was highly expressed in PDAC. Transcriptome profiling through siRNA knockdown of DDX60L uncovered its functional roles in the PDAC in humans.

Peer Review reports

Introduction

Pancreatic ductal adenocarcinoma (PDAC) is an extremely aggressive and deadly cancer. Approximately 80% of patients with PDAC are diagnosed with locally advanced or metastatic disease, and their prognosis remains dismal [1]. Genomic analysis of PDAC tissues showed that mutated cancer-related genes were significantly enriched, such as KRAS, TP53 [2]. These mutated genes aggregate into multiple signaling pathways, including KRAS, TGF-β, WNT, NOTCH [3]. Mutations in these genes will significantly change these signaling pathways, which play vital roles in regulating DNA repair, cell proliferation, cell survival and death, and further promote the progress of PDAC, including migration and metastasis [4, 5].

Currently, next-generation sequencing technologies are widely used for the identification of genetic alterations in cancer cells. Biological functions and molecular mechanisms of oncogenic genes and their interactions with cell signaling pathways involving cellular behaviors could also be revealed by advanced sequencing data analyses [6]. Thus far, a list of gene mutations and biomarkers, including serologic patterns, aberrant overexpressed mRNAs, miRNAs, proteins, and epigenetic signatures, have been associated with PDAC states. These could potentially be used as future early diagnostic and therapeutic strategies. In recent years, such strategies utilizing high-throughput next-generation sequencing discovered a series of novel somatic mutations and differentially expressed genes (DEGs) in solid tumors and circulating tumor cells [7,8,9]. However, the information obtained from genomic and transcriptomic sequencing data have not yet improved the care of patients with PDAC. The approach through which high-throughput genomic sequencing techniques can be applied to targeted-gene therapies for PDAC remains undetermined. Thus, there is an urgent need to identify novel oncogenes that can alter disease progression and provide guidance on therapeutic options for PDAC.

In the present study, we analyzed gene expression profiles using our RNA sequencing (RNA-seq) data and public microarray data of specimens obtained from patients with PDAC to identify highly interconnected genes that may serve as potential oncogenes for targeted therapy against PDAC. Furthermore, we investigated the function of these possible oncogenes in cell lines and revealed their associated signaling pathways that may be potentially used for targeted gene therapies.

Materials and methods

Cell culture

Human pancreatic cancer cell lines (MiaPaCa2 and PANC-1; American Type Culture Collection, Manassas, VA, USA) were maintained in Dulbecco’s modified Eagle’s medium supplemented with 2 mM glutamine, 1 mM Na-pyruvate, 100 units/ml penicillin, 100 μg/ml streptomycin, and 10% fetal bovine serum (all from Gibco, Thermo Fisher Scientific, Inc., Waltham, MA, USA) at 37 °C in a humidified atmosphere containing 10% CO2.

Tissue collection, cell collection, and RNA extraction

Pancreatic cancer tissues were collected at The First Affiliated Hospital of Kunming Medical University (Kunming, China). Following resection, the tumor and adjacent tissues were examined by a pathologist, placed in cryotubes with RNAlater reagent, and frozen in liquid nitrogen. Total RNA was isolated using RNeasy Mini kit (Qiagen, Germantown, MD, USA) according to the instructions provided by the manufacturer. PANC-1 and MiaPaCa2 cells were transfected with siRNAs for 48 h (three replicates per sample) and collected. Total RNAs were extracted using the RNeasy Mini kit (Qiagen, Germantown, MD, USA). The quantity and quality of extracted RNAs were assessed using a Nanodrop spectrophotometer and Agilent Bioanalyzer 2100 (Agilent Technologies, Inc., Santa Clara, CA, USA), respectively. Sequence libraries were prepared using a TruSeq Stranded mRNA Library Prep kit for NeoPrep, according to the instructions provided by the manufacturer, and sequenced using an Illumina HiSeq 2000 platform. All the RNA-seq data in the present study are available from the Gene Expression Omnibus repository (GSE171485 and GSE171486).

RNA-seq and bioinformatics analyses

Sequencing libraries were constructed and sequenced using the Illumina HiSeq2000 platform with SE50. A total of 21.0 ± 6.9 million reads were generated for each sample. RNA-seq data were aligned to the human reference genome (GrCH37, Ensembl build 74) using Tophat version 2.0.12 [10]. Gene expression levels, represented as fragments per kilo-base per million mapped reads (FPKM), were obtained for 63,783 genes/transcripts annotated using Ensemble GrCH37 database release 74. Functional and pathway enrichments were assessed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) bioinformatics resources. Only functional/pathway enrichments meeting a false-discovery rate < 5% are presented.

Public databases used in this study

For identifying the potential oncogenes in PDAC, we downloaded Affymetrix CEL files from the Gene Expression Omnibus database with accession number GSE28735, GSE16535, and GSE15471. For gene expression in cancer tissues, we downloaded the expression data in Gene Expression Profiling Interactive Analysis (GEPIA, http://gepia.cancer-pku.cn/). For gene expression in cancer cell lines, we downloaded the expression data in Cancer Cell Line Encyclopedia (CCLE, https://depmap.org/portal/download/) Expression 21Q2 Public.

Microarray data analysis

The microarray gene expression data were processed and analyzed using R and Bioconductor. Affymetrix U133 plus 2.0 whole genome microarrays and Affymetrix Gene Chip Human Gene ST1.0 microarrays were analyzed using R with affy and oligo packages, respectively. The microarray data were subsequently normalized using the Robust Multi-Array Average RMA method. The normalized expression values from all samples were log2 transformed. Genes differentially regulated between the biologic groups were identified using limma.

Quantitative real-time polymerase chain reaction (PCR)

For quantitative real-time PCR, 1 μg of total RNA was reverse-transcribed and triplicate PCR reactions were performed on an ABI 7500 Real-Time PCR System (Foster City, CA, USA). Glyceraldehyde-3-phosphate dehydrogenase mRNA was used as internal control. The PCR reaction was conducted according to the recommended protocol. The primers used for quantitative PCR are described in Supplementary Table S1.

Cell viability assay

Cell viability was determined using the Cell Counting Kit-8 (CCK-8) assay. Briefly, fresh culture medium containing 10 μl CCK-8 (Dojindo Molecular Technologies, Inc.) was placed in each well. The culture plates were incubated for an additional 4 h at 37 °C, and absorbance was measured at 450 nm using a microplate spectrophotometer (Molecular Devices Corp., Sunnyvale, CA, USA).

Cell Transwell assay

The chambers were placed in a 24-well plate; serum-free medium (100 μl) was added to the upper chamber, and the plate was placed in an incubator at 37 °C for 1 h. A serum-free cell suspension was prepared, and 100 μl (105 cells) were placed on the plate. Next, 30% fetal bovine serum medium (600 μl) was added in the lower chamber, and the cells were incubated at 37 °C. The cells were fixed with 4% paraformaldehyde. The cells were counted using microscope photos (original magnification × 200), the data were analyzed, and differences in cell migratory ability between the experimental group (siRNAs) and the control group (siControl) were determined.

The Kaplan–Meier plotter analysis

To explore the potential prognostic value of the genes in cancer patients, the Kaplan–Meier Plotter analysis (http://kmplot.com/analysis/) pan-cancer RNA-seq database was used to perform overall survival analysis (PDAC, n = 177) [11]. In accordance with the instructions, the prognostic value of two genes were also analyzed with the Kaplan–Meier Plotter (PDAC), and GraphPad Prism 8.3.0 was used to further display the data.

Statistical analysis

All experiments in this study were independently performed in triplicate. The data are presented as the mean with standard error of the mean or standard deviation. All graphs were plotted and analyzed using the GraphPad Prism 8 software (San Diego, CA, USA). A p-value < 0.05 denoted statistically significant difference.

Results

PDAC DEGs identified from RNA-seq and microarray analyses

For the identification of specific genes expressed in PDAC, we constructed sequencing libraries from polyadenylated-RNA extracted from six PDAC specimens, three non-tumor adjacent tissues, and three pancreatic tissues from normal individuals. Approximately 21 M raw reads were obtained. We aligned the RNA-seq data, yielding an average mapping rate of 97.0 ± 0.6%. The quality of the RNA-seq data was determined using the Pearson correlation of the associations between the transcriptome data obtained from different specimens. RNA-seq data indicated that the Pearson correlation coefficients of transcriptome data with same phenotypes (PDAC versus PDAC: 0.87 and Control versus Control: 0.92) were higher than those noted with different phenotypes (PDAC versus Control: 0.66) (Figs. 1A–C). These findings suggested that non-tumor pancreas and PDAC specimens show distinct global gene expression patterns. The heatmap of a total of 1371 unique DEGs identified from the comparison between six PDAC and six non-tumor control pancreas specimens distinguished PDAC versus non-tumor controls (Fig. 1D). Of those, 607 genes were upregulated and 763 genes were downregulated in PDAC specimens using cutoff criteria of |log2(fold-change [FC])| > 0.5 (log2(FC) > 0.5) and p < 0.05 (the full list of DEGs is provided in Supplementary Table S2A).

Fig. 1
figure 1

Identification of differentially expressed genes (DEGs) in PDAC using RNA-seq and microarray data. (A–C) Scatter plots of FPKM values for gene expression of PDAC versus PDAC (A), non-tumor pancreas versus non-tumor pancreas (B), and PDAC versus matching non-tumor pancreas (C). Pearson correlation coefficients are shown above the corresponding scatter plots. (D) Heatmap of 1371 DEGs identified from the comparison between six PDAC and six non-tumor control specimens (p < 0.05 and |log2FC| > 0.5). (E, F) Venn diagram of DEGs identified from previous published microarray studies of PDAC versus non-tumor controls using two types of Affymetrix microarray platforms. Overlapping upregulated DEGs (n = 795) (E) and downregulated DEGs (n = 252) (F) are shown. (G) Heatmaps and hierarchical clustering of all PDAC specimens versus non-tumor controls from three independent groups clearly show that these 1047 DEGs can distinguish PDAC from non-tumor pancreas specimens FC, fold change; FPKM, fragments per kilo-base per million mapped reads; PDAC, pancreatic ductal adenocarcinoma; RNA-seq, RNA sequencing.

AffymetrixHG-U133 Plus 2.0 and Human Gene 1.0 ST microarray platforms have been commonly used in previous studies [12,13,14] to study the gene expression profiles of the whole genome. This approach resulted in the discovery novel tumor biomarker genes and investigated the molecular mechanisms for PDAC [15]. We sought to better understand the biology of cancer and compare the numbers and percentages of overlapping DEGs between RNA-seq and other independent data sets on microarray platforms. For this purpose, we revisited the gene expression profiles of PDAC versus non-tumor pancreatic tissues by reprocessing the gene expression microarray data from three independent studies [12,13,14]. We combined two data sets on the same Affymetrix HG-U133 Plus 2.0 platform and processed the combined data sets (GSE15471 and GSE16515) [12, 13] using R and Bioconductor packages “affy” and “affycoretools”. In the comparison between 75 PDAC and 55 non-tumor pancreatic tissue gene expression profiles, 7776 probes which reflect a total of 5598 unique genes (log2(FC) > 0.5 and p < 0.001) were revealed to be differentially expressed between PDAC and non-tumor pancreas tissues. The identified DEGs consist of 3636 unique upregulated genes (5388 probes) and 1962 unique downregulated genes (2388 probes) in PDAC (the full list of DEGs is provided in Supplementary Table S2B). Meanwhile, using R and the Bioconductor package “oligo”, we obtained a list of unique DEGs by reanalyzing the dataset GSE28735 of 90 human PDAC and adjacent non-tumor tissues (45 PDAC versus 45 non-tumor specimens) reported by Zhang [14] using the Affymetrix Human Gene 1.0 ST microarray platform. A total of 1546 probes (1473 unique DEGs) were identified; 1006 probes (961 unique genes) were upregulated and 540 probes (512 unique genes) were downregulated in PDAC with cutoff criteria of |log2(FC)| > 0.5 and p < 0.001 (the full list of DEGs is provided in Supplementary Table S2C).

To determine which overrepresented signaling pathways are shared among all PDAC data sets, we compared the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways identified from all lists of DEGs. DAVID-KEGG analysis of these DEGs revealed that the shared signaling pathways for upregulated genes across all data sets at three platforms were the signaling pathways of the pathways in cancer, p53, extracellular matrix-receptor, apoptosis, and cell cycle (Supplementary Table S3). However, there were no shared KEGG pathways identified for downregulated genes in the DEG lists in these data sets.

Identification of potential oncogenes responsible for the progression of PDAC

We subsequently compared the lists of DEGs identified from the data sets on two independent Affymetrix microarray platforms. Approximately 82% (795/961) of the upregulated DEGs and 49% (252/512) of the downregulated DEGs from Human Gene 1.0 ST microarray data were also identified in the HG-U133 Plus 2.0 microarray platform (Fig. 1E and F). Heatmaps and hierarchical clustering of all PDAC versus non-tumor pancreatic specimens obtained from three independent groups clearly showed that these 1047 DEGs can distinguish between PDAC and non-tumor pancreas specimens (Fig. 1G). With higher statistical stringency, the comparison between two lists of DEGs from microarray data revealed 204 upregulated and 39 downregulated overlapped DEGs in PDAC using cutoff criteria of |log2(FC)| > 1 and p < 0.0001 (the full lists of overlapped DEGs are provided in Supplementary Table S4A).

Next, we compared the list of DEGs from our PDAC RNA-seq data with the lists of DEGs identified from the above microarray data sets. The cross-platform comparison revealed that approximately 43% (260/607) and 64% (387/607) of the upregulated DEGs identified through RNA-seq overlapped with the lists of DEGs obtained from the HG 1.0 ST and HG-U133 Plus 2.0 platforms, respectively (Fig. 2A). Conversely, only 6% (47/764) and 15% (116/764) of the downregulated DEGs recorded from the RNA-seq data were identified from the HG 1.0 ST and HG-U133 Plus 2.0 data sets, respectively (Fig. 2B). The cross-platform analysis identified 227 upregulated and 32 downregulated genes in PDAC that were shared in both RNA-seq data and two other microarray data sets (the full list of DEGs is provided in Supplementary Table S4B).

Fig. 2
figure 2

Identification of overlapped genes responsible for the progression of PDAC. (A-B) Venn diagram of DEGs identified in the RNA-seq analysis of this study and in previous microarray studies of PDAC versus non-tumor controls. Overlapping upregulated (A) and downregulated (B) DEGs across multiple platforms are shown. (C) Venn diagram of the identified DEGs and the Top 100 overall survival-related genes in PDAC. Five most significant potential oncogenes were identified: MYEOV, KCNN4, FAM83A, S100A16, and DDX60L. (D) Expression of five significant potential oncogenes in PDAC tissues of the GEPIA database. One-way ANOVA is used to analyze significant differences in the gene expression levels between PDAC normal tissues and cancer tissues. *represents significant differences in RNA expression, p < 0.01. DDX60L, DExD/H-box 60 like; DEGs, differentially-expressed genes; FAM83A, family with sequence similarity 83 member A; GEPIA, Gene Expression Profiling Interactive Analysis; KCNN4, potassium calcium-activated channel subfamily N member 4; MYEOV, myeloma overexpressed; PDAC, pancreatic ductal adenocarcinoma; RNA-seq, RNA sequencing; S100A16, S100 calcium binding protein A16

For the identification of the most significant potential oncogenes in PDAC, we overlapped these 259 identified genes with the top 100 significant survival genes in PDAC (Gene Expression Profiling Interactive Analysis [GEPIA] database); five genes (myeloma overexpressed [MYEOV], potassium calcium-activated channel subfamily N member 4 [KCNN4], family with sequence similarity 83 member A [FAM83A], S100 calcium binding protein A16 [S100A16], and DExD/H-box 60 like [DDX60L]) were identified (Fig. 2C, Supplementary Table S4C). Kaplan–Meier Plotter analysis showed that the expression levels of these five genes are also closely related to the prognosis of patients with PDAC. Of note, patients with pancreatic cancer expressing high levels of these potential oncogenes have a very poor prognosis (p < 0.001) (Supplementary Fig. S1A). In the GEPIA database, these potential oncogenes were highly expressed in the cancer tissues of patients with pancreatic cancer (Fig. 2D), and their expression levels were positively correlated with the cancer stages (Supplementary Fig. S1B). Compared with the early-stage (stages I–II), pancreatic cancer tissues obtained from patients with stage III–IV disease showed markedly higher expression levels of these potential oncogenes. The results indicated that these genes may be involved in the progression of pancreatic cancer (Supplementary Fig. S1B). In summary, based on our RNA-seq and public microarray data, we identified five potential oncogenes that may play vital roles in the occurrence and progression of PDAC.

Expression and siRNA knockdown of novel oncogene DDX60L in PDAC

We investigated the expression of these potential oncogenes (MYEOV, KCNN4, FAM83A, S100A16, and DDX60L) in the pancreatic cancer cell lines of the Cancer Cell Line Encyclopedia (CCLE) database. The results showed that MYEOV, KCNN4, S100A16, and DDX60L are highly expressed in most pancreatic cancer cell lines, but FAM83A is not (Fig. 3A). Further quantitative real-time-PCR validation performed using the MiaPaCa2 cell line was consistent with the CCLE data (Fig. 3B). Recent studies have shown that MYEOV [16,17,18], KCNN4 [19, 20], and S100A16 [21,22,23] act as oncogenes and play vital roles in cancer progression. In this study, we found a novel potential oncogene – DDX60L. To validate the potential functions of DDX60L in pancreatic cancer (MYEOV was used as positive control), we designed three siRNAs to knockdown its expression in cancer cells (sequences are shown in Supplementary Table S1B). The knockdown efficiency of these siRNAs was tested using the MiaPaCa2 cell line. The analysis showed that siMYEOV-3 and siDDX60L-2 resulted in effective inhibition of target gene expression (Fig. 3C). Next, we used these validated siRNAs to inhibit the expression of target genes in MiaPaCa2 and PANC-1 cells, demonstrating excellent effectiveness (> 60%) (Fig. 3D and E).

Fig. 3
figure 3

Knockdown of MYEOV, KCNN4, S100A16, and DDX60L expression through siRNA. (A) Expression of MYEOV, KCNN4, FAM83A, S100A16, and DDX60L in pancreatic cancer cells using data from the CCLE database (CCLE-PAAD). (B) Expression of MYEOV, KCNN4, FAM83A, S100A16 and DDX60L in MiaPaca2 pancreatic cancer cell (compared with the expression of GAPDH). (C) Screening for siRNAs that could knock down the expression of MYEOV and DDX60L in MiaPaca2 cells. (D, E) Knockdown of MYEOV and DDX60L expression with indicated siRNAs in MiaPaca2 cells (D) and PANC-1 cells (E)

One-way ANOVA is used to analyze significant differences in the gene expression levels between control cancer cells and siRNA treated cancer cells. *p < 0.05; **p < 0.01; ***p < 0.001; ***p < 0.0001. CCLE, Cancer Cell Line Encyclopedia; DDX60L, DExD/H-box 60 like; FAM83A, family with sequence similarity 83 member A; GAPDH, glyceraldehyde-3-phosphate dehydrogenase; KCNN4, potassium calcium-activated channel subfamily N member 4; MYEOV, myeloma overexpressed; PAAD, pancreatic adenocarcinoma; S100A16, S100 calcium binding protein A16.

Knockdown of DDX60L inhibited the migration of PDAC cells, but not their proliferation

DDX60L is highly expressed in pancreatic cancer and related to the overall survival of patients with cancer. Therefore, it may play a vital role in the progression of pancreatic cancer. Using siRNA knockdown of DDX60L, we further investigated its potential functions in MiaPaCa2 and PANC-1 cells. The CCK-8 assay of siRNA knockdown target genes showed that siMYEOV significantly reduced the proliferation of MiaPaCa2 and PANC-1 cells (p < 0.001), whereas siDDX60L had no effect on proliferation (Fig. 4A and B). The effect of DDX60L on the migration of pancreatic cancer cells was investigated using Transwell migration assay. The results showed that siRNA knockdown of MYEOV and DDX60L could significantly inhibit the migration of MiaPaCa2 and PANC-1 cells (Fig. 4C and D).

Fig. 4
figure 4

Knockdown of DDX60L inhibited the migration of PDAC cells, but not their proliferation. (A, B) Viability of siRNA-transfected PDAC cells. CCK-8 assay analyzed the siControl-, siMYEOV-, and siDDX60L-transfected PANC-1 (A) and MiaPaca2 (B) cells. (C, D) The invasive ability of siRNA-transfected PDAC cells. Transwell assay analyzed the siControl-, siMYEOV-, and siDDX60L-transfected PANC-1 (C) and MiaPaca2 (D) cells. One-way ANOVA is used to analyze significant differences. *p < 0.05; **p < 0.01; ***p < 0.001. CCK-8, Cell Counting Kit-8; DDX60L, DExD/H-box 60 like; MYEOV, myeloma overexpressed; PDAC, pancreatic ductal adenocarcinoma

Knockdown of DDX60L induced apoptosis of PDAC cells, but not cell cycle arrest

Since DDX60L did not affect the proliferation of pancreatic cancer cells, we also investigated the rate of apoptosis and cell cycle in pancreatic cancer cells. The analysis was performed in siRNA-treated or -untreated pancreatic cancer cells through flow cytometry. The apoptosis analysis of siRNA-treated MiaPaCa2 and PANC-1 cells showed that, compared with control siRNA, siDDX60L (p < 0.001 in MiaPaCa2; p < 0.01 in PANC-1 cells) and siMYEOV (p < 0.05 in MiaPaCa2; p < 0.05 in PANC-1 cells) significantly induced cell apoptosis (Fig. 5A and B). The cell cycle analysis of siRNA-treated MiaPaCa2 and PANC-1 cells showed that, compared with control siRNA, siDDX60L and siMYEOV did not have any significant effects on the cell cycle (Fig. 5C and D). In summary, these results indicated that DDX60L could inhibit the migration of pancreatic cancer cells and induce their apoptosis, which may play a vital role in the progression of pancreatic cancer.

Fig. 5
figure 5

Knockdown of DDX60L induced apoptosis of PDAC cells, but not cell cycle arrest. (A, B) Apoptosis of siRNA-transfected PDAC cells. FACS analyzed the apoptosis of the siControl-, siMYEOV-, and siDDX60L-transfected PANC-1 (A) and MiaPaca2 (B) cells. (C, D) Cell cycle of siRNA-transfected PDAC cells. FACS analyzed the cell cycle of the siControl-, siMYEOV-, and siDDX60L-transfected PANC-1 (C) and MiaPaca2 (D) cells. One-way ANOVA is used to analyze significant differences. Ns, no significant differences, *p < 0.05; **p < 0.01; ***p < 0.001. DDX60L, DExD/H-box 60 like; FACS, fluorescence-activated cell sorting; MYEOV, myeloma overexpressed; PDAC, pancreatic ductal adenocarcinoma

Transcriptome analyses revealed alterations of genes and signaling pathways in PDAC after knockdown of DDX60L

The transcriptome alterations induced by siDDX60L were also investigated through RNA-seq. For the transcriptome analysis, we acquired approximately 26 M clean reads from siRNA-transfected and control MiaPaCa2 and PANC-1 cells with Q30 > 90%, respectively. After bioinformatics analysis of these RNA-seq data, the upregulated and downregulated DEGs were identified.

For the positive control gene, MiaPaCa2 and PANC-1 cells were transfected with siMYEOV; 618 downregulated genes and 513 upregulated genes were identified (Fig. 6A). The KEGG_PATHWAY category of these DEGs showed that the upregulated genes were clustered into seven subcategories (p < 0.05; three subcategories with p < 0.01), including alcoholism (p = 0.0014), systemic lupus erythematosus (p = 0.0022), and viral carcinogenesis (p = 0.0043) (Fig. 6B, Supplementary Table S5A). The downregulated genes were clustered into 27 subcategories (p < 0.05; 10 subcategories with p < 0.01), including the tumor necrosis factor (TNF) signaling pathway (p = 3.92*E-05), small cell lung cancer (p = 0.0012), p53 signaling pathway (p = 0.003), nuclear factor-kappa B (NF-κB) signaling pathway (p = 0.005) and cell cycle (p = 0.017) (Fig. 6B, Supplementary Table S5A). In particular, siMYEOV significantly inhibited the expression of chemokine C-X-C motif chemokine ligand 1 (CXCL1), chemokine C-X-C motif chemokine ligand 2 (CXCL2), and chemokine C-X-C motif chemokine ligand 3 (CXCL3) in the TNF signaling pathway. These results suggested that siRNA knockdown of MYEOV repressed the expression of genes involved in the cell cycle, and significantly influenced the immune response of PDAC cells.

Fig. 6
figure 6

Transcriptome analyses showed alterations of genes and signaling pathways in siMYEOV- and siDDX60L-transfected PDAC cells. (A) Heatmap and Venn diagram of DEGs identified in siMYEOV-transfected MiaPaca2 and PANC-1 cells. (B) DAVID-KEGG analysis of the DEGs in siMYEOV-transfected MiaPaca2 and PANC-1 cells (p < 0.05). (C) Heatmap and Venn diagram of DEGs identified in siDDX60L-transfected MiaPaca2 and PANC-1 cells. (D) DAVID-KEGG analysis of the DEGs in siDDX60L-transfected MiaPaca2 and PANC-1 cells (p < 0.05). DAVID-KEGG, Database for Annotation, Visualization and Integrated Discovery-Kyoto Encyclopedia of Genes and Genomes; DDX60L, DExD/H-box 60 like; DEGs, differentially-expressed genes; PDAC, pancreatic ductal adenocarcinoma; MYEOV, myeloma overexpressed

In siDDX60L-transfected MiaPaCa2 and PANC-1 cells, 214 downregulated genes and 186 upregulated genes were identified (Fig. 6C). These DEGs were subjected to DAVID analysis (a web-based high-throughput functional genomics analysis tool) for a systematic clustering of these genes. In the KEGG_PATHWAY category, we found that the upregulated genes were clustered into 16 subcategories (p < 0.05; one subcategory with p < 0.01), including the pathways in cancer (p = 0.011) and purine metabolism (p = 0.012) (Fig. 6D, Supplementary Table S5B). The downregulated genes were clustered into 15 subcategories (p < 0.05; eight subcategories with p < 0.01), including the TNF signaling pathway (p = 8.91*E-07), NF-κB signaling pathway (p = 0.0001), insulin signaling pathway (p = 0.0018), and nucleotide-binding oligomerization domain (NOD)-like receptor signaling pathway (p = 0.0055) (Fig. 6D, Supplementary Table S5B). Interestingly, siDDX60L significantly inhibited the expression of chemokine CXCL2 in the TNF signaling pathway. These results indicate that siRNA knockdown of DDX60L significantly altered the TNF and NF-κB signaling pathways in PDAC cells. Hence, the downregulated chemokines may influence the immune response of PDAC cells.

Knockdown of DDX60L and MYEOV inhibited the expression of CXCL2 and potentially affected the prognosis in PDAC

Interestingly, the RNA-seq data of DDX60L and MYEOV showed some similar alterations of the signaling pathways in PDAC cells, such as the TNF signaling pathway and NF-κB signaling pathway. These pathways were significantly inhibited in PDAC cells treated with siDDX60L and siMYEOV. Similarly, we found that the expression of neutrophil chemotactic chemokine CXCL2 was significantly inhibited in PDAC cells treated with siDDX60L and siMYEOV (Fig. 7A). With the Kaplan–Meier Plotter analysis, the expression of CXCL2 in pancreatic cancer significantly affects the prognosis of patients (n = 177), and high expression levels are associated with a poor prognosis (Fig. 7B). Furthermore, by combining the expression levels of DDX60L and MYEOV, the cohorts of CXCL2lowDDX60Llow and CXCL2lowMYEOVlow showed a better prognosis in patients with PDAC from TCGA database (Fig. 7C and D). In summary, these data indicated that the four genes identified in this study affect the cellular functions of PDAC. These genes may also affect the microenvironment of cancer cells through transcriptomic alterations of chemokines.

Fig. 7
figure 7

Knockdown of DDX60L and MYEOV inhibited the expression of CXCL2 and led to good prognosis in PDAC. (A) Expression of CXCL2 in siDDX60L- and siMYEOV-transfected MiaPaCa2 and PANC-1 cells. Two-way ANOVA is used to analyze significant differences in the gene expression levels between control cancer cells and siRNA treated cancer cells. ****p < 0.0001. (B) The Kaplan–Meier Plotter analysis showed that patients from the TCGA-PDAC dataset (n = 177) [11] who expressed high levels of CXCL2 had a poor prognosis. (C-D) The Kaplan–Meier Plotter analysis showed that PDAC patients with low CXCL2 and DDX60L expression in tumors (C) and those with low CXCL2 and MYEOV expression in tumors (D) had a significantly improved overall survival. The p-values are shown. CXCL2, C-X-C motif chemokine ligand 1; DDX60L, DExD/H-box 60 like; MYEOV, myeloma overexpressed; OS, overall survival; PDAC, pancreatic ductal adenocarcinoma; TCGA, The Cancer Genome Atlas

Discussion

In this study, we integrated a large number of transcriptome profiling data from PDAC and non-tumor pancreatic tissue specimens to identify significant differentially expressed genes. A list of aberrantly overexpressed genes was associated with PDAC states. Among these identified genes, we found five genes related to cancer prognosis, and further verified their cellular functions and potential signaling pathways in pancreatic cancer cells. Therefore, it is reasonable to hypothesize that these genes could be served as immunohistochemical prognostic markers and potential therapeutic targets for PDAC.

Interestingly, the five genes identified in this study (DDX60L, FAM83A, KCNN4, MYEOV, and S100A16) showed some similar characteristics; they were highly expressed and corelated with poor prognosis in pancreatic cancer. Previously, MYEOV [16,17,18], KCNN4 [19, 20], S100A16 [21,22,23], and FAM83A [24, 25] had been reported as oncogenes in various types of cancer, which could promote the proliferation and invasion of cancer cells. However, the novel oncogene DDX60L we identified did not influence the proliferation of cancer cells, but significantly influenced their invasive ability. Further RNA-seq analysis showed that siRNA knockdown of DDX60L significantly inhibited the genes in the TNF signaling and NF-κB signaling pathway, particularly the chemotactic cytokine CXCL2.

Chemotactic cytokines, also termed chemokines, play a key role in mediating the recruitment of immune cells to tumor sites [26, 27]. CXC chemokine CXCL2 possesses potent neutrophil chemotactic activity [28, 29]. The involvement of tumor monocyte-derived chemokines and cytokines in modulating neutrophil accumulation and functions was previously studied [30, 31]. The knockdown of oncogenes MYEOV and DDX60L inhibited the expression of CXCL2 in pancreatic cancer cells. These findings suggested that repression of these oncogenes could potentially reduce the metastasis of cancer cells and repress the immunosuppressive ability of cancer cells. However, these results should be further validated in mouse xenograft models of pancreatic cancer.

Conclusions

The present study integrated a systematic approach to identify key genes associated with PDAC tumorigenesis based on gene expression profiling data. Several key genes including MYEOV, KCNN4, FAM83A, S100A16, and DDX60L were identified. Further cellular experiments validated the function of the novel oncogene DDX60L. These genes could potentially serve as targets and for tumor imaging to guide therapeutic selection in PDAC.

Availability of data and materials

All the RNA-seq data in the present study are available from the Gene Expression Omnibus repository GSE171485 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE171485) and GSE171486 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE171486).

References

  1. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2020. CA Cancer J Clin. 2020;70(1):7–30. https://doi.org/10.3322/caac.21590.

    Article  PubMed  Google Scholar 

  2. Collisson EA, Maitra A. Pancreatic Cancer genomics 2.0: profiling metastases. Cancer Cell. 2017;31(3):309–10. https://doi.org/10.1016/j.ccell.2017.02.014.

    Article  CAS  PubMed  Google Scholar 

  3. Bailey P, Chang DK, Nones K, Johns AL, Patch AM, Gingras MC, et al. Genomic analyses identify molecular subtypes of pancreatic cancer. Nature. 2016;531(7592):47–52. https://doi.org/10.1038/nature16965.

    Article  CAS  PubMed  Google Scholar 

  4. Sanchez-Vega F, Mina M, Armenia J, Chatila WK, Luna A, La KC, et al. Oncogenic signaling pathways in the Cancer genome atlas. Cell. 2018;173(2):321–37. https://doi.org/10.1016/j.cell.2018.03.035.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Ding L, Bailey MH, Porta-Pardo E, Thorsson V, Colaprico A, Bertrand D, et al. Perspective on oncogenic processes at the end of the beginning of Cancer genomics. Cell. 2018;173(2):305–20. https://doi.org/10.1016/j.cell.2018.03.033.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Thomas JK, Kim MS, Balakrishnan L, Nanjappa V, Raju R, Marimuthu A, et al. Pancreatic Cancer database: an integrative resource for pancreatic cancer. Cancer Biol Therapy. 2014;15(8):963–7. https://doi.org/10.4161/cbt.29188.

    Article  CAS  Google Scholar 

  7. Biankin AV, Waddell N, Kassahn KS, Gingras MC, Muthuswamy LB, Johns AL, et al. Pancreatic cancer genomes reveal aberrations in axon guidance pathway genes. Nature. 2012;491(7424):399–405. https://doi.org/10.1038/nature11547.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Jones S, Hruban RH, Kamiyama M, Borges M, Zhang X, Parsons DW, et al. Exomic sequencing identifies PALB2 as a pancreatic cancer susceptibility gene. Science. 2009;324(5924):217. https://doi.org/10.1126/science.1171202.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Wang L, Tsutsumi S, Kawaguchi T, Nagasaki K, Tatsuno K, Yamamoto S, et al. Whole-exome sequencing of human pancreatic cancers and characterization of genomic instability caused by MLH1 haploinsufficiency and complete deficiency. Genome Res. 2012;22(2):208–19. https://doi.org/10.1101/gr.123109.111.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and cufflinks. Nat Protoc. 2012;7(3):562–78. https://doi.org/10.1038/nprot.2012.016.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Nagy A, Munkacsy G, Gyorffy B. Pancancer survival analysis of cancer hallmark genes. Sci Rep. 2021;11(1):6047. https://doi.org/10.1038/s41598-021-84787-5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Badea L, Herlea V, Dima SO, Dumitrascu T, Popescu I. Combined gene expression analysis of whole-tissue and microdissected pancreatic ductal adenocarcinoma identifies genes specifically overexpressed in tumor epithelia. Hepato-gastroenterology. 2008;55(88):2016–27.

    CAS  PubMed  Google Scholar 

  13. Pei H, Li L, Fridley BL, Jenkins GD, Kalari KR, Lingle W, et al. FKBP51 affects cancer cell response to chemotherapy by negatively regulating Akt. Cancer Cell. 2009;16(3):259–66. https://doi.org/10.1016/j.ccr.2009.07.016.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Zhang G, Schetter A, He P, Funamizu N, Gaedcke J, Ghadimi BM, et al. DPEP1 inhibits tumor cell invasiveness, enhances chemosensitivity and predicts clinical outcome in pancreatic ductal adenocarcinoma. PLoS One. 2012;7(2):e31507. https://doi.org/10.1371/journal.pone.0031507.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Zhang G, He P, Tan H, Budhu A, Gaedcke J, Ghadimi BM, et al. Integration of metabolomics and transcriptomics revealed a fatty acid network exerting growth inhibitory effects in human pancreatic cancer. Clinical Cancer Res. 2013;19(18):4983–93. https://doi.org/10.1158/1078-0432.CCR-13-0209.

    Article  CAS  Google Scholar 

  16. Zhang R, Ma A. High expression of MYEOV reflects poor prognosis in non-small cell lung cancer. Gene. 2020;770:145337. https://doi.org/10.1016/j.gene.2020.145337.

    Article  CAS  PubMed  Google Scholar 

  17. Liang E, Lu Y, Shi Y, Zhou Q, Zhi F. MYEOV increases HES1 expression and promotes pancreatic cancer progression by enhancing SOX9 transactivity. Oncogene. 2020;39(41):6437–50. https://doi.org/10.1038/s41388-020-01443-4.

    Article  CAS  PubMed  Google Scholar 

  18. Lawlor G, Doran PP, MacMathuna P, Murray DW. MYEOV (myeloma overexpressed gene) drives colon cancer cell migration and is regulated by PGE2. J Exp Clin Cancer Res. 2010;29(1):81. https://doi.org/10.1186/1756-9966-29-81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Jiang SH, Zhu LL, Zhang M, Li RK, Yang Q, Yan JY, et al. GABRP regulates chemokine signalling, macrophage recruitment and tumour progression in pancreatic cancer through tuning KCNN4-mediated ca (2+) signalling in a GABA-independent manner. Gut. 2019;68(11):1994–2006. https://doi.org/10.1136/gutjnl-2018-317479.

    Article  CAS  PubMed  Google Scholar 

  20. Wen J, Lin B, Lin L, Chen Y, Wang O. KCNN4 is a diagnostic and prognostic biomarker that promotes papillary thyroid cancer progression. Aging (Albany NY). 2020;12(16):16437–56. https://doi.org/10.18632/aging.103710.

    Article  CAS  Google Scholar 

  21. Ou S, Liao Y, Shi J, Tang J, Ye Y, Wu F, et al. S100A16 suppresses the proliferation, migration and invasion of colorectal cancer cells in part via the JNK/p38 MAPK pathway. Mol Med Rep. 2021;23(2):164. https://doi.org/10.3892/mmr.2020.11803.

    Article  CAS  PubMed  Google Scholar 

  22. Fang D, Zhang C, Xu P, Liu Y, Mo X, Sun Q, et al. S100A16 promotes metastasis and progression of pancreatic cancer through FGF19-mediated AKT and ERK1/2 pathways. Cell Biol Toxicol. 2021;37(4):555–71. https://doi.org/10.1007/s10565-020-09574-w.

    Article  CAS  PubMed  Google Scholar 

  23. Li T, Ren T, Huang C, Li Y, Yang P, Che G, et al. S100A16 induces epithelial-mesenchymal transition in human PDAC cells and is a new therapeutic target for pancreatic cancer treatment that synergizes with gemcitabine. Biochem Pharmacol. 2020;189:114396. https://doi.org/10.1016/j.bcp.2020.114396.

    Article  CAS  PubMed  Google Scholar 

  24. Zhou F, Wang X, Liu F, Meng Q, Yu Y. FAM83A drives PD-L1 expression via ERK signaling and FAM83A/PD-L1 co-expression correlates with poor prognosis in lung adenocarcinoma. Int J Clin Oncol. 2020;25(9):1612–23. https://doi.org/10.1007/s10147-020-01696-9.

    Article  CAS  PubMed  Google Scholar 

  25. Zheng YW, Li ZH, Lei L, Liu CC, Wang Z, Fei LR, et al. FAM83A promotes lung Cancer progression by regulating the Wnt and hippo signaling pathways and indicates poor prognosis. Front Oncol. 2020;10:180. https://doi.org/10.3389/fonc.2020.00180.

    Article  PubMed  PubMed Central  Google Scholar 

  26. Nagarsheth N, Wicha MS, Zou W. Chemokines in the cancer microenvironment and their relevance in cancer immunotherapy. Nat Rev Immunol. 2017;17(9):559–72. https://doi.org/10.1038/nri.2017.49.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Yang Z, Feng J, Xiao L, Chen X, Yao Y, Li Y, et al. Tumor-derived peptidoglycan recognition protein 2 predicts survival and antitumor immune responses in hepatocellular carcinoma. Hepatology. 2020;71(5):1626–42. https://doi.org/10.1002/hep.30924.

    Article  CAS  PubMed  Google Scholar 

  28. Greer A, Irie K, Hashim A, Leroux BG, Chang AM, Curtis MA, et al. Site-specific neutrophil migration and CXCL2 expression in periodontal tissue. J Dent Res. 2016;95(8):946–52. https://doi.org/10.1177/0022034516641036.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Mohanty SK, Ivantes CA, Mourya R, Pacheco C, Bezerra JA. Macrophages are targeted by rotavirus in experimental biliary atresia and induce neutrophil chemotaxis by Mip2/Cxcl2. Pediatr Res. 2010;67(4):345–51. https://doi.org/10.1203/PDR.0b013e3181d22a73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Girbl T, Lenn T, Perez L, Rolas L, Barkaway A, Thiriot A, et al. Distinct compartmentalization of the chemokines CXCL1 and CXCL2 and the atypical receptor ACKR1 determine discrete stages of neutrophil diapedesis. Immunity. 2018;49(6):1062–76. https://doi.org/10.1016/j.immuni.2018.09.018.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Peng ZP, Jiang ZZ, Guo HF, Zhou MM, Huang YF, Ning WR, et al. Glycolytic activation of monocytes regulates the accumulation and function of neutrophils in human hepatocellular carcinoma. J Hepatol. 2020;73(4):906–17. https://doi.org/10.1016/j.jhep.2020.05.004.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

The authors gratefully acknowledge the excellent technical assistance of Dr. Limei Cao and Dr. Yongjin Zhang (NHC Key Laboratory of Drug Addiction Medicine, The First Affiliated Hospital of Kunming Medical University).

Funding

The present study was supported by grants from the China Postdoctoral Science Foundation (grant no. 2020M673596XB); National Natural Science Foundation of China (grant nos. 81903046 and 31860306); and Science and Technology Department of Yunnan Province (grant nos. ZX2019-03-04). The funders had no role in the study design, data collection, and analysis, decision to publish, or preparation of the manuscript.

Author information

Authors and Affiliations

Authors

Contributions

JY, HW and XZ conceived the study. HW and XT performed all experiments. XL, ZW, JS and ZL prepared the reagents and samples. JY and HW analyzed the data and wrote the manuscript. JY, HW and WT reviewed the results and participated in the discussion regarding the content of the manuscript. HW and WT supplemented the experiment and revised the manuscript. All authors read and approved the manuscript.

Corresponding authors

Correspondence to Juehua Yu, Zhihua Wang or Xiaosong Zhu.

Ethics declarations

Ethics approval and consent to participate

All sample collections, patient consent procedure, and patient documentation review and approval were performed using protocols approved by the institutional review board of The First Affiliated Hospital of Kunming Medical University. All patients provided written informed consent. This study was conducted in accordance with the tenets of the Declaration of Helsinki.

Consent for publication

Not applicable.

Competing interests

We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, there is no professional or other personal interest of any nature or kind in any product, service and/or company that could be construed as influencing the position presented in the manuscript.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Fig. S1

RNA expression levels of the five key genes and their relationship with overall survival in PDAC. (A) Relationship of the five key genes with overall survival in PDAC. (B) RNA expression of the five key genes in PDAC. PDAC, pancreatic ductal adenocarcinoma

Additional file 2: Fig. S2.

Influence of MYEOV and DDX60L in the cell cycle of PDAC. DDX60L, DExD/H-box 60 like; MYEOV, myeloma overexpressed; PDAC, pancreatic ductal adenocarcinoma

Additional file 3: Supplementary Table S1.

The qRT-PCR primers and siRNAs used in this study. qRT-PCR, quantitative real-time-polymerase chain reaction

Additional file 4: Supplementary Table S2.

DEGs of human PDAC and adjacent non-tumor tissues identified using our RNA-seq data (2A), as well as the GSE15471 and GSE16515 (2B), and GSE28735 (2C) data sets. DEGs, differentially-expressed genes; RNA-seq, RNA sequencing

Additional file 5: Supplementary Table S3.

DAVID-KEGG analyses of upregulated genes in our RNA-seq and published microarray data. DAVID-KEGG, Database for Annotation, Visualization and Integrated Discovery-Kyoto Encyclopedia of Genes and Genomes; RNA-seq, RNA sequencing

Additional file 6: Supplementary Table S4.

Venn diagram of overlapped DEGs in the microarray and RNA-seq platform analyses of PDAC and the top 100 survival OS-related genes of PDAC. (4A) Overlapped DEGs in the microarray platform of PDAC using cutoff criteria of |log2(FC)| > 1 and p < 0.0001. (4B) Overlapped DEGs in our RNA-seq data and three published microarray data. (4C) Top 100 OS-related genes of PDAC. DEGs, differentially-expressed genes; OS, overall survival; PDAC, pancreatic ductal adenocarcinoma; RNA-seq, RNA sequencing

Additional file 7: Supplementary Table S5.

DAVID analysis of upregulated and downregulated genes in siMYEOV- (5A) and siDDX60L- (5B) transfected MiaPaCa2 and PANC-1 cells. DAVID, Database for Annotation, Visualization and Integrated Discovery; DDX60L, DExD/H-box 60 like; MYEOV, myeloma overexpressed

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Wu, H., Tian, W., Tai, X. et al. Identification and functional analysis of novel oncogene DDX60L in pancreatic ductal adenocarcinoma. BMC Genomics 22, 833 (2021). https://doi.org/10.1186/s12864-021-08137-5

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-021-08137-5

Keywords