Pancreatic ductal adenocarcinoma (PDAC) is one of the deadliest malignant tumors worldwide and has poor prognosis. New biomarkers can highlight biological differences among PDAC samples and help predict survival outcomes. However, new biomarkers with high sensitivity and specificity have not been discovered to date. Although various potential biomarkers of PDAC have been predicted using bioinformatics methods, they have not been identified in large samples. DDX31, a member of the DEAD box protein family, participates in cellular processes involving RNA secondary structure changes. However, the functions of DDX31 in PDAC remain to be elucidated.

The key gene DDX31 was identified by using a combination of a risk model and weighted gene co-expression network analysis. The biological functions of DDX31 in PDAC were investigated through bioinformatics analysis and in vitro experiments. The results of gene set enrichment analysis showed that DDX31 was correlated with cell invasive metastasis and proliferation. We found that the inhibition of DDX31 inhibited the invasive migration of PDAC cells. Survival analysis showed that DDX31 overexpression predicted poor prognosis in patients with PDAC.

The findings of our research indicate that DDX31 may be a potential prognostic biomarker for PDAC. The inhibition of DDX31 may be a potential way to treat PDAC. Combining DDX31 with other well-known biomarkers could be used to predict functional outcomes in patients with PDAC.

The 3-level mRNA expression data and corresponding phenotypic data of 160 primary PC specimens were downloaded from TCGA database. Gene expression was transformed by log2 (normalized RSEM count 1). Standardized RNA data for GSE129455 and GSE62452, GSE102238 come from the GEO database. GSE129455 contains 7500 cells of four KPC mouse samples (Elyada et al., 2019), GSE62452 contains 69 PDAC samples (Yang et al., 2016), and GSE102238 contains 50 PDAC samples (Yang et al., 2020). The GSE62452 and GSE102238 datasets were used as validation sets. Quality control was implemented using the relative expression (RLE) and standardized scale-free standard error (NUSE) in the affyPLM package provided by Bioconductor (Wang et al., 2021).

DEGs were screened by limma (Ritchie et al., 2015) package. TCGA and GEO datasets contained DEGs of metastatic and non-metastatic patients as well as cancer and normal tissues. The results were displayed by volcano plot, which were drawn using ggplot2 (Ito and Murphy, 2013).

The clusterProfile (Yu et al., 2012) package was used to perform Gene Ontology (GO) enrichment analysis, including biological processes (BPs), molecular functions (MFs), and cellular components (CCs), and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000) enrichment analysis. An adjusted p value <0.05 was considered statistically significant. GSEA local software was used for enrichment analysis and visualization. The threshold of gene set enrichment analysis (GSEA) (Subramanian et al., 2005) was set to adjusted p value <0.05 and FDR <0.25 after correction. “c2. cp. kegg. v7.0. symbols. gmt” was selected as in review the reference gene set.

Lasso-Cox regression analysis was conducted using glmnet (Friedman et al., 2010) package. The expression of the selected genes was combined with multivariate Cox regression coefficient to establish a prognostic model. The risk score of each patient is equal to the sum of gene expression and regression coefficient. The median risk score was used as the threshold to evaluate the model’s prediction effect to stratified PC patients. Patients with a risk score higher than the threshold were assigned to the high-risk group, and the rest were assigned to the low-risk group. Log-rank test was used to evaluate the difference of survival rate between the two groups. In addition, on the basis of the increase in risk value, the distribution of patients’ death events is shown by the point diagram. The heat map was used to observe the expression distribution of each characteristic gene in two different risk groups. The time-dependent receiver operating characteristic (ROC) curve was used to evaluate the risk score, and the specificity and sensitivity of predicting the survival rate of PC in 1-, 3-, and 5-year follow-up were observed. TNM staging (III–IV/I–II) and histological grade (4/3/2/1) were transformed into classification variables. The clinical features and risk scores of PDAC patients in TCGA cohort were analyzed by univariate and multivariate Cox regression analysis. The clinical features and risk scores were used to establish a nomogram model. The calibration analysis (Guy et al., 2019) was performed aligned with nomogram. Decision curve analysis (DCA) (Van Calster et al., 2018) was conducted to verify the risk model.

WGCNA (Guy et al., 2019) was conducted using the WGCNA package. It is a systematic way for effectively acquiring the expression patterns of multiple genes in different samples, which can obtain a gene cluster with the same expression pattern.

We were allowed by the Ethics Committee of Tianjin Cancer Institute and Hospital (Tianjin, China) to acquire the paraffin sections of four patients. A total of 86 patients received radical pancreaticoduodenectomy from January 2017 to September 2020. We randomly selected clinical patients who had been diagnosed as TNM3, TNM2, and TNM1. IHC of DDX31 was performed on tumor tissues of patients and matched normal pancreas tissues. Fourteen cases of fresh clinical tumor and normal pancreas tissue samples were collected to identify the expression of DDX31 by Western blot analysis. IHC analysis of the PDAC tissue for DDX31 (NOVUS; NBP1-21322, 1:400) was performed using a DAB substrate kit (ORIGENE, ZLI-9019). Three representative images of PC tissues IHC stain (100×, 200× magnification) and representative paired normal and tumor tissues IHC stain (100×, 200× magnification) were evaluated under a light microscope.

Human PC cell lines (BxPC-3, MIA-PaCa2, SW 1990, and L3.7) were purchased from the ATCC (Rockville, MD). All PC cell lines were cultured in a 5% CO2 incubator at 37°C. PC cell lines were cultured in RPMI-1640 medium and DMEM (GIBCO) with 10% fetal bovine serum (FBS).

DDX31 overexpression in PC cell lines and lentivirus-mediated plasmid was conducted using the pCDH-cDNA system (Biosettia) following the manufacturer’s instructions. Lentiviruses were produced in 293T cells for the stable transfection of cell lines. Human DDX31 cDNA was cloned into a pCDH plasmid expression vector (pCDH-DDX31), and the pCDH vector was used as the control. Stable cell lines were generated using puromycin. The overexpression efficiency was confirmed by Western blot analysis. Stable knockdown PC cell lines and shRNA were designed using http://biosettia.com/support/shrna-designer. PLVi-shRNA-bsd vectors were purchased from Biasatti. Three shRNA sequences for DDX31 were synthesized and cloned into the plasmid. Detailed information of the shRNA sequence for DDX31 is listed in Supplementary Table S1. The most effective one was used for the next experiments. The most effective shRNA was identified by Western blot analysis.

Five-week-old female nude NU/NU mice were purchased from SiPeiFu Biotechnology Co. All mice were maintained in a barrier facility on HEPA-filtered racks. All animal studies were conducted under an approved protocol (Zhao et al., 2020). Tumor cells were harvested by trypsinization, washed with ice PBS, and resuspended at 1 × 10⁷ cells per milliliter in PBS. Subsequently, 1 × 10⁶ cells were used to establish every subcutaneous xenotransplant tumor model of human PC in nude mice. In the log phase, BxPC-3 was implanted subcutaneously in nude mice and observed three times a week.

Statistical analysis was performed with GraphPad Prism version 8.0 (San Diego, CA, United States), R software 3.4.0.3, and SPSS version 26.0 (IBM SPSS, Armonk, NY, United States). Lasso-Cox regression analysis was performed using the R software package glmnet. In addition, the survdiff (Tan et al., 2020) function in the survival package was used for the logarithmic rank test. Time-dependent ROC was analyzed by timeROC (Blanche et al., 2013) package. Volcano plot was drawn with ggplot2 package. The establishment and application of line graph were realized with package rms (Sheng Zhang et al., 2019). The enrichment analysis was conducted with the R package cluster filer. Each experiment was conducted in triplicate, and data were presented as the mean ± SD unless otherwise stated. The variance between groups was statistically compared. Student’s t-test was conducted to compare mean values. Correlations between DDX31 expression level and patients’ survival time after surgery were determined by Kaplan–Meier method. The categorical data were analyzed by Chi-square test. *p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001 indicated significant differences, and NS meant nonsignificant.

Workflow of identification of DEGs among the sequencing data of TCGA, KPC-scRNAseq, and GSE36668 in PC was shown in Figure 1A. The concrete work flowchart was presented in Supplementary Figure S1. TCGA and GSE36668 datasets were divided into “with metastasis” and “without metastasis” group; KPC mouse single-cell transcriptome sequencing datasets were divided into “ductal cells” and “normal cells” group; The flowchart showed our study aimed to screen the differentially expressed genes based on screening conditions that are highly expressed and prone to metastasis in cancer. The intersection of upregulated DEGs among the sequencing data of TCGA, KPC-scRNAseq, and GSE36668 in PC was composed of 109 genes (Figure 1B), and the DEGs are shown in the volcano plot in Supplementary Figures S2A–C. The unsupervised clustering viable cells from KPC-scRNAseq were presented in Supplementary Figures S2D,E. All the differentially expressed genes were showed in Supplementary Table S1, and the results indicated that 104 genes such as ADAMTS2, BIRC5, MYC, ID1, VIM, CDH2, SMAD7, AXIN2, FAP, CAP1, SNAI1, CD59, ANXA8, YAP1, DUSP4, SOX2, KRT9, KRT13 et al. were enriched in “metastasis group” and had high expression in pancreatic ductal adenocarcinoma. Part of these genes have been reported in basic studies of pancreatic cancer and are closely related to the invasion and metastasis of pancreatic cancer according to the literature, including MYC (Sodir et al., 2020), CDH2 (Shintani et al., 2008), SNAI1 (Dai et al., 2017), YAP1, SOX2 (Zhang et al., 2017). Besides, these genes, ADAMTS2 (Dey et al., 2020), KRT9 (Andolino et al., 2018), KRT13 (Nguyen et al., 2021), LY6D (Barros-Silva et al., 2018), had been reported in a variety of cancers. GO enrichment analysis (Table 1) was performed on the intersections. The results showed that the upregulated DEGs were mainly enriched in BPs, such as Cadherin Binding involved in cell-cell Adhesion, Ephrin Receptor Binding, and Proline-Rich Region Binding, and CCs, such as tight junctions, Dendritic Shaft, and Cell-Cell Junction. Moreover, they mainly play the molecular function (MF) of Negative Regulation of transferase Activity and Regulation of Protein Tyrosine Kinase Activity (Figures 1C,D).

Lasso is another data dimensionality reduction method, which is not only suitable for linear cases, but also suitable for nonlinear cases. Lasso selects the variables of sample data based on the penalty method. By compressing the original coefficients, the originally small coefficients are directly compressed to 0, so that the variables corresponding to these coefficients are regarded as non significant variables and the non-significant variables are directly discarded. In our study, Lasso-Cox regression analysis was conducted on the selected DEGs according to the coefficients and log(λ) value. Four variables (including CAP1, CD59, YAP1, and DDX31) were screened through the change of regression parameters and the final regression to zero variables (Figures 2A,B). Lasso regression showed that the four variables contributed significantly to the penalty coefficient of pancreatic cancer patients’ outcomes, which could be used as candidate variables for further analysis. In order to further determine whether these four variables affect the survival of patients with pancreatic cancer, univariate cox regression analysis and multivariate cox regression analysis were conducted on these four variables and clinical factors (Figure 2C). These clinical factors contained T stage (T3&T4 vs T1&T2), N stage (N1 vs N0) and M stage (M1 vs M0), which were often used to reflect important clinical risk level for tumor progression. The results were summarily displayed in the forest plot, including CAP1 (HR = 1.526, p = 0.049), CD59 (HR = 1.758, p = 0.009), DDX31 (HR = 1.517, p = 0.048), YAP1 (HR = 1.904, p = 0.003), T stage (HR = 2.023, p = 0.002), N stage (HR = 1.957, p = 0.003), and M stage (HR = 0.756, p = 0.701) in univariate analysis and YAP1 (HR = 1.831, p = 0.005), CD59 (HR = 1.669, p = 0.024), DDX31 (HR = 1.229, p = 0.0453) T stage (HR = 1.280, p = 0.0466), and N stage (HR = 1.580, p = 0.041) in multivariate analysis (Figure 2C). The T, N, M stage can reflect the ability of tumor invasion and metastasis, progression and survival of patients. It is commonly used to assess the survival status of cancer patients. Thus, we found that these four parameters combined with T, N, M staging may be a pathogenic factor in pancreatic cancer. Therefore, we used these four variables to jointly construct a prediction model to assess the invasion and metastasis and survival status of patients with pancreatic cancer.

We still found in the analysis results that a small number of patients with high expression of DDX31 had HR less than 1 (DDX31′s HR value partially crossed 1), but the analysis results of most patients showed that it might be a risk factor of pancreatic cancer patients. But we still need to explore the reasons for its low risk in subsequent experiments.

A metastasis-related prognostic model was constructed on the basis of the four selected prognostic genes (CAP1, CD59, YAP1, DDX31). In our model, the probability of survival of patients was positively correlated with the risk score accompanying with the expression of CAP1, CD59, YAP1, and DDX31 (Figure 3A). the mRNA expression of CAP1, CD59, YAP1 and DDX31 was elevated obviously in high risk groups, and ultimately the high-risk group had more deaths than low-risk groups. A total of 178 PDAC samples in TCGA cohort study were used to evaluate the risk model. The results showed that the survival rate of the high-risk group was significantly poorer than that of the low-risk group (p < 0.0001), and the median survival time was significantly shortened (Figure 3B). In clinical epidemiological research, disease status and markers change over time (time-to-event outcomes). The early disease-free individuals may have a late onset due to the long follow-up time, and their markers may change from baseline during the follow-up. If the traditional ROC is used, the disease status or the time dependence of markers will be ignored. At this time, the time-dependent ROC is more appropriate. So, we performed a time-dependent ROC diagnosis and prediction analysis. We found that the 1-year survival rate prediction rate (0.602), 3-year survival rate prediction rate (0.745), and 5-year survival rate prediction rate (0.819) of the model was reliable, which were more than 0.55, further proving the effectiveness of the model (Figure 3C). In order to further evaluate the prediction performance of the risk model, we selected different clinical parameters as subgroups between high- and low-risk groups to verify the model. The median overall survival (OS) of the high-risk group was significantly shorter than that of the low-risk group in clinical stage (stage I and stage II, HR = 1.80, p = 0.007, Figure 3D), pathological stage (stage I and stage II, HR = 1.62, p = 0.026, Figure 3E), residual tumor recurrence (R0 and R1, HR = 2.22, p = 0.001, Figure 3F), and primary tumor treatment efficacy (PD and CR, HR = 1.83, p = 0.018, Figure 3G). In our model, the higher the risk score, the earlier the patient event occurred. The above results showed that we had established an effective metastasis-related model in PC. To exclude the possibility of overfitting of the model in TCGA, the risk model was verified in the GEO datasets. We analyzed the survival rate in the risk model between the high- and low-risk groups in the GEO datasets. The results showed that the probability of survival in the high-risk group was poor (HR = 2.04, p = 0.019), and the median survival time in the high-risk group was significantly shorter than that in the low-risk group (Supplementary Figure S3A). We also performed a time-dependent ROC diagnostic prediction analysis (Supplementary Figure S3B) and found that the 1-year survival rate prediction (0.661), 3-year survival rate prediction (0.703), and 5-year survival rate prediction (0.779) all had valuable predictive effects, which were greater than 0.55, indicating the effectiveness of the model. In addition, we carried out DCA (Supplementary Figure S3C) and found that the risk score has good decision-making effect. The nomogram has been regarded as a reasonable tool to create a easy intuitive graph for a statistical predictive model that quantifies the risk of a clinical event. The nomogram prediction model was also constructed (Supplementary Figure S3D) on the basis of risk factors and clinical characteristics. The 1-, 3-, and 5-year survival rate prediction of patients can match well. Furthermore, the prediction effect of the risk score was higher than that of clinical parameters, reflecting the prognostic value of the risk model. Calibration analysis was also conducted with risk scores. The results showed that the calibration of 1-, 3-, and 5-year survival corresponded to the expected prediction effect (Supplementary Figure S3E). The above results indicate that the risk model can be applied on different platforms. The percentage of copy number variation (CAP1, 4%; CD59, 4%; DDX31, 4%; YAP1, 1.8%) was calculated in TCGA samples (Supplementary Figures S4A,C). In addition, we further analyzed the clinical pathological factors and risk scores (p < .05, HR > 1) by univariate and multivariate Cox regression analysis (Supplementary Figures S3F,G), which further showed the prognostic value of the risk model.

WGCNA is a integrated way for effectively acquiring the expression patterns of multiple genes in different samples, which can obtain a gene group with the same expression pattern. The association between modules and phenotype of samples such as clinical characteristics can be studied. A total of 170 samples with clinical characteristics were included in WGCNA (Supplementary Figure S5A). In this study, the power of β = 9 (scale-free R2 = 0.8Figure 50) was selected as the soft threshold to establish the scale-free network (Supplementary Figures S5B,C). Therefore, three co-expressed modules were identified after removing the gray modules by combined dynamic tree cutting (Supplementary Figure S5D). The TOM was mapped to 2 Figures, 544 genes selected in the analysis, indicating that each module was independently verified (Supplementary Figure S5F). We found that the blue module containing DDX31 was correlated with PDAC liver metastasis (R2 = 0.38, p = 0.04). A scatter plot was mapped between GS and MM (blue module and liver metastasis) in PDAC patients (correlation = 0.470, p = 0.70). Therefore, DDX31 was considered to be one of the, critical genes in the blue module (Figure 4A). DDX31 was identified as the common target gene between the risk model and weighted co-expressed network through Venn diagram. Then, we determined DDX31 is a prognostic target gene (Figure 4B). We compared the probability of survival between patients with high DDX31 expression and patients with low DDX31 expression. The results showed that the OS of patients with high DDX31 expression was significantly decreased, and the median survival time was significantly shortened (Figure 4C). ROC diagnostic prediction analysis (Figure 4D) was conducted, and the AUC was 0.920, indicating that DDX31 had a good predictive effect on PC metastasis. Then, 135 significant genes that were positively correlated with DDX31 were selected and shown in a volcano plot (Figure 4E). Then, we conducted GO analysis on these genes and found that they were mainly enriched in ribosome biogenesis, nuclear export, and rRNA metabolic process. These genes mainly play a critical role in the production of nuclear material components, signal transmission, and structural shape (Figure 4F). In addition, KEGG pathway enrichment analysis was performed, and the results showed that gap junction, extracellular structure organization, cell junction assembly, invasion and metastasis, and cell adhesion were significantly enriched (Figure 4G).

To our knowledge, no study has investigated the biological function of DDX31 in PDAC. To investigate the expression of DDX31 in PDAC specimens, we performed IHC analysis in 86 tumor tissues and paired normal pancreas tissues (Figures 5A,B). According to the results of IHC, 62.8% (54/86) of patients, the expression of DDX31 was higher in tumor tissues than in paired normal tissues. Interestingly, in this series of 86 patients, 16.7% (9/86) of patients, the expression of DDX31 was lower in tumor tissues than in paired normal tissues and the expression was the approximate level in normal and tumor tissues in 26.7% of (23/86) these patients. According to the results of IHC, we speculate that the reasons may be as follows : First, it may be the existence of heterogeneity among pancreatic cancer patients in our center ; furthermore, there may be some possible copy number variations in this part of the patients resulting in a decline in their expression levels. Moreover, of these tissues, 14 paired fresh tumor tissues and normal pancreatic tissues were selected randomly to perform Western blot analysis (Figures 5C,D). We found that the expression of DDX31 was increased in tumor tissues. The 86 tumor tissues were divided into two groups according to the expression of DDX31 (DDX31-Low and DDX31-High) (Figure 5E). According to the results of correlation analysis of DDX31 expression and patients’ clinical pathological features, we found that high DDX31 expression was strongly correlated with the tumor size (χ² = 3.9681, p = 0.0464), lymph node metastasis (χ² = 9.9801, p = 0.0016), and TNM grade (χ² = 5.8091, p = 0.0159) of PDAC patients and histological grade of tumor tissues (χ² = 4.3641, p = 0.0369) (Table 1). At the same time, the Kaplan–Meier analysis of TMA data showed that patients with high DDX31 expression had significantly lower OS rate and relapsed-free survival (RFS) rate compared with patients with low DDX31 expression (Figures 5F,G). Next, univariate and multivariate analyses of clinical follow-up data of PDAC patients were performed. The results indicated that the expression of DDX31 was negatively correlated with OS and RFS in PDAC patients (Table 2). Together, our results indicated that high expression of DDX31 might be a risk factor of poor prognosis in PDAC occurrence and progression.

To explore the potential roles of DDX31 in PDAC cells, the basic protein expression of DDX31 in all cell lines was determined (Figure 6A). We found that most PDAC cell lines showed higher protein expression of DDX31 than the immortalized normal ductal epithelial cell line HPDE6c7. MIA-PaCa2 with low endogenous expression of DDX31 and SW1990 with high endogenous expression of DDX31 were chosen to confirm the cellular location of DDX31 in PDAC cells (Supplementary Figure S5A). The results of immunofluorescence analysis revealed that DDX31 was mainly localized to the nucleus in PDAC cells. GSEA was conducted according to the expression level of DDX31 (Figure 6B). Finally, we found that the pathways that were positively correlated with DDX31 were mainly enriched in cell adhesion molecules and cell cycle pathways. Then, we performed loss and gain-of-function studies according to the results of Western blot analysis (Figure 6C). BxPC-3-vector/DDX31-OE and MIA-PaCa2-vector/DDX31-OE cell lines were established, and related functional experiments were performed to elucidate the role of DDX31 in PDAC. First, the protein expression levels of DDX31 in BxPC-3-vector/DDX31-OE and MIA-PaCa2-vector/DDX31-OE cell lines were validated by Western blot analysis. To explore the function by which DDX31 promotes PDAC cell migration and proliferation, we performed Western blot analysis with some confirmed EMT-related proteins and proliferation-related protein antibodies. As shown in Figure 6C, the proteins N-cadherin, Snail, ZEB1, Ki67, and PCNA were positively correlated with DDX31 overexpression (Figure 6C, BxPC-3, MIA-PaCa2, L). The opposite results were observed in DDX31-knockdown PDAC cells (Figure 6C, SW 1990, L3.7, R). These results revealed that DDX31 positively regulated these classical EMT and proliferation markers. To determine the role of DDX31 in cellular migration, wound healing assay and transwell migration assay were performed. Compared with the WT PDAC cells and vector control PDAC cells, the results of transwell assay showed that the migration rate of PDAC cells significantly increased upon DDX31 overexpression (Figure 6D, BxPC-3, MIA-PaCa2), and the overexpression of DDX31 also significantly promoted cells’ wound closure (Figure 6F, BxPC-3, MIA-PaCa2). Next, we focused on the effects of DDX31 knockdown on cellular mobility capacities in SW1990 and L3.7 cell lines (Figures 6E,G; the sequence of DDX31 shRNA was placed in Table 3). By contrast, the opposite results were obtained in DDX31-knockdown PDAC cell lines (Figures 6E,G, SW1990 and L3.7). According to the results of transwell assay and wound-closing procedure, the depletion of DDX31 significantly decreased the capacity of cellular migration. Furthermore, the functional verification work of DDX31 had also conducted in HPDE6c7 (normal pancreatic ductal epithelial cells) (Supplementary Figures S7).

Crystal violet staining was conducted to evaluate the clone formation ability. DDX31-overexpressing PDAC cell lines demonstrated a greater number of colonies than the vector control and WT control (Figure 7A). EdU staining assays were also performed. As shown in Figure 7B, PDAC-DDX31-OE showed a higher percentage of EdU-positive cancer cells. This finding indicated that DDX31 overexpression led to an increased cell proliferation capability in DDX31-OE PDAC cells (BxPC-3/MIA-PaCa2) compared with the WT-CTRL/pCDH-VEC control group. In addition, the overexpression of DDX31 promoted the cell viability of BxPC-3 and MIA-PaCa2 cells by CCK-8 assay (Figure 7C). To further confirm the proliferation oncogene function of DDX31, clone formation and EDU stain assay was performed. The results revealed that DDX31 downregulation suppressed the proliferation ability of DDX31-knockdown PDAC cells (SW1990/L3.7) (Figures 7D,E). However, the downregulation of DDX31 obviously inhibited the cell viability of SW1990/L3.7 cells compared with the scramble group and WT PDAC cells by CCK-8 assay (Figure 7F). To further identify the role of DDX31 in PDAC proliferation, in vivo the mouse xenograft experiment were used. Subcutaneous tumors in nude mice were formed by BxPC-3 infected with pCDH-DDX31, pCDH-VECTOR lentivirus, and WT-CTRL PDAC cell lines, and we isolated tumors after 19 and 5 days. As shown in Figure 7G, the results revealed that DDX31 overexpression significantly promoted tumor growth.

To explore how the DDX31 affects invasive metastasis and proliferation in PDAC cells, the GSEA was carried out according to the expression level of DDX31 (Figure 8A). We found that MAPK signaling pathway was significantly enriched in high expression of DDX31. And then, correlation analysis between DDX31 and the core genes of the MAPK signaling pathway were shown with heatmap (Figure 8B and Supplementary Figures S8A,B); the results indicated that DDX31 was significantly positively correlated with these core genes (MAPK9, AKT1, ELK1,- ELK4, RAF1, KRAS, EGFR, MYC, FOS, ATF2 and so on) which played vital roles in MAPK signaling pathway (Figure 8B); for further figuring out the link between DDX31 and MAPK signaling pathway; DDX31-OE PDAC cells (BxPC-3/MIA-PaCa2) and DDX31-knockdown PDAC cells (SW1990/L3.7) were chosen to performed the western blot assays; the results of western blot assays revealed that P-ERK1/2, P-P38, EMT-related markers and proliferation markers (N-cadherin, Snail, ZEB1, and PCNA) were elevated with high expression of DDX31. To further validate the fact, the inhibitor of MAPK signaling pathway (selumetinib) was used for blocking to reverse the function of DDX31; finally, we found that the proteins P-ERK1/2, P-P38, EMT-related markers and proliferation markers (N-cadherin, Snail, ZEB1, and PCNA) in PDAC cells were consistently changed with the use of inhibitor of MAPK signaling pathway. The results indicated that DDX31 exerted the function of invasive metastasis and proliferation were probably dependent on MAPK signaling pathway; for further validating the results, wound healing assay and transwell migration assay were performed. Compared with the WT PDAC cells and vector control PDAC cells, the results of transwell assay showed that the migration rate of PDAC cells did not significantly change upon DDX31 overexpression by using the inhibitor of MAPK signaling pathway (Figure 8E, BxPC-3, MIA-PaCa2), and the PDAC cells’ wound closure showed the same results (Figure 8F, BxPC-3, MIA-PaCa2). According to the results of transwell assay and wound-closing procedure, the inhibitor of MAPK signaling pathway significantly blocked the function of DDX31. Crystal violet staining was conducted to evaluate the clone formation ability. DDX31-overexpressing PDAC cell lines demonstrated no obvious difference of colonies than the vector control and WT control by using the inhibitor of MAPK signaling pathway (Figure 8G). EdU staining assays were also performed. As shown in Figure 8H, PDAC-DDX31-OE did not show a higher percentage of EdU-positive cancer cells. This finding indicated that DDX31 overexpression did not increased cell proliferation capability in DDX31-OE PDAC cells (BxPC-3/MIA-PaCa2) compared with the WT-CTRL/pCDH-VEC control group in the case of simultaneous use of the inhibitor of MAPK signaling pathway. In addition, the overexpression of DDX31 did not change the cell viability of BxPC-3 and MIA-PaCa2 cells by CCK-8 assay (Figure 8I). To further identify the role of DDX31 in PDAC proliferation, in vivo mouse xenograft experiment were used. Subcutaneous tumors in nude mice were formed by BxPC-3 infected with pCDH-DDX31, pCDH-VECTOR lentivirus, and WT-CTRL PDAC cell lines, and the inhibitor of MAPK signaling pathway was intraperitoneally injected three times a week; finally, we isolated tumors after 19 and 5 days. As shown in Figure 8J, the results revealed that DDX31 overexpression did not show significant differences of promoting tumor growth in the case of simultaneous use of the inhibitor of MAPK signaling pathway.