Identification of Spindle and Kinetochore-Associated Family Genes as Therapeutic Targets and Prognostic Biomarkers in Pancreas Ductal Adenocarcinoma Microenvironment

Aim The role of spindle and kinetochore-associated (SKA) genes in tumorigenesis and cancer progression has been widely studied. However, so far, the oncogenic involvement of SKA family genes in pancreatic cancer and their prognostic potential remain unknown. Methods Here, we carried out a meta-analysis of the differential expression of SKA genes in normal and tumor tissue. Univariate and multivariate survival analyses were done to evaluate the correlation between SKA family gene expression and pancreas ductal adenocarcinoma (PDAC) prognosis. Joint-effect and stratified survival analysis as well as nomogram analysis were used to estimate the prognostic value of genes. The underlying regulatory and biological mechanisms were identified by Gene set enrichment analysis. Interaction between SKA prognosis-related genes and immune cell infiltration was assessed using the Tumor Immune Estimation Resource tool. Results We find that SKA1–3 are highly expressed in PDAC tissues relative to non-cancer tissues. Survival analysis revealed that high expression of SKA1 and SKA3 independently indicate poor prognosis but they are not associated with relapse-free survival. The prognostic value of SKA1 and SKA3 was further confirmed by the nomogram, joint-effect, and stratified survival analysis. Analysis of underlying mechanisms reveals that these genes influence cancer-related signaling pathways, kinases, miRNA, and E2F family genes. Notably, prognosis-related genes are inversely correlated with several immune cells infiltrating levels. Conclusion We find that SKA1 and SKA3 expression correlates with prognosis and immune cell infiltration in PDAC, highlighting their potential as pancreatic cancer prognostic biomarkers.


INTRODUCTION
Pancreatic cancer is one of the most aggressive malignancies in the world and is associated with a high rate of metastasis and mortality. Its 5-year survival rate is estimated to be 5% globally, and 11.7% in China (1)(2)(3). It is estimated that in 2018, there were 458,918 new pancreatic cancer cases and 432,242 patient deaths worldwide (4). Pancreatic cancer management is mainly by surgical resection, which is thought to improve the 5-year survival to 20-30% (5). Non-etheless, surgical resections are recommended for <20% of pancreatic cancer cases since the disease is often diagnosed at the advanced stage (5)(6)(7). Therefore, strategies for early detection and treatment of pancreatic cancer are urgently needed.
It is well-known that mitosis is a common biological process in eukaryotic cells. During mitosis, the spindle ensures that sister chromatids are correctly distributed between daughter cells (8). SKA1-3 family of genes are essential for the accurate timing of late mitosis. Spindle and kinetochore-associated (SKA) complex produced by the SKA genes maintain stability of metaphase plate or spindle checkpoint silencing (9,10). Upregulation of SKA1 triggers nucleation of interphase microtubules, while depletion of the SKA1 complex results in abnormal mitosis (11,12). The association between SKA1 and cancer has been widely investigated. It has been reported that SKA1 overexpression may lead to the development of pancreatic cancer in mouse models (13). SKA1 upregulation has also been observed in multiple malignancies, including gastric, oral, bladder, non-small cell lung, hepatocellular, and prostate cancer. Elevated levels of SKA1 have been shown to promote cancer cell proliferation and influence (14)(15)(16)(17)(18)(19)(20).
Given the limited treatment options for pancreas ductal adenocarcinoma (PDAC), novel biomarkers are needed to enhance treatment and prognosis. Previous reports show that SKA family genes may influence cancer treatment response and prognosis. We therefore speculated that SKA family genes might be prognostic in PDAC. To test this possibility, we used bioinformatics to investigate their expression in PDAC and how it correlates with prognosis.

Acquisition of Public RNA-Seq and Gene Microarray Data
A schematic representation of our study outline is shown in Figure 1. Violin plots and RNA-seq data of SKA gene expression in normal tissues were obtained from Genotype-Tissue Expression (GTEx) 1 (21,22). RNA-seq raw data for 149 samples (145 PDAC tissue and 4 non-tumor pancreatic tissue) were obtained from The Cancer Genome Atlas (TCGA) 2 (23). Matched clinical information of the PDAC patients was obtained from University of California Santa Cruz Xena Platform (UCSC Xena). 3 RNA-seq raw data were normalized by DESeq package (24). Pancreas ductal adenocarcinoma SKA gene expression profile microarrays, normalized using limma package (25), were retrieved from Gene Expression Omnibus (GEO). 4 All the expression profile data were log2-transformed before further analysis.

Meta-Analysis
The following search terms were used to retrieve relevant datasets from GEO: "pancreatic cancer, " "mRNA, " and "pancreas." The inclusion criteria for PDAC microarrays were as follows: (1) it must contain both normal and tumor samples, (2) expression profile data are available, (3) samples are all human, and (4) ≥10 samples available. A schematic of the search process is shown in Supplementary Figure 1. Details of the included microarrays are shown in Supplementary Table 1. Eleven previously published datasets were retrieved (GSE71729, GSE28735, GSE15471, GSE62165, GSE62452, GSE16515, GSE32676, GSE1542, GSE74629, GSE91035, and GSE101462) (26)(27)(28)(29)(30)(31)(32)(33)(34)(35). The clinical information of GSE62452 and GSE28735 was also downloaded to explore the potential prognostic value of SKA genes. After adding the TCGA PDAC cohort, a total of 715 PDAC samples and 297 non-cancer pancreatic samples from 12 datasets were combined and meta-analysis was used to compare gene expression between PDAC and normal tissue.

Bioinformatic Analysis Using SKA Family Genes
The online tools of Gene Multiple Association Network Integration Algorithm (GeneMANIA) 5 and the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) 6 were used to structure interaction networks in gene-gene and protein-protein, respectively (36)(37)(38)(39). A co-expression matrix of SKA mRNA was constructed using the TCGA, GTEx, and GSE62452 datasets. Stratified survival analysis with multiple clinicopathological features was also used to evaluate the predictive value of the genes in PDAC patients.

Nomogram
All PDAC patient data from TCGA and its matched clinicopathological features were combined with prognosisrelated genes to construct a nomogram. Each prognosis-related gene was classified into the high or low expression groups based on median gene expression. The nomogram can calculate the total score of each patient based on existing information and predict survival probability. In addition, the nomogram can evaluate how each parameter affects the probability of survival.

Prognostic Value Evaluation and Clinical Relevance of SKA1 and SKA3
Area under the curve (AUC) of the receiver operating characteristic (ROC) curve was done to assess the accuracy of prognostic-related genes in predicting patient survival using survivalROC package. We compared the correlation between expression of prognosis-related genes and clinicopathology features of PDAC using the datasets with adequate clinical information.

Potential Mechanism and Regulatory Factors of Genes
Samples from TCGA were divided into two categories according to median expression values of prognosis-related SKA genes. Gene set enrichment analysis (GSEA) 7 was then done to explore potential biological mechanisms (40,41). Gene sets of Kyoto Encyclopedia of Genes and Genomes (KEGG) (c2.all.v6.2.symbols.gmt) and gene ontology (GO) (c5.all.v6.2.symbols.gmt) were used in this study (42,43). The criteria for statistically significance were a nominal P value < 0.05 and FDR < 0.25. (44). Next, kinase-target, transcription factortarget, and miRNA network enrichment analyses were done on the LinkedOmics database 8 online tool (45). The top five outcomes from each gene set are shown.

Analysis of Genomic Alterations and Methylation Level of Prognosis-Related Genes
To elucidate the potential regulatory mechanism of gene expression, gene mutations, and copy number variations (CNVs), data were accessed on cBioportal. 9 The calculation method of CNVs was Genomic Identification of Significant Targets in Cancer 2.0 (GISTIC2) (46). Next, DNA methylation data for the genes with matched RNA-seq expression profile (log2(count + 1)) were accessed from UCSC Xena and analyzed by Infinium Human Methylation 450 BeadChip. The methylation assessment for the genes was expressed as β-values. Correlation analysis and KM survival analysis were used to explore potential methylation sites regulating gene expression.

Exploration of Tumor Immune Infiltration
RNA-seq data were analyzed using ESTIMATE (estimation of stromal and immune cells in malignant tumor tissues using expression data) 10 algorithm to calculate cellular component in tumor tissues and estimate tumor microenvironment (TME) scores for stromal and immune cells (47). The immune cells analyzed to assess tumor immunity are B cells, CD4+ T cells, CD8+ T cells, dendritic cells (DCs), neutrophils, and macrophages. Tumor Immune Estimation Resource (TIMER) 11 was then used to evaluate infiltration by the immune cells (48). Next, the correlation between immune cell infiltration and prognosis-related gene expression in PDAC was evaluated.

Statistical Analysis
Stata (version 12.0) was used to plot forest plots and summary ROC (sROC) curve of the meta-analysis. The standard mean difference (SMD) and 95% confidence interval (CI) were used to identify differentially expressed SKA genes in cancer and noncancer tissues. Summary ROC was used to evaluate the capacity of genes to discriminate between cancer and non-cancer tissues. All statistical analysis was carried out on SPSS (version 22.0). Survival differences were identified by hazard ratio (HR) and 95% CI. Two or multiple groups of continuous variables were compared by Student t test and one-way ANOVA, respectively. All the correlation analysis method used Spearman's correlation. P value < 0.05 was considered statistically significant.

Bioinformatics Analysis
Gene and protein interaction analysis SKA1-3 using GeneMANIA and STRING, respectively, revealed a significant degree interaction, as well as protein homology (Figures 3A,B). Furthermore, there existed strong co-expression of SKA genes not only in the normal dataset of GTEx but also in tumor datasets of TCGA and GSE62452 (Figures 3C-E), especially the correlation between SKA1 and SKA3 (r = 0.62, 0.6, and 0.69).

Survival Analysis
We first analyzed TCGA dataset and removed a case lacking the survival time so that 145 eligible PDAC patients were enrolled at last. After performing KM survival analysis with the clinicopathological indicators of PDAC, the P value less than 0.05 clinical features including radical resection, radiation therapy, and targeted molecular therapy were selected out, which were further involved in the Cox regression analysis for adjustment ( Table 1). Then, the KM method was also used for SKA1-3 gene survival analysis, indicating that all the genes were significantly associated with prognosis in PDAC patients except for SKA2, and high expression suggest a shorter survival time (Figures 4A-C). The median survival time (MST) for low and high expression groups of SKA genes were SKA1 (695 vs 485 days), SKA2 (568 vs 518 days), and SKA3 (698 vs 485 days) ( Table 2). Similarly, the Cox regression analysis outcomes showed that SKA1 (adjusted P = 0.04; adjusted HR = 1.656, 95% CI: 1.024-2.677) was an independent prognostic factor for PDAC patients as well as in SKA3 (adjusted P = 0.034; adjusted HR = 1.688, 95% CI: 1.040-2.742), while it was not with SKA2 (adjusted P = 0.837; adjusted HR = 0.952, 95% CI: 0.592-1.529). Survival analysis from GEO datasets of GSE62452 and GSE28735 showed that high expression of SKA1 and SKA3 was significantly correlated with a poor OS, except for SKA1 in GSE28735 due to the limited sample (Figures 4D-I). To estimate the role of SKA genes in recurrence of PDAC, relapse-free survival (RFS) time was used to structure the KM curve, which suggests that the high expression group of SKA1 or SKA3 was not associated with the recurrence rate of PDAC compared with the low expression group (Figures 4J-L).
We further explored the effect of the gene combination on the prognosis of patients and divided samples into three groups. KM survival curves showed that group 3 (MST: 394 days) with high expression of SKA1 and SKA3 genes had the worst prognosis, while group 1 (MST: 695 days) with low expression of the two genes had the best prognosis, suggesting that their high expression represents a high risk of death ( Table 2 and

Nomogram
Patients were grouped into three categories as mentioned previously, and patient mRNA expression data with matched clinical features were used to build a nomogram. Each variable had a score that was denoted by a line length. The higher the score, the greater the effect of the gene on prognosis. This analysis reveals that that mRNA expression of SKA1 and SKA3 can significantly affect patient survival ( Figure 6C).
Subsequently, we evaluated correlation between prognosisrelated genes and clinicopathological parameters. Analysis of the GSE62452 dataset revealed SKA1 upregulation in stage II patients relative to stage I patients (Figures 7F,G). In terms of histologic grade, analysis of the GSE62452 dataset showed that SKA1 and SKA3 expression were significantly higher in G3/G4 PDAC relative to G1/G2 disease stage and SKA1 and SKA3 only showed an upregulated trend in G3 samples in TCGA (Figures 7H,I).
Moreover, relative to tumor-free survival patients, tumor-bearing patients significantly expressed high levels of SKA1 and SKA3 MST, median survival time; HR, hazard ratio; CI, confidence interval.

Gene Enrichment Analysis
Gene set enrichment analysis analysis of the potential mechanism by which prognosis-related genes influence PDAC revealed that in the high SKA1 expression group, there was significant enrichment for cell cycle-related biological processes (GO: cell division, cell cycle checkpoint, and cell cycle phase transition, Figures 8A-C) and tumor-related signaling pathways (KEGG: cell cycle, P53 signaling pathway, and DNA replication, Figures 8D-I). GO term analysis of the high SKA3 expression group revealed enrichment for processes involved in negative regulation of T cell proliferation and positive regulation of the WNT pathway (Figures 8J,K). KEGG annotation showed the high SKA3 expression group participated in tumor-related signaling pathways (pancreatic cancer, pathways in cancer, toll-like receptor pathway, and adherens junction) (Figures 8L-O). The enrichment results of network analyses on the LinkedOmics database indicated that SKA1 and SKA3 are mainly regulated by the same transcription factors (V$E2F1_Q6, V$E2F_Q4, V$E2F_Q6, V$E2F_Q4_01), networks, and kinases, including cyclindependent kinase 1 (CDK1), polo-like kinase 1 (PLK1), cyclin-dependent kinase 2 (CDK2), and aurora B kinase (AURKB). The miRNA-target network for SKA1 was related to MIR-185, MIR-512-3p, MIR-507, MIR-218, and MIR-96, while SKA3 was associated with MIR-507, MIR-119A, MIR-513, MIR-338, and MIR-137 ( Table 5).

Genomic Alterations and DNA Methylation Level of Prognosis-Related Genes
Analysis of SKA gene mutation in the 153 PDAC patients showed that 12 (7%) of them were mutation carriers ( Figure 9A). Furthermore, box plot analysis of CNV data from the PDAC patients revealed that increased copy number of the genes poorly correlates with higher expression of prognosis-related genes (Figures 9B,C). First, We found poor correlation between DNA methylation and gene expression in pure PDAC (data not show). Next, we explored whether DNA methylation status influenced gene expression in 185 pancreatic cancer samples by analyzing 11 and 15 CpG sites in the SKA1 and SKA3 DNA locus, respectively. Five of the 11 sites in the SKA1 DNA locus had high methylation with β values > 0.6, while two sites (cg18558188: β ± SD = 0.84 ± 0.04, | r| = 0.254; cg18742986: β ± SD = 0.89 ± 0.06, | r| = 0.19) exhibited significant negative correlation with gene expression, indicating that DNA methylation may regulate SKA1 expression ( Figure 9D). When samples were divided into two groups based on median methylation level, loci cg18558188 and cg18742986 methylation were not associated with PDAC OS, while higher total CpG methylation levels of SKA1 correlated with improved patient survival ( Supplementary  Figures 3A-C). All SKA3 CpG sites exhibited very low methylation (β-value < 0.4) and were not explored further ( Figure 9E).

Prognosis-Related Genes Correlated With Tumor Immune Microenvironment and Key Gene Mutations
GSEA revealed that SKA3 modulates T cell activity. Next, we explored relationship between prognosis-related genes and   tumor immunity, including immune cell infiltration and immune scores. Results from TIMER analysis show that CD8+ T cells, CD4+ T cells, and macrophages exhibit a significant negative correlation with SKA1 mRNA expression and its CNVs significantly impact immune cell infiltration levels (Figures 10A,B). In particular, SKA3 expression significantly exhibits a strongly negative correlation with all immune cell infiltration levels except for B cell. Infiltration of natural killer (NK) cells is not related to SKA gene expression ( Supplementary  Figures 4A-D). Additionally, infiltration levels for B cell and CD4+ T cell can be affected by SKA3 CNVs (Figures 10C,D). Then, scatter plot analysis was used to visualize the distribution of CD4+ T and CD8+ T cell between the groups expressing high and low levels of prognosis-related genes. This analysis revealed that CD8+ T cells were significantly fewer in the high expression group, but no difference in CD4+ T cell distribution was observed between the groups (Figures 10E,F). So, we speculated that high expression levels of SKA1 and SKA3, to some extent, may mediate tumor escape and inhibit the infiltration levels of immune cells. To test this possibility, we calculated the tumor immune scores using the RNA-seq data from the TCGA cohort. Consistently, immune score negatively correlated with SKA3 expression but did not correlate with SKA1 expression, indicating that SKA3 has a great influence on immune cell infiltration in PDAC tissues ( Figure 10G). Tumor cells have the capacity to evade clearance by macrophages through the upregulation of antiphagocytic surface proteins including cluster of differentiation 47 (CD47), programed death-ligand 1 (PD-L1), programmed cell death protein 1 (PD1), cytotoxic T-lymphocyteassociated protein 4 (CTLA4), and beta-2-microglobulin (B2M) (49)(50)(51)(52). We explored whether SKA1 and SKA3 had correlation with these antiphagocytic molecules. SKA1 and SKA3 were positively correlated with the expression of CD47, PD-L1, and B2M in TCGA dataset, respectively. Similar results were displayed in the GSE625452 dataset (Supplementary Table 2). These results may explain why high expression of SKA1 and SKA3 could result in reduced infiltration levels of immune cells in PDAC. Mutation analysis revealed that high SKA1 expression highly correlates with the KRAS mutation group but not with TP53 mutations. High SKA3 expression in the KRAS and TP53 mutant groups was significantly associated with the mutant type (Figures 10H,I).

DISCUSSION
In this study, we combined PDAC cohorts' data from GEO and TCGA and used bioinformatics to evaluate the potential prognostic value of SKA genes in PDAC. Our meta-analysis found that all SKA genes are significantly upregulated in PDAC and exhibited medium diagnostic value for PDAC. Survival analysis showed that overexpression of SKA1 and SKA3 is associated with poor prognosis. The prognostic value of the two genes was further studied by ROC, joint-effect, and stratified analysis among other strategies. Our data highlight SKA1 and SKA3 as potential therapeutic targets against PDAC. Additionally, we explored the potential mechanisms regulating SKA1 and SKA3 expression and found that DNA methylation may influence SKA1 expression and patients' OS. Possible mechanisms include cancer-and immune-related pathways, and potential regulators include cancer-related kinases, miRNA, and the E2F family. Finally, we explored the association between prognosis-related genes and tumor immune microenvironment. Biochemically, SKA1 is known to directly bind microtubules through its C-terminal domain to stimulate oligomerization (11,53). Spindle and kinetochore-associated complex deficiency is associated with chromosome congression failure and cell death (54)(55)(56). Multiple studies have reported that SKA1 promotes cancer progression in a variety of tumors (13,19,(57)(58)(59)(60)(61). An in vitro study found that SKA1 accelerates cell proliferation and cancer progression in glioma via tumor-associated signaling pathways, including cell cycle and Wnt/β-catenin signaling (62). An immunohistochemical study of 126 hepatocellular carcinoma patients revealed that SKA1 expression is significantly elevated in tumor tissues, and it can regulate the hepatocellular carcinoma cell cycle and contributes to poor prognosis (63). Qin et al. used targeted small interfering RNA to knockdown SKA1 and observed hepatocellular carcinoma cell cycle arrest in the G0/G1 phase (17). Besides, SKA1 has been reported to contribute to chemotherapy resistance in lung carcinoma through prevention of cisplatin-induced apoptosis (19). Here, we find that high SKA1 expression predicts poor PDAC prognosis and participates in cell cycle, cell cycle checkpoint, P53 signaling pathway, and DNA replication, processes that significantly correlate with cancer progression (64)(65)(66)(67). Our data also associated high SKA1 expression with advanced cancer phenotypes such as stage II, G3 + G4, and patients survived with tumor. SKA2, a novel cell cycle gene, has been proposed as a biomarker and therapeutic target against cancer (68). It has been reported that suppressed SKA2 expression triggers kinetochore fiber instability, leading to mitotic failure (9). A study by Ren et al. found that SKA2 overexpression induces epithelial-mesenchymal transition in breast cancer, promoting cell invasion and metastasis (69). SKA2 also influences proliferation, migration, and invasive capacity of gastric and lung cancer cells (70,71). While SKA2 was also highly expressed in our study, no correlation was found with PDAC prognosis. A probable reason for this may be intertumor heterogeneity, which requires additional investigation. SKA3 mediates appropriate mitotic exit by interacting with the NDC80 complex, which regulates meiotic spindle migration and anaphase spindle stability (55,72,73). SKA3 upregulation in lung adenocarcinoma cells correlates with increased metastases and tumor growth (74). Similarly, it has been reported that high SKA3 expression promotes lung cancer cell proliferation and predicts patient outcomes. Through bioinformatics analysis, Tang et al. found that SKA3 is associated with elevated susceptibility to breast cancer brain metastasis and negatively correlates with breast cancer survival (75). Our data further validate the correlation between high SKA3 expression and advanced clinical features of G3/G4. Our findings show that SKA3 may be involved in Wnt signaling, pancreatic cancer, pathways in cancer, toll-like receptor signaling pathway, and adherens junction, which are known to play a crucial role in tumor progression (76)(77)(78).
Kinases regulate various processes such as genomic stability, mitosis, and the cell cycle. CDK1, a member of the cyclindependent kinase protein family, participates in mitosis, cell cycle, cell differentiation, and somatic reprogramming (79). Previous investigators have shown that dysregulation of CDK1 leads to G2 phase arrest and promotes tumor progression, making it an ideal biomarker and therapeutic target (80,81). CDK1 binds Ndc80, thereby phosphorylating SKA3 and recruiting SKA to kinetochores to facilitate mitotic progression (H,I) The expression distribution of SKA1 and SKA2 genes in different pathological T grade and cancer status in TCGA cohort. *p < 0.05, **p < 0.01, ***p < 0.001. (82). Another study by Hou et al. showed that SKA3 interacted with CDK2 and inhibited P53 phosphorylation, thereby regulating proliferation of liver cancer cell (83). PLK1, a tumorigenic factor, is a mitotic cyclin-independent serine threonine kinase that has been proposed as a potential therapeutic target for pancreatic cancer (84,85). Another kinase, AURKA, participates in pancreatic carcinogenesis via the MAPK1/ERK2 signaling pathway (86). The same regulatory network of kinases (CDK1, CDK2, PLK1, and AURKA) was identified in our study as potential regulators of SKA1 and SKA3 in PDAC. The E2F family is part of the transcription factors that regulate gene expression. It participates in the cell cycle process, and activated E2F initiates oncogenic signaling in several cancer types (87,88). Several approaches have been developed to directly or indirectly target E2F1 with the aim of modifying malignant phenotypes of PDAC (89)(90)(91)(92). Herein, E2F1 was found to work with SKA1 and SKA3 to jointly regulate cell cycle and aggravate PDAC. MicroRNA, an endogenous small RNA with a length of about 20-24 nucleotides, has been implicated in human carcinogenesis (93). The miRNAs identified in our study have been associated with deteriorated neoplastic malignant phenotype such as proliferation, cell cycle, invasion, drug resistance, and angiogenesis (94)(95)(96)(97)(98). In fact, miR-185, miR-96, miR-218, miR-137, and miR-338 have been proven to exhibit therapeutic and prognostic value in PDAC, indicating that the SKA gene may be one of the target genes for these miRNAs (99)(100)(101)(102)(103).
A research team found that genomic alteration is a common phenomenon in human tumors (104). Copy number variations, induced by genomic rearrangement, disrupt genes and alter genetic content, causing different phenotypes. In this study, CNVs of SKA1 and SKA3 were found to be weakly correlated with gene expression. DNA methylation is a crucial epigenetic mechanism, which maintains genome stability, chromatin structure, and pluripotency in human cells. Changes in DNA methylation levels often accompany neoplasm development (105,106). The role of DNA methylation in SKA genes has not been investigated. In this preliminary study, we found that two sites in SKA1 were significantly involved with gene expression. High total DNA methylation level of SKA1 predicts favorable prognosis in PDAC. We thus hypothesize that DNA methylation may be an important regulatory mechanism for SKA1 that influences the OS of patients with PDAC. However, no DNA methylation sites were found to regulate SKA3 gene expression, which may need more work to test it. In recent years, immunotherapy has shown great promise for cancer patients, especially those with refractory cancer. The complexity of immune cells in tumor tissue influences the host biological behavior and the outcome of immunotherapy (107). The TME, composed of non-cancerous cells in and around the tumor, also performed an important role in the genomic analysis of various tumors (108). CD8+ cytotoxic T lymphocytes (CTLs) can specifically recognize major histocompatibility complex (MHC) antigens, which are widely used in targeted therapy (109). A pioneer study revealed that infiltration of CD4+/CD8+ T cells correlated with good prognosis in PDAC (110). Likewise, a study by Masugi et al. (111) measured densities of CD8+ T cells in different locations of tumor. They found that high CD8+ T cell density was significantly associated with prolonged survival of 214 patients with pancreatic cancer (111). Wu et al. also reported that the low infiltration level of CD4+ T lymphocytes was associated with poor prognosis of pancreatic cancer patients (112). GSEA suggested that SKA3 negatively regulated T cells. Thus, we assessed the association between immune cell infiltration and SKA prognosis-related genes in PDAC. Results showed that SKA1 was negatively correlated with CD8+ T cell and macrophages, and the infiltration level of CD8+ T cells was significantly lower in the high SKA1 expression group, indicating an immunosuppressive state. Similarly, the abundance of CD4+ T cells, CD8+ T cells, and macrophages was decreased in the SKA3 high expression group. There is a widespread belief that CD4+ T cells compromise the majority of T cells in pancreatic cancer and are positively associated with tumor metastasis and negatively associated with OS (113). Some subsets of CD4+ T cells may also be needed for antitumor immunity. CD4+ helper T cells may promote and maintain CTL memory, amplify T and B cells, and help CTL resist negative regulation (114). CD4+ T lymphocytes may inhibit tumor cell growth by cytolysis or by regulating the TME (115). More detailed studies are necessary to illustrate the specific role of each CD4+ T lymphocyte subset in pancreatic cancer. Macrophages, participating in the production, mobilization, activation, and regulation of immune effector cells, have at least three major functions: antigen presentation, phagocytosis, and immunomodulation (116). Cytokines such as IL-1, IL-6, TNF-α, interferon (IFN)-α/β, IL-10, IL-12, and IL-18 released from macrophages could participate in the regulation of immune/inflammatory responses. IL-12 stimulates proliferation of activated T and NK cells, enhances NK and lytic activity of CTLs, and induces IFN-γ production by T and NK cells. In addition, they produce chemokines that stimulate lymphocyte movement and regulation of migration lymphocytes from the blood to tissues (117). Similarly to macrophages, DCs have long been established as indispensable antigen-presenting cells (APCs), which act as systemic sentinels capable of responding to endogenous and exogenous "danger" signals to initiate and propagate immune responses to inciting antigens or induce immune tolerance (118,119). On sensing of appropriate cues, DCs mature and express chemokine receptors and costimulatory molecules under normal conditions. DCs promote immunity or tolerance by sampling and presenting antigens to T cells and by providing immunomodulatory signals through cell-cell contacts and cytokines (120,121). DCs are often associated with superior cross-presentation of antigens, which results in stronger CD8+ T cell immunity, and DCs can additionally support T helper 1 cell polarization of CD4+ T cells (122,123). In tumor patients, DCs acquire, process, and present tumor-associated antigens on MHC molecules and provide costimulation and soluble factors to shape T cell responses. However, a number of active mechanisms in the TME perturb DC functions, resulting in insufficient T cell activation and, potentially, the induction of T cell tolerance to tumor-associated antigens (119). Lymphocytes and macrophages decreased significantly in patients with high expression of SKA1 and/or SKA3 in this study, leading to reduction of activation of immune cells, including CD8+ T and CD4 + T cells and NK and lytic activity of CTLs, suggesting that they might have an immune-excluded phenotype where CD4+ T/CD8+ T cells were maintained in the stroma, restricting cancer immunity. Interestingly, the infiltration of cytotoxic CD8+ T cells might be modified by immunotherapy. The above mechanism may be potentially responsible for short-term survival in PDAC patients correlates with increased SKA1 and/or SKA3 expresses. Analysis of immune scores in TME yielded similar outcomes as those of immune infiltration. In addition, SKA1 and SKA3 were positively correlated with the expression of antiphagocyticrelated genes (CD47, PD-L1, and B2M). In general, these findings indicate that SKA1 and SKA3 play a significant role in the recruitment and regulation of immune-infiltrating cells in PDAC, which may eventually influence patients' survival time. Thus, we hypothesized that patients with high expression of SKA1 or SKA3 might benefit from immunotherapy than those with low expression. Further research is required to address this hypothesis.
Finally, some limitations exist in this study. First, this study was performed using retrospectively collected data, which may contain selection bias and recall bias. Secondly, our results are based on bioinformatics analysis and thus underlying biological mechanisms remain undefined. Thirdly, the protein expression levels of SKA1 and SKA3 and their involvement in the pathogenesis and progression of PDAC deserve further studies. Despite these limitations, this is the first study to reveal the potential correlation between SKA genes and PDAC tumor immune escape and comprehensively explore the prognostic value of SKA genes in patients with pancreatic cancer. Our results may be informative for future research and clinical management of PDAC patients. Despite these limitations, this is the first study to reveal the association of SKA genes with immune function regulation of PDAC patients. There have been no reports demonstrating that SKA genes could regulate the immune infiltration of PDAC.
Collectively, high expression of SKA1 and SKA3 predicts poor prognosis of PDAC and may therefore be potential biomarkers for this disease. These genes could regulate cancer-related signaling pathways and inhibit immune infiltration within the tumor in PDAC. Further prospective studies are required to verify these molecular mechanisms.

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.

ETHICS STATEMENT
Ethical review and approval was not required for the study on human participants in accordance with the local legislation and institutional requirements. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.

AUTHOR CONTRIBUTIONS
XH, YC, and YL performed the data analysis work and aided in writing the manuscript. QL and ZJ designed the study and assisted in writing the manuscript. YL and QL edited the manuscript. All authors read and approved the final manuscript.

ACKNOWLEDGMENTS
This is a short text to acknowledge the contributions of specific colleagues, institutions, or agencies that aided the efforts of the authors. We would like to thank Sheng-xin self-study network for providing technical support.