Identification of Novel Molecular Therapeutic Targets and Their Potential Prognostic Biomarkers Among Kinesin Superfamily of Proteins in Pancreatic Ductal Adenocarcinoma

Background Kinesin superfamily of proteins (KIFs) has been broadly reported to play an indispensable role in the biological process. Recently, emerging evidence reveals its oncogenic role in various cancers. However, the prognostic, oncological, and immunological values of KIFs have not been comprehensively explored in pancreatic ductal adenocarcinoma (PDAC) patients. We aimed to illustrate the relationship between KIFs and pancreatic ductal adenocarcinoma by using bioinformatical analysis. Methods We use GEPIA, Oncomine datasets, cBioPortal, LOGpc, TIMER, and STRING bioinformatics tools and web servers to investigate the aberrant expression, prognostic values, and oncogenic role of KIFs. The two-gene prognostic model and the correlation between KIFs and KRAS and TP53 mutation were performed using an R-based computational framework. Results Our results demonstrated that KIFC1/2C/4A/11/14/15/18A/18B/20B/23 (we name it prognosis-related KIFs) were upregulated and associated with unfavorable clinical outcome in pancreatic cancer patients. KIF21B overexpression is associated with better clinical outcome. The KIFC1/2C/4A/11/14/15/18A/18B/20B/23 profiles were significantly increased compared to grade 1 and grade 2/3. Besides, KIFC1/2C/4A/11/14/15/18A/18B/20B/23 was significantly associated with the mutation status of KRAS and TP53.Notably, most prognosis-related KIFs have strong correlations with tumor growth and myeloid-derived suppressor cells infiltration (MDSCs). A prognostic signature based on KIF20B and KIF21B showed a reliable predictive performance. Receiver operating characteristic (ROC) curve was employed to assess the predictive power of two-gene signature. Consequently, the gene set enrichment analysis (GSEA) showed that KIF20B and KIF21B’s overexpression was associated with the immunological and oncogenic pathway activation in pancreatic cancer. Finally, real-time quantitative PCR (RT-qPCR) was utilized to investigate the expression pattern of KIF20B and KIF21B in pancreatic cancer cell lines and normal pancreatic cell. Conclusions Knowledge of the expression level of the KIFs may provide novel therapeutic molecular targets and potential prognostic biomarkers to pancreatic cancer patients.

Methods: We use GEPIA, Oncomine datasets, cBioPortal, LOGpc, TIMER, and STRING bioinformatics tools and web servers to investigate the aberrant expression, prognostic values, and oncogenic role of KIFs. The two-gene prognostic model and the correlation between KIFs and KRAS and TP53 mutation were performed using an R-based computational framework.

INTRODUCTION
Pancreatic adenocarcinoma (PDAC) accounts for 85% of all pancreatic cancer, and it is one of the most lethal digestive duct malignancies with an increasing incidence in recent years (1). It is estimated to be the second deadliest cancer in the United States by 2025 (2). Surgical resection offers only possibility for a cure at present. However, most of the patients have lost the opportunity for surgery when they were diagnosed due to the inability of early detecting of tumors, high malignancy, and high risk of metastasis. At present, the principal treatment option for advanced patients is gemcitabine-based chemotherapy, but it only has a limited effect on long-term survival. Although immunotherapy has shown unprecedented response rates in some recalcitrant cancers, the clinic benefit of immunotherapy in PDAC such as PD-1/L1 and CTLA-4 is still dismal (3). Hence, uncovering more potential functions of known molecules will provide new therapeutic targets and indicators of prognosis.
Kinesin superfamily of proteins (KIFs) are known as molecular motors with microtubule (MT) class-binding protein superfamily, which plays a crucial role in biological processes, including cell division (4), intracellular transport (5), microtubule stabilization (6), and microtubule depolymerizers (6). KIFs include 14 families (kinesin-1 to kinesin-14) and 45 members, which exert various biological functions in our body, especially in the brain (7)(8)(9). For instance, the kinesin 4 family motor KIF4A is required for activitydependent neuronal survival (10). KIF20A was reported to play an indispensable role in cell division and Golgi-derived vesicle transportation (11). In addition, KIF21A and KIF26A are responsible for microtubule stabilizers (12,13), while KIF2A and KIF19A are implicated in microtubule depolymerizers (14). In the past few years, numerous studies have indicated that kinesin's aberrant expression is involved in development and progress many kinds of cancers (15)(16)(17)(18)(19). As a member of KIFs, KIF18A can independently predict unfavorable prognosis in lung adenocarcinoma; overexpression of KIF18A can promote cell proliferation and inhibit apoptosis (20). Moreover, KIF21B is upregulated in NSCLC and acts as an oncogene, which promotes the growth and metastasis of NSCLC (16). Besides, aberrant expressions of KIF1B (15), KIF20A (11), KIFC (21), and KIF2C (22) proteins were detected in gastric cancer, which exhibits diverse mechanisms in carcinogenesis and cancer progression.
In our study, we carried out an integrated analysis to better understand the oncogenic role of KIFs in the development of pancreatic cancer. A guided flowchart of the strategy and methodology is shown in Supplementary Figure 3. We first screened KIFs genes, all of which have the characteristic of high expression in PDAC and are associated with prognosis as our candidates based on both GEPIA and Oncomine databases. Finally, 12 prognosis-related KIFs were included for further study. Genomic alterations in prognosis-related genes were analyzed to determine the mutation status of prognosis-related KIF and overall patient survival. We performed a series of correlation analyses to explore potential roles of 12 prognosisrelated KIFs. We then continued to examine the correlations between mutation status of KRAS and TP53 and the parameters, including tumor grade, immune infiltration, and cell growth. To further explore the prognosis value of prognosisrelated KIFs, we utilized statistical methods to construct a twogene prognostic model. Protein-protein interaction (PPI) of two gene products, and some additional related gene products, was also performed via GSEA computational method to investigate its precise mechanism involved in tumor progression.

Identification of Aberrantly Expressed
Kinesin Superfamily Genes in PDAC GEPIA (http://gepia.cancer-pku.cn/index.html) is an online tool used for integrating analysis of gene expression data from TCGA (23). In our study, GEPIA was employed to analyze the aberrant gene expression of 45 kinesin superfamily members in PDAC. All the results were presented as box plots by comparing TCGA pancreatic cancer and GTEx data. The significantly up-or downregulated genes are marked with "*", which means p < 0.05.
Oncomine datasets (www.oncomine.org) is a publicly accessible online cancer microarray database used to validate the expression patterns of the prognosis-related KIFs (24,25). Student's t-test was utilized to compare the transcription levels of KIFs in pancreatic cancer specimens with those in normal controls. The cutoffs of p-value and fold-change were defined as 0.05 and 1.5, respectively. superfamily genes' prognostic values. The Cox proportional hazard ratio and 95% confidence interval information were also included in the survival plot.
LOGpc (http://bioinfo.henu.edu.cn/DatabaseList.jsp) is an online server that encompasses 209 expression datasets for survival analysis. It provides 27 types of malignant tumors for 31,310 cancer patients. Its patient samples were mainly derived from TCGA and GEO cohorts (26). The LOGpc datasets were then performed by combined overall survival analysis of prognosisrelated kinesin superfamily genes to reaffirm the genes' predictive values in PDAC. When the p-value of the combined overall survival is <0.05, it is regarded as a significant gene.

Tumor Grade Correlation With Prognosis-Related KIFs
UALCAN (http://ualcan.path.uab.edu) is an online platform aimed to perform silicon validation of potential genes of interest and clinical-pathological parameters, such as tumor grade (27). Box-whisker plot was then employed to virtualize the correlation between the tumor grade and kinesin superfamily expression patterns. The statistical significance of the two samples was accomplished using the Student's t-test, and a pvalue like to or <0.05 was deemed significant.
Kinesin Superfamily Genomic Alterations in PDAC cBioPortal (https://www.cbioportal.org/) is an online platform applying to visualize analysis and download large-scale cancer genomics datasets (28). In our study, cBioPortal was utilized to analyze genetic alterations containing alterations rate and the categories of genetic alterations of each prognosis-related KIFs in PDAC.

KRAS, TP53 mutation, and Prognosis-Related KIFs
Data of mutations about TP53 and KRAS were downloaded from cBioPortal database (http://www.cbioportal.org/) (29). We then combined these mutation data with the gene expression data of prognosis-related data downloaded from UCSC Xena (https:// xena.ucsc.edu/) to do analysis. Wilcoxon rank-sum test was employed to clarify the relation between KRAS, TP53 mutation, and kinesin superfamily members. Violin plot was archived using the R package of "ggpubr" in an R environment (R version: 4.02). A p<0.05 is regarded as statistical significance.

Prognosis-Related KIFs and Tumor Proliferation
The well-known proliferation marker ki67 was employed to reflect tumor proliferation in PDAC samples downloaded from UCSC Xena. We then assessed the association between individual prognosis-related KIFs and proliferation ability by Spearman's correlation, and |Rs| > 0.2, and FDR < 0.05 was considered as statistical significance (30). CCLE database (https://portals.broadinstitute.org/ccle/about) is applied to reaffirm the proliferation promotion role of prognosis-related KIFs in various pancreatic cell lines (31).

Immune Infiltration and Prognosis-Related KIFs
TIMER (https://cistrome.shinyapps.io/timer/) is a web resource for systematical evaluations of the clinical impact of six immune cell types: B cell, CD4 T cell, CD8 T cell, neutrophil, macrophage, and dendritic cell in diverse cancer types (32). We then utilized the correlation module of TIMER to illustrate the relationships between prognosis-related kinesin superfamily and immune cell infiltration. A purity-corrected partial Spearman's rho value and statistical significance were presented in expression scatter plots. We performed a similar analysis between the CD274 (also known as PD-L1) and kinesin superfamily. TIMER2 (http:// timer.cistrome.org/), the latest version of TIMER, was applied to clarify the relationship between immune cell infiltration and myeloid-derived suppressor cells (33).
TISIDB database (http://cis.Hku.hk/TISIDB/) is a web portal for comprehensive investigation of tumor-immune interactions, which integrate multiple heterogeneous data types (34). We analyze the correlation between prognosis-related KIFs expression and tumor-infiltrating lymphocytes in this platform.

PPI Analysis and the Similar Genes of KIF20B and KIF21B in PDAC
Cytoscape is an open-source software platform for visualizing complex networks and integrating these with any type of attribute data (35). We first retrieved 40 similar genes of KIF20B and KIF21B through expression analysis of GEPIA (http://gepia.cancer-pku.cn/detail.php?gene=ATM). Those repeated genes were deleted, and we further analyzed the protein-protein interaction of KIF20B and KIF21B and their similar genes in STRING, an online tool for performing proteinprotein interaction networks analysis. Finally, the interaction network was then imported to Cytoscape to virtualize proteinprotein interaction.

The Construction of Prognostic Model
The TCGA expression data and clinic data were downloaded from UCSC Xena, datasets that collect public databases, including TCGA, TARGET, GTEx, ICGC, and CCLE (36). To ensure the quality of our analysis, the annotations from PanCanAtlas Publications were employed to filter the unqualified samples. Eight neuroendocrine samples, 11 samples <1% neoplastic cellularity, 2 IPMN, 1 acinar cell carcinoma sample, 1 systemic treatment given to the prior/other malignancy sample, and 1 sample arise from ampulla were excluded from this study. The patients whose overall survival time <1 month were excluded from our further investigation. Then, univariate cox regression analysis was performed to identify several genes as candidate markers correlated with PDAC prognosis. LASSO regression using 10-fold crossvalidation analysis were employed to refine these genes. The refined genes were extracted to perform multivariate regression analysis. The risk score was calculated by the following formula: (1.067*KIF20B exp.) + (−0.765*KIF21B exp.).
TCGA-PAAD was divided into low-and high-risk groups according to the median risk score. The survival differences between two groups were compared and visualized with the survival status plot, risk heatmap, and K-M plotter. AUC of the 95% confidence interval was calculated based on the ROC curve. The accuracy and specificity of the two-gene model were evaluated. Diagnosis of age, gender, histological grade, stage, and risk score were all included for univariate and multivariate Cox regression analysis, which determined that risk score based on KIF20B and KIF21B was the independent prognostic factor for PDAC.

Gene Set Enrichment Analysis
To further explore underlying mechanisms of the two prognosisrelated genes, GSEA (https://www.gsea-msigdb.org/gsea/index. jsp) was performed in this study (37). We then divided patient samples into two groups with high and low expression levels, respectively, according to median expression value of prognosisrelated genes.  The amplification reaction included the following steps: 95°C for 1 min, followed by 39 cycles of 95°C for 20 s and 60°C for 1 min. GAPDH was used as an internal control for mRNA, and the relative expression level of mRNAs was calculated by the relative quantification (2 −DDCt ) method.

KIFs and Patient Survival
We then investigated the association between expression of KIFs genes and patient survival using GEPIA tool. It was noticed that poorer patient survival is significantly associated with the high expression of KIFC1 (HR = 1. We also conducted a similar analysis using the LOGpc tool, an online web tool integrating multiple GEO and TCGA datasets, to reaffirm the association between high expression of KIFs genes and overall patient survival. A similar association between high expression of KIF23, KIF2C, KIF18B, KIF20A, KIF4A, KIFC1, KIF18A, KIF14, KIF11, KIF15, and KIF20B and shorter patient survival was conformed (Supplementary Figure 2). However, the association between KIF13B and KIF21B and patient survival does not seem significant. Next, we selected genes upregulated in PDAC and their association with prognosis in patients. KIFC1/ 2C/4A/11/14/15/18A/18B/KIF20A/20B/23 expression pattern was included in our further study. The KIF21B is a protective factor detected when using GEPIA tool, so it will be also included in our further analysis. We have named the 12-genes expression pattern as prognosis-related KIFs.

Correlation Between Prognosis-Related KIFs and Tumor Grade in PDAC
Given that kinesin superfamily may be involved in the development of PDAC, we then assessed the correlation of kinesin superfamily expression and tumor grade in PDAC. We found that KIFs are significantly correlated with different grades of tumor, including KIFC1/2C/4A/11/14/15/18A/18B/ 20A/20B/23. The expression of KIFC1/2C/4A/11/14/15/18A/ 18B/20A/20B/23 increased in the comparison between grade 1 and grade 2/3 ( Figure 3). We observed that the expression of KIF21B had decreased slightly between grade 1 and grade 2/3, although not statistically significant. In brief, we concluded that some of the kinesin superfamily members are closely tied to different grades of tumor and may play potential roles in tumorigenesis and development of PDAC.

Kinesin Superfamily Genomic Alterations in Pancreatic Adenocarcinoma
Since the expression levels is the premise of function, we aimed to explore the kinesin superfamily genes highly expressed in pancreatic adenocarcinoma and their association with overall patient survival. The cBioPortal tool was employed to investigate genomic alterations of prognosis-related KIFs in TCGA pancreatic adenocarcinoma samples. The mutation landscape of prognosis-related KIFs is depicted in Figure 4. Ten percent of PDAC samples (84/850) harbored at least one kinesin superfamily gene alteration event, with amplification, mutation, and deep deletion being the most frequent alteration event (5.28%, 3.48%, and 1.93%, respectively). No significant statistical difference in survival was found between the altered and unaltered groups (p = 0.421; Figure 4). We also explored the different survival rates for individual kinesin superfamily genes in which KIF14/21B alterations were associated with poorer overall survival (Supplementary Figures 3A, B). We then performed co-occurrence and mutually exclusive analysis for prognosis-related KIFs mutations. Most of the genes were cooccurrent with each other, in which KIF20B was co-occurrent with KIF2C, KIF11, KIF14, KIF15, KIF18B, and KIF20A. However, no statistical significance was observed when we performed mutually exclusive analysis among prognosis-related KIFs (Supplementary  Figure 4C).

Prognosis-Related KIFs and KRAS and TP53 Mutations
The mutations of two major driver genes, KRAS and TP53, are associated with malignant characteristics in pancreatic cancer. We then sought to clarify the relationship between the expression pattern of the kinesin superfamily and KRAS and TP53 mutations. We generated violin plots to visualize results. We identified that the expression of most kinesin superfamily members have significantly increased in those who have KRAS mutations (KIF23, KIF2C, KIF18B, KIF20A, KIF4A, KIFC1, KIF18A, KIF14, KIF11, KIF15, and KIF20B) ( Figure 5A). On the contrary, KIF21B has significantly decreased in patients who have KRAS mutations. Only KIF20B did not harbor any significant difference between KRAS mutation and wild-type groups. A similar analysis was conducted to illustrate the relationship between the kinesin superfamily expression pattern and TP53 mutation. We noted that the expression of KIF23, KIF2C, KIF20A, KIF4A, KIF18A, KIF11, and KIF15 were lower in the TP53 mutation group ( Figure 5B). However, no significant expression pattern was detected between KIF14, KIF18B, KIF20B, KIFC1, and TP53 mutation. Interestingly, we noticed an increased expression between KIF21B and TP53 mutation, indicating that it may be a protective indicator in pancreatic cancer.

Correlation Between Prognosis-Related KIFs and Tumor Proliferation
To characterize the functional roles of the kinesin superfamily in cell proliferation, we calculated the correlation between prognosis-related KIFs and the well-known proliferation marker ki67 in PDAC. We identified a total of 11 significant associations, all of which were positive, suggesting the tumorprompting role of prognosis-related KIFs. Multiple prognosisrelated kinesin superfamily genes showed a strong correlation with Mki67, with KIF11, KIF4A, and KIF18B ranking the highest correlation coefficients ( Figure 6). To further confirm the functional roles of prognosis-related kinesin superfamily genes in cell proliferation, we performed expression profiling analysis for 41 pancreatic cancer cell lines from CCLE and observed similar expression patterns ( Table 2). To our surprise, KIF21B was negatively correlated with cell proliferation in the 41 pancreatic cancer cell lines, with the p-value of KIF21B being 0.0702. To conclude, we found that prognosis-related KIFs may play an indispensable role in cell proliferation.

Correlation Between Prognosis-Related KIFs and Immune Infiltration
To further explore the roles of prognosis-related KIFs in cellenvironment interaction, we performed a correlation of prognosis-related kinesin superfamily genes and ICB therapyrelated genes CD274 [also known as programmed death-ligand 1 (PDL1)] was also performed to assess the possible roles of those genes in the immunotherapy of ICB in PDAC. KIF23/2C/20A/ 4A/18A/14/11/15/20B were positively related to CD274, while KIF18B/C1/13B have no apparent relevance. Our findings may uncover a possible role of kinesin superfamily genes in ICB therapy ( Figure 7C).

The Construction of a Two-Gene Prognostic Model
To improve the prognostic ability of prognosis-related KIFs genes, a signature was constructed based on TCGA-PAAD cohort. One hundred fifty-four pancreatic cancer patients with prognosis information were enrolled in this study. According to univariate cox regression analysis, 11 genes (KIFC1 KIF2C KIF4A KIF11 KIF14 KIF18A KIF18B KIF20A KIF20B KIF21B KIF23) were correlated with patient overall survival. The hazard ratio and pvalue are shown in Table 3. Next, least absolute shrinkage and selection operator (LASSO) regression was carried out using 10fold cross-validation on 11 selected genes ( Figures 9A, B). Then, five genes (KIF20A, KIF4A, KIF14, 20B, and KIF21B) produced by LASSO regression were then subjected to multivariate cox regression analysis. Ultimately, KIF20B and KIF21B were used to build a two-gene panel for predicting the survival of pancreatic cancer patients. Risk score was calculated as follows: (1.067*KIF20B exp.) + (−0.765*KIF21B exp.) According to the FIGURE 6 | Correlation between prognosis-related KIFs expression and tumor growth. Scatter plot was presented to display the correlation between prognosisrelated KIFs expression and Mki67. Pearson test was employed to evaluate the correlation coefficients. Pearson r > 0 mean positive correlation, Pearson r <0 means negative correlation. Pearson r = 0 means no correlation between two genes. p < 0.05 is considered to be statistically different. median of the risk score, we divided patients into a high-and a low-risk group. The distribution of survival status and expression profile of KIF20B and KIF21B between subgroups is presented in Figure 9E. Kaplan-Meier survival analysis revealed that patients with high scores have a dismal overall survival, indicating that the two-gene scores may be important indicators for pancreatic patients' survival ( Figure 9D). What is more, receiver operating characteristic (ROC) curves were utilized to evaluate the predictive power, and area under the curves (AUCs) for 1-, 2-, and 3-year OS were 0.730, 0.665, and 0.669, respectively ( Figure 9C).

Identification of Two-Gene Model as an Independent Prognostic Factor
To identify whether the risk score could serve as the independent prognostic factor, TCGA-PAAD cohort with complete clinical information was used to perform univariate Cox regression analysis and multivariate Cox regression analysis. The results demonstrated that the two-gene signature was an independent prognostic factor (p < 0.001). The hazard ratios for OS were 1.776 and 1.547, respectively ( Figure 10).

The Protein-Protein Interaction of KIF20B, KIF21B, and Their Similar Proteins in PDAC
After knowing the possible roles of KIF20B, KIF21B genes, we began to explore their similar genes in PDAC, which may play synergistic roles in PDAC development. The results of PPI networks demonstrated that those similar genes were closely tied to each other, which affirmed a synergistic function of comparable genes ( Figure 11A). Gene Ontology (GO) analysis showed that those similar genes were significantly associated with cell division, cell cycle regulation, microtubule cytoskeleton organization, etc. (Figures 11B-D). Our results revealed that KIF20B and KIF21B and their similar genes play an essential role in cellular biological function.

Gene Set Enrichment Analysis of KIF21B and KIF20B in PDAC
Gene set enrichment analysis was carried out to explore the possible mechanisms of different KIF20B and KIF21B expression levels affected clinical prognosis in PDAC patients. The MSigDB C2 (the curated gene sets), C5 (G.O. gene sets), C6 (oncogenic gene sets), and C7 (immunologic signature) were analyzed in our study. Enrichment of C2 demonstrated that high expression of KIF20B was involved in P53 signaling pathway, Hedgehog signaling pathway, and ERBB signaling pathway ( Figure 12A). Enrichment of C5 showed that high expression of KIF20B was involved in DNA replication initiation, regulation of chromosome segregation, regulation of DNA replication, sister chromatid segregation ( Figure 12B). Enrichment of C6 showed that high expression of KIF20B was involved in various oncogene signatures such as EGFR, VEGF, and MYC ( Figure 12C). Enrichment of C2 demonstrated that increased expression of KIF21B was involved in the JAK-STAT signaling pathway, chemokine signaling pathway, T receptor signaling pathway, and MAPK signaling pathway ( Figure 12E). Enrichment of C5 showed that high expression of KIF21B was involved in positive regulation of cell killing, T cell receptor  complex, T cell receptor signaling pathway, WNT ( Figure 12F). Enrichment of C6 showed that high expression of KIF21B involved in various oncogene signatures such as P53, KRAS, MTOM, and WNT ( Figure 12G). Enrichment of C7 demonstrated that high expression of KIF20B and KIF21B was both involved in the immunologic process ( Figures 12D, H). Accidentally, we found that KIF2C/4A/11/14/15/16B/20A/22/23/25 were core enriched in immunologic signature (Supplementary Table 1).

The Expression Patterns of KIF20B and KIF21B in Pancreatic Cancer Cell Lines
Next, we investigated the gene expression of KIF20B and KIF21B in pancreatic cancer cell lines and normal pancreatic cell. We found relatively higher expressions of KIF20B in ASPC-1, BxPC-3, Capan-1, Capan-2, CFPAC-1, MIA PACA-2, and PANC-1 compared that in hTERT-HPNE ( Figure 13A). The expression of KIF21B was relatively higher in ASPC-1 Capan-1, Capan-2, CFPAC-1 compared to that in hTERT-HPNE ( Figure 13B). However, the expression level of KIF21B in BxPC-3 and MIA PACA-2 presented a relatively lower expression. In conclusion, KIF20B and KIF21B were highly expression in most pancreatic cancer cell lines compared to normal pancreatic cell.

DISCUSSION
Including PDAC, numerous studies have indicated that KIFs participate in biological functions such as intracellular transportation, cell division, occurrence, and development of many tumors (4,5). For example, KIF23 promotes cell proliferation in PDAC and is a potent therapeutic target (18). KIF18B promotes the proliferation of PDAC by activating the expression of CDCA8. Although some KIFs genes have been investigated in PDAC (43), there were no reports of different KIFs expression and their role in PDAC. As far as we know, we are the first to systematically investigated prognostic values, genetic variation, and tumor-promoting potential of different KIFs in PDAC. Moreover, for the first time, we demonstrate that the overexpression of KIF20B and KIF21B may promote immune suppression of PDAC by activating MDSCs. More importantly, we revealed the possible molecular mechanism involving the immunosuppression of PDAC and may provide new targets for immune therapy. Although previous study has investigated the prognostic value of Kinesin-4 family genes, our study uncovered the immune infiltration roles of prognosisrelated KIFs and comprehensively investigated the tumorpromoting potential of different prognosis-related KIFs in PDAC.
Tumor grade is one of the most important parameters used to judge the degree of malignancy. In this investigation, we evaluated the prognosis values of the 12 prognosis-associated KIFs. Interestingly, most of those genes were strikingly lower in grade 1 than those in grade 2/3 but not in grade 4. It indicated that the expression levels of KIFs could be good indicators of patient prognosis. Interestingly, the association between higher expression of KIF21B and a better clinical outcome indicated that KIF21B might be a good indicator in PDAC prognosis. The KIFs have been widely reported to play fundamental roles in intracellular transportation, such as in organelle, macromolecule, and messenger RNA (mRNA) (5). When tumor displays increased invasive and metastatic potential, the metabolic demand gets improved accordingly by upregulation of metabolic-related genes, such as KIFs. There is another explanation that insufficient samples of normal tissue and grade 4 result in their non-significance with other stages. Thus, further studies should be carried out to elucidate the potential roles of prognosis-associated KIFs in PDAC.
Based on the correlation between Ki67 and prognosisassociated KIFs, our study indicated that most of the genes are highly implicated in tumor growth, which is consistent with previous findings. Gu et al. suggested that knockdown of KIF26B suppressed breast cancer cell growth and invasion (17). KIFC1 was implicated in the progression of NSCLC by regulating the cell proliferation and cell cycle (44). KIF20B was not only related to tumor sizes and T stage but also promoted the progression of ccRCC by stimulating cell proliferation (19). GSEA results demonstrated that high expression of KIF20B was significantly associated with cell division, ERBB signaling pathway, and Hedgehog signaling pathway, which is in line with their role in regulating the cell proliferation. However, up to now, most of KIFs have not been systematically assessed regarding their potential abilities in supporting PDAC cell proliferation. Thus, our study provided deeper insight into the potential mechanisms of KIFs in pancreatic carcinogenesis.
The mutations of two major driver genes, KRAS and TP53, were associated with malignant characteristics in pancreatic cancer. The KRAS mutations lead to constant activation of KRAS and persistent stimulation of its downstream signaling pathways that drive many of the hallmarks of cancer, sustained proliferation, metabolic reprogramming, antiapoptosis, remodeling of the tumor microenvironment, evasion of the immune response, cell migration, and metastasis (45). Our study found that KIFs are significantly upregulated in patients with KRAS mutations, indicating that there are some complex interactions between KRAS signaling pathway and KIFs. Our GSEA results also showed that KIF20B was implicated in the MAPK signaling pathway, an essential downstream KRAS signaling pathway. On the contrary, most genes are significantly downregulated in patients with TP53 mutation. Interestingly, KIF21B was upregulated in patients with TP53 mutations while downregulated in patients with KRAS mutations. The results from GSEA also showed that high expression is strikingly enriched in the TP53 signaling pathway. An experimental study will be performed to explore the role of KIFs in the KRAS signaling pathway and TP53 signaling pathway during pancreatic carcinogenesis in our future work.
Despite rampaging evidence has shown the clinical benefits of immune-targeted approaches, checkpoint-based immunotherapy has failed to elicit responses in the vast majority of patients with pancreatic cancer (3). The immunological role of prognosis-related KIFs was investigated in some studies. For examples, KIF20A has shown its roles in PDAC immune therapy (46). However, studies on the immune suppression of KIFs in pancreatic cancer have not been performed to date. In this study, we explored the immune infiltration role of the 12 prognosis-related KIFs. The results indicated that 8/13 genes were positively correlated with B-cell infiltration. Increasing evidence revealed that B-cell subsets in PDAC tumor microenvironment upregulate associated immunosuppressive cytokines (most notably, IL-10, IL-18, and IL-35) and immune checkpoint ligands (particularly PD-L1), contributing to oncogenesis and the inhibition of T cell-mediated tumor immunity (47). Interestingly, we also found a positive correlation between prognosis-related KIFs and PD-L1. We also investigated the immune infiltration roles of prognosis-related KIFs. The results proved that prognosis related KIFs were correlated with immune infiltration in pan-cancer.
Myeloid-derived suppressor cells (MDSCs) are among the most important immunosuppressive cells in stroma of PDAC, which directly induces the T-cell suppression by secreting factors and indirectly by inducing tumor-cell-specific PD-L1 expression to inhibit spontaneous antitumor immunity (48). In our study, the correlation between infiltration levels of MDSCs and prognosis-associated KIFs was also explored. Most of those genes were significantly correlated with the infiltration level of MDSCs. Numerous studies have examined the pathways that regulate MDSCs. Thomas Welte et al. indicated that mTOR signaling recruits myeloid-derived suppressor cells to promote tumor initiation (49). Nikolaos Svoronos et al. found that tumorcell-independent estrogen signaling enhanced pSTAT3 activity through transcriptional upregulation of JAK2 and increased total STAT3 expression in myeloid progenitors and therefore drove disease progression through mobilization of myeloid-derived suppressor cells (50). GSEA also affirmed the correlation between KIF20B and MDSCs.
Moreover, GSEA showed that KIF20B was significantly enriched in the mTOR oncogenic and estrogen signaling pathway. Similar results showed that KIF21B was improved considerably in the chemokine signaling pathway, JAK-STAT signaling pathway, and T-cell signaling pathway. All these results indicated that the KIFs might regulate MDSCs in stroma by activating the mTOR oncogenic and estrogen signaling pathway and therefore mediate the immune suppression of PDAC. Consequently, we propose that KIF20B and KIF21B play an indispensable role in the immune suppression of PDAC, which could be promising therapeutic targets complementary to current immunotherapies.
In conclusion, we integrated analysis of the 12 prognosisassociated KIFs by utilizing multiple datasets and confirmed their involvements in the development and progression of PDAC. The expression pattern was first identified, followed by a series of correlation analysis. The results showed that 13 prognosis-associated KIFs were strongly correlated with tumor stage, immune infiltration, cell growth, and mutation status of KRAS and TP53. A two-gene prognostic model was used to predict the prognosis of PDAC based on each patient's risk score. ROC curves evaluated the predictive power of two-gene signature. PPI of two genes and their similar genes and GSEA demonstrate that KIF20B and KIF21B were both implicated in immune regulation and oncogenic processes. KIF20B and The expression pattern of KIF21B in pancreatic cancer cell lines and normal pancreatic cell. ns, not significant. "*" denote p < 0.05; "**" denote p < 0.01; "***" denote p < 0.001; "****" denote p < 0.0001.
KIF21B could be promising immunological therapeutic targets. However, there do exist some limitations in our study. The first limitation of this study is the insufficiency of validating external clinical cohort and experimental validation. Second, enough normal and stage 4 samples should be utilized to evaluate the prognostic value of KIFs.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.