We used the term “type 2 diabetes mellitus” or “pancreatic cancer” to search for T2DM and PC patients’ gene expression profiles in Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo/) database. The obtained datasets were screened by the following criteria. First, profile information should include both case and control groups. Second, all samples should be derived from pancreas tissue. Third, pancreatic cancer was not treated with chemoradiotherapy before surgical resection. Fourth, these datasets must provide raw data that can be further analyzed. Last, GEO datasets, GSE38642, GSE20966, GSE25724, GSE91035, GSE55643, and GSE32688, were chosen for next research. In addition, the somatic mutation data and transcriptome data were acquired from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) database. HTSeq-FPKM value needs to be transformed into transcripts per million (TPM).

WGCNA is an algorithm for analyzing gene expression patterns in numerous samples. It is capable of clustering genes and constructing modules based on similar gene expression patterns, as well as analyzing the associations between modules and biological traits (16). In this study, we used the R package “WGCNA” to create the gene co-expression networks of T2DM and PC. Firstly, we built the adjacency matrix using soft-threshold b (7 for T2DM, 8 for PC) and gene–gene correlation matrix to describe the degree of association between the nodes. The adjacency matrix was then converted into the topological overlap matrix (TOM). Next, the gene hierarchical clustering dendrogram was performed to identify co-expression modules. Finally, we calculated the module eigengene (ME), as well as the correlation between ME and clinical traits, to identify clinical-related modules.

The shared genes in T2DM and PC modules with positive correlation coefficients were overlapped using Jvenn (17). ClueGO, a plug-in of Cytoscape, could classify non-redundant GO terms and visualize the functionally related genes in a clustered network (18). To understand the biological function of these shared genes, we performed a GO analysis by utilizing the ClueGO. The p-value <0.05 was identified as the significant GO term.

To verify the shared gene signatures in T2DM and PC, the protein–protein interaction (PPI) network was next created using the “MCODE” algorithm with default settings in Cytoscape. The hub clusters were further enriched by GO. Disease Ontology (DO) analysis was conducted using R package “DOSE” (19).

To confirm hub shared genes in T2DM and PC, we conducted the differentially expressed gene (DEG) analysis on validation datasets (GSE20966, GSE25724, GSE55643, and GSE32688) using the R package “limma.” p < 0.05 was considered as the cutoff value. The overlap of above shared genes in test cohorts and DEGs in validation cohorts was represented by a Venn diagram using Jvenn tool.

Based on S100A6 median expression, PC patients were categorized into high- and low-expression groups. GSEA is a method based on functional categories that could calculate the enrichment score of gene sets and discover different functional phenotypes (20). We used GSEA to compare the biological pathways between the two groups. The h.all.v7.2.symbols.gmt gene set was downloaded as a reference. Kyoto Encyclopedia of Genes and Genomes (KEGG) significant pathways with false discovery rate (FDR) <0.05 were included.

The tumor immune microenvironment (TIME) is assessed by single-sample Gene Set Enrichment Analysis (ssGSEA) method using R package “GSVA.” ssGSEA, an extension of GSEA, could calculate separate enrichment scores (ESs) for each sample (21). Geneset variation analysis (GSVA) is an algorithm that converts an expression matrix with a single gene into an expression matrix with a specific gene set and calculates its ES (22). The levels of 16 immune cells and 13 immune-related functions between high- and low-S100A6 expression groups were then compared by ssGSEA scores. ESTIMATE is a tool to compute each sample’s immune score and stromal scores (23). Spearman’s correlation was used to investigate the relationship between S100A6 and immune cell infiltrates.

The PC cells were cultured in Dulbecco’s Modified Eagle Medium (DMEM) with 10% Fetal bovine serum (FBS) at 37°C in a 5% CO₂ incubator. The cells in logarithmic phase were chosen for functional experiments. S100A6 knockdown was achieved by transfecting cells with indicated S100A6-siRNA using Lipofectamine^® 3000 reagent (Invitrogen, USA) as directed by the manufacturer’s instructions. The target sequences for siS100A6-1 and siS100A6‐2 were 5’-GCUCACCAUUGGUGCUAAG-3′ and 5′-CCUCUCUGAGUCAAAUCCATT-3′.

Total RNA was isolated from cells using TRIzol reagent (Invitrogen, America), and cDNAs were generated using HiScript Synthesis kit (Vazyme, China). Next, real‐time quantitative PCR (qRT-PCR) analysis was carried out using the Fast SYBR Green Master Mix (Roche, America) on the StepOnePlus Real-Time PCR system (Applied Biosystems, CA, USA). Primers were as follows: S100A6: forward, 5’- GGGAGGGTGACAAGCACAC-3’; reverse, 5’- AGCTTCGAGCCAATGGTGAG-3’.

Cell counting kit-8 (CCK-8) and colony formation were conducted to evaluate PC cells’ proliferative abilities. For CCK-8 assay, cells were seeded at 2,000 cells per well in 96-well plates overnight; cell growth was then monitored at various time points using the CCK-8 kit (C0038, Beyotime, China), and 2-h absorbance at 450 nm was measured on an enzyme labeling. For colony formation, about 1,000 cells were placed into each well in six-well plates. When colonies appeared, 4% methanol and crystal violet were successively used to fix and stain the cells.

The 24-well Transwell chambers precoated with Matrigel (Corning, USA) were used to perform invasion assays. About 2 × 10⁴ MIA PaCa-2 cells were seeded into the upper chambers of the Transwell in serum-free DMEM medium, while the lower chamber was supplemented with DMEM containing 10% FBS. After 24 h, the penetrated cells were fixed with 4% methanol, then stained with crystal violet, and finally photographed and counted under an inverted microscope.

A total of six GEO datasets numbered GSE38642, GSE20966, GSE25724, GSE91035, GSE55643, and GSE32688 were selected in this study; details were summarized in Table 1 . We coupled the GSE38642 and GSE91053 as discovery cohorts for the WGCNA analysis and the remaining sets as validation cohorts. TCGA database was used for copy number variation (CNV) analysis.

We identified 16 modules in GSE38642 by WGCNA, and each module is represented by a different color. Based on Spearman correlation coefficient, a heat map about module–trait relationships was drawn to assess the relationships between modules ( Figures 1A, C ). Three modules “grey,” “darkorange,” and “purple” had a highly positive association with T2DM, which were chosen as T2DM-related modules (darkorange module: r = 0.43, p = 4e−04, genes = 55; purple module: r = 0.41, p = 0.001, genes = 107; cyan module: r = 0.39, p = 0.002, genes = 290). Additionally, 14 modules were discovered in GSE91035, with the module “yellow” (r = 0.83, p = 4e−13, genes = 444) and “darkgreen” (r = 0.67, p = 2e−07, genes = 1037) being highly positively linked with PC ( Figures 1B, D ).

T2DM and PC shared 44 genes in positivity related modules ( Figure 2A ), which were considered to be extremely associated with the pathogenesis of T2DM and PC. Most of these shared genes were correlated with each other and showed a strong and significant degree of correlation in TCGA database ( Figure 2B ). Interestingly, we also found CNVs in all 44 shared genes in TCGA PC database. We observed the alterations of these genes with CNVs on the chromosome ( Figure 2C ). There were 19 CNVs that were copy number gain and 25 CNVs that were copy number loss ( Figure 2D ).

To investigate potential roles of the shared genes in T2DM and PC, we used GlueGo to run a GO analysis. The results revealed that these genes were enriched in various biological activities including the regulation of endodermal cell fate specification, regulation of dendritic cell, Interleukin (IL)-1β-related signaling pathway, ion channel activities, and regulation of timing of anagen. Regulation of endodermal cell fate specification accounted for about 40% of total GO terms ( Figures 3A, B ), indicating great importance of this pathway both in T2DM and PC. The Set4 and Set5 modules were shown to be strongly related to PC, and the genes shared with T2DM were primarily found in Set5 module. Therefore, we continued to construct the PPI network at protein level from Set5 module. Figures 4A–D displayed four clusters by MCODE method ( Figures 4A–D ). Enrichment analysis suggested that they were mostly involved in pathways of cell division, cell cycle, nuclear division, DNA replication, transcription, translation, and symbiotic process ( Figure 4E ). These results were also visualized with a GO circle diagram ( Figure S1 ). DO analysis showed that diseases enriched by the shared genes were linked with T2DM and PC ( Figure S2 ).

To validate our findings, we conducted differential genes analysis between cases and controls on validation cohorts (GSE20966, GSE25724, GSE55643, and GSE32688). The shared genes in T2DM and PC were then intersected with these differential genes. Interestingly, we obtained a hub shared gene S100A6, which may be very significant both in T2DM and PC ( Figure 5A ).

The increased expression of S100A6 was firstly validated in 176 PC tissues in comparison to 171 normal pancreatic tissues (p < 0.001) from TCGA and GTEx database ( Figure 5B ). Then, the values of S100A6 expression in 4 tumor samples were considerably higher than that in matching normal samples in TCGA cohort (p = 0.029; Figure 5C ). In addition, S100A6 was expressed significantly differently in PC tissues with different tumor grades ( Figure 5D ). To further examine S100A6 protein expression, we retrieved representative immunohistochemistry (IHC) stainings from Human protein atlas (HPA) and showed that S100A6 was remarkably overexpressed in PC samples compared to normal tissue. ( Figure 5E ). We also performed an investigation on the prognostic outcomes of S100A6 in PC patients. Kaplan–Meier (K-M) plot showed that high S100A6 expression was closely associated with poor OS (p = 0.005; Figure 5F ). Univariate Cox analysis revealed that age, grade, and S100A6 expression were markedly correlated with OS ( Table S1 ). Multivariate Cox analysis showed a significant relationship between S100A6 overexpression and low OS [hazard ratio (HR): 1.00, p = 0.0034; Figure 5G ], indicating that the shared gene S100A6 may be an independent predictor for patients with PC.

GSEA results showed that S100A6 overexpression was associated with cancer- and glycometabolism-related pathways such as pathway in cancer, p53 signaling pathway, Notch signaling pathway, and cell cycle, as well as in glycan biosynthesis, pentose phosphate pathway, sugar metabolism, galactose metabolism, and sucrose metabolism, demonstrating that S100A6 was possibly implicated in cancer growth and glucose metabolism ( Figure 6 ).

Based on the ssGSEA algorithm, we compared TIME landscape between S100A6-high and -low groups ( Figure 7A ). In addition, S100A6 was shown to be positively linked with macrophages, natural killer (NK) cells, regulatory T (Treg), T effective Memory (Tem), T helper 1 (Th1) cells, neutrophils, B cell, clusters of differentiation 8 (CD8) T cell, dentritic cell (DC), eosinophils, T gamma delta (Tgd), T central memory (Tcm), T cells, T helper cells, cytotoxic cells, mast cells, and T follicular helper (TFH) while being negatively correlated with NK CD56bright cells ( Figure 7B ).

We discovered that low-S100A6 group had higher EstimateScore, ImmuneScore, and StromalScore levels than those in high-S100A6 group (p < 0.001) ( Figure 8A ). Moreover, S100A6 was significantly negatively correlated with immune score and StromalScore but positively associated with DNAss and RNAss ( Figure 8B ). We also found that Tfh cell, CD8+ T cell, DC cell, and NK cell were increased in low-S100A6 group ( Figures 8C ). The fractions of cytolytic activity, Type II IFN, and T cell function were higher in low-S100A6 group than those in the high group ( Figures 8D ). We further discovered substantial differences in immune checkpoints such as PDCD1, CTLA4, ICOS, and LAG3 between the two groups ( Figure 8E ).

To confirm the value of S100A6 in PC, S100A6 expression levels in five PC cell lines (HPDE6-C7, BxPC-3, Mia Paca-2, PANC-1, CFPAC-1) were examined. Mia Paca-2 cell was found to have the highest S100A6 expression ( Figure 9A ), which was therefore selected for functional analysis. Then, we successfully knocked down the S100A6 expression in Mia Paca-2 cell by si-S100A6 transfection ( Figure 9B ). CCK-8 and colony formation assays showed that S100A6 knockdown cause a decreased proliferation ability of Mia Paca-2 cell ( Figures 9C, D ). Transwell assay showed that inhibiting S100A6 with indicated siRNA resulted in a significant reduction of invaded cells ( Figure 9E ).