The mRNA expression, copy number alteration, and clinical information of more than 30 cancer types were obtained from The Cancer Genome Atlas (TCGA). The mRNA expression data from normal tissue sites across nearly 1,000 people were obtained from the Genotype-Tissue Expression (GTEx) to supplement the normal tissue mRNA sequencing data lacking in TCGA. Siglec-9 mRNA expression data of tumor cell lines were obtained from the Broad Institute Cancer Cell Line Encyclopedia (CCLE). The mRNA expression data and clinical information of glioma were obtained from the Chinese Glioma Genome Atlas (CGGA) to verify the result in TCGA-LGG. We used the “mRNAseq_693” dataset (Wang et al., 2015; Liu et al., 2018). WHO grade II and WHO grade III cases were defined as LGG to explore the role of Siglec-9 in LGG (Liu et al., 2018).

GEPIA2 (http://gepia2.cancer-pku.cn/#index) (Tang et al., 2019) was used to analyze the correlation between tumor stages and Siglec-9 mRNA expression using “major stage” and “log2 (TPM + 1) for log-scale”. The association between Siglec-9 mRNA expression and overall survival (OS), disease-specific survival (DSS), disease-free interval (DFI), and progression-free interval (PFI) depended on TCGA databases and was analyzed on the DriverDBv3 (http://driverdb.tms.cmu.edu.tw/) (Liu et al., 2020) with mean of Siglec-9 expression as cutoff value. The online database UALCAN network (http://ualcan.path.uab.edu/) (Chandrashekar et al., 2017) was also used to verify the OS analysis with default settings. Siglec-9 expression in diverse molecular subtypes was analyzed on The Integrated Repository Portal for Tumor-Immune System Interactions (TISIDB, http://cis.hku.hk/TISIDB/index.php) (Ru et al., 2019).

TIMER (https://cistrome.shinyapps.io/timer/) (Li et al., 2017) and TIMER2.0 (http://timer.comp-genomics.org/) (Li et al., 2020) were used to assess the correlation between the Siglec-9 expression and immune cell infiltration for TCGA tumors and CGGA-LGG dataset with TIMER score and CIBERSORT score (Newman et al., 2015) using immune infiltrates query function and estimation function. The results of estimation were visualized by R package “psych”. TISIDB, a user-friendly web instrument to explore comprehensive investigation of tumor-associated immunity, was used to analyze the association between Siglec-9 mRNA expression and tumor immune subtypes, and specific types of immune cell infiltration. We estimated Siglec-9 expression in diverse molecular subtypes, and immune subtypes involving C1 (wound healing), C2 (IFN-γ dominant), C3 (inflammatory), C4 (lymphocyte deplete), C5 (immunologically quiet), and C6 (TGF-β dominant) subtypes. Specific types of immune cell infiltration of 28 TIL types were inferred by using gene set variation analysis (GSVA) based on the Siglec-9 expression profile on TISIDB.

GSEA based on WEB-based Gene SeT AnaLysis Toolkit (WebGestalt, http://www.webgestalt.org/) (Liao et al., 2019) was performed with the TCGA-LGG dataset on LinkedOmics (Vasaikar et al., 2018) (http://www.linkedomics.org/login.php), a publicly available portal that includes multi-omics data from all 32 TCGA Cancer types. CGGA mRNA expression data were downloaded and LGG cases were selected to be analyzed with GSEA on WebGestalt. The parameters for the enrichment analysis were as follows: minimum number of IDs in the category: three; maximum number of IDs in the category: 2000; significance level: Top 25; number of permutation: 1,000. The top ten positive related categories are shown in the main figure.

The expression of Siglec-9 in different tissue was used by Kruskal–Wallis test, and between tumor tissues and normal tissues were estimated by t-test. Additionally, the expression of Siglec-9 in different grades of glioma was utilized by t-test and ANOVA test. Kaplan–Meier curves were visualized to compare the survival patients stratified based on different levels of expression of Siglec-9. The relationship between Siglec-9 and TMB (tumor mutation burden), MSI (microsatellite instability), MMR gene mutation, immune checkpoints, DNMT, immune score, stromal score, ESTIMATE score, and immune cells was evaluated by Pearson and Spearman correlation analyses. p < 0.05 was recorded as statistically significant for all analyses unless otherwise specified.

The expression of Siglec-9 mRNA in normal tissue was firstly analyzed using the GTEx database. It showed that the level of Siglec-9 expressions was highest in tissues of blood, spleen, and lung and the lowest in muscle and bone marrow among all types of tissues assessed (Kruskal–Wallis test, p = 0) (Supplementary Figure S1).

Then, the expression of Siglec-9 in tumor and normal samples from various tissue types was assessed using TCGA database. Differential expression analysis was performed with normal tissues used as control. Siglec-9 was upregulated in tumor tissues from cancers including breast invasive carcinoma (BRCA), glioblastoma multiforme (GBM), head and neck squamous cell carcinoma (HNSC), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), brain lower-grade glioma (LGG), stomach adenocarcinoma (STAD), thyroid carcinoma (THCA), and uterine corpus endometrial carcinoma (UCEC) (Figure 1A). Siglec-9 was downregulated in cancers including colon adenocarcinoma (COAD), liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), and pancreatic adenocarcinoma (PAAD) (Figure 1A).

For the tumors lacking normal sample controls in TCGA, the corresponding normal samples from GTEx databases were included for further differential expression analysis. The result showed that Siglec-9 was upregulated in BRCA, cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), esophageal carcinoma (ESCA), acute myeloid leukemia (LAML), LGG, ovarian serous cystadenocarcinoma (OV), pancreatic adenocarcinoma (PAAD), skin cutaneous melanoma (SKCM), STAD, testicular germ cell tumors (TGCT), THCA, UCEC, and uterine carcinosarcoma (UCS), and was downregulated in adrenocortical carcinoma (ACC) (Figure 1B).

In addition, the expression pattern of Siglec-9 in the tumors of different clinal stages and molecular subtypes in pan-cancer was also investigated. It showed that Siglec-9 expression was significantly different among clinal stages in bladder urothelial carcinoma (BLCA), ESCA, KICH, THCA, and SKCM (Figure 2A). The different expression pattern was also observed among molecular subtypes in ACC, BRCA, COAD, HNSC, KIRP, LGG, LIHC, LUSC, OV, PCPG, PRAD, and STAD (Figure 2B). No different expression was found among tumor stages or molecular subtypes in other tumors such as ACC, BRCA, and CESC (Supplementary Figure S2A ).

The prognostic value of Siglec-9 expression was assessed by investigating the relationship between Siglec-9 expression and clinical characteristics including OS, DSS, DFI, and PFI using DriverDBv3. Higher Siglec-9 expression level was found to be associated with shorter OS and 5-year OS (HR > 1) in LUSC, THYM, and LGG, but longer OS and 5-year OS (HR < 1) in ACC (Table 1; Figures 3A,B; Supplementary Figure S3A). DSS analysis was performed to further verify these correlations and found similar results (Figures 3C,D; Supplementary Table S1 ). High Siglec-9 expression correlated with the poor outcome of PFI in GBM, PRAD, and LGG, while it correlated with a good outcome of PFI in ACC (Supplementary Table S2; Figure 3E; Supplementary Figure S3B). DFI analysis found that increased Siglec-9 expression was associated with shorter OS and 5-year survival in ESCA only (Supplementary Table S3; Supplementary Figure S3C ). OS analysis was also performed using UALCAN for verification. It showed that Siglec-9 was a poor prognostic factor in LGG, SKCM, UVM, and a good one in READ (Supplementary Table S4). A significant correlation of Siglec-9 with the prognosis in LGG was observed from both DriverDBv3 and UALCAN analysis (Table 1; Supplementary Table S4).

Previous research has identified six molecular immune subtypes containing wound healing (C1), IFN-γ dominant (C2), inflammatory (C3), lymphocyte depletes (C4), immunologically quiet (C5), and TGF-β dominant (C6) subtypes, which were related to tumor molecular characterization and prognosis of patients (Thorsson et al., 2018). We analyzed the expression pattern of Siglec-9 in six immune subtypes in TCGA pan-cancer. The Siglec-9 was differentially expressed in different immune subtypes in ACC, BLCA, BRCA, CESC, COAD, KICH, KIRC, LGG, LUSC, LIHC, LUAD, Mesothelioma (MESO), OV, PAAD, PCPG, PRAD, READ, SARC, SKCM, STAD, TGCT, THCA, UCEC, UCS, and UVM (Figure 4; Supplementary Figure S4A). No difference among different immune subtypes was observed in other tumors including CHOL, ESCA, GBM, HNSC, and KIRP (Supplementary Figure S4B).

Then, the association between Siglec-9 expression and immune cell infiltration in tumor tissues was analyzed using TIMER. A positive correlation between Siglec-9 expression and the infiltration of B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, macrophages, and dendritic cells was observed in most TCGA tumors including LGG, ACC, and LUSC (Figures 5A–C; Supplementary Figure S5). In THYM, a negative correlation was observed in B-cell, CD4⁺ T-cell, and CD8⁺ T-cell infiltration and a positive correlation was observed in neutrophil infiltration (Figure 5D ). Besides, Siglec-9 expression was negatively correlated with CD8⁺ T cells in GBM, with B cells in STAD and UVM, and with macrophages in DLBC (Supplementary Figure S5). ESTIMATE analysis was applied to further investigate the correlation between Siglec-9 expression and the tumor microenvironment. A strong positive correlation was observed in immune scores (Supplementary Figure S6), stromal scores (Supplementary Figure S7), and ESTIMATE scores (Supplementary Figure S8) in all 33 tumor types in TCGA including ACC, LGG, LUSC, and THYM (Figures 6A–C). To explore the association between Siglec-9 expression and specific types of immune cell infiltration, the relative abundance of 28 TIL types was inferred by using gene set variation analysis (GSVA) based on the Siglec-9 expression profile on TISIDB. The result showed that Siglec-9 expression was strongly correlated with follicular helper T cells (Tfh), regulatory T cells (Treg), myeloid-derived suppressor cells (MDSC), and macrophages across 33 TCGA tumors (Figure 6D). Moreover, the correlation between Siglec-9 expression and the expression of 47 immune checkpoint genes was also analyzed. Siglec-9 was positively correlated with immune checkpoints including LAIR1, HAVCR2, CD80, PDCD1 (PD1), PDCD1LG2 (PD-L2), VSIR, and CD86 across most TCGA tumors (Figure 6E). Siglec-9 had a low correlation with immune checkpoints in tumors including ACC, CHOL, DLBC, LAML, MESO, THYM, and UCS. In particular, LAIR1, HAVCR2, LGALS9, and CD86 expression were strongly positively correlated with Siglec-9 expression, and only CD200 expression was negatively correlated with Siglec-9 expression in LGG. In LUSC, Siglec-9 expression was strongly positively correlated with LAIR1, HAVCR2, and CD86, and negatively correlated with VTCN1 and TNFRSF18.

TMB, MSI, DNA MMR, and DNA methylation were considered immune-related and prognosis-related factors in tumors (Duffy and Crown, 2019). Therefore, we analyzed the correlation between the Siglec-9 expression and TMB, MSI, DNA MMR genes, and DNMT. Siglec-9 expression was found to be positively correlated with TMB in CESC, COAD, LGG, OV, PRAD, SARC, and THYM and negatively correlated with TMB in GBM and THCA (Figure 7A). Siglec-9 expression was negatively related to MSI in HNSC, LGG, LUAD, LUSC, PAAD, SKCM, STAD, and TGCT and positively related to MSI in COAD (Figure 7B). Five DNA MMR genes (MLH1, MSH2, MSH6, PMS2, and EpCAM) were included for analyzing the association between MMRs and Siglec-9. MLH1 was positively correlated with Siglec-9 in BLCA, BRCA, ESCA, HNSC, KIRC, LIHC, LUAD, LUSC, PAAD, PRAD, STAD, UCEC, and negatively correlated in GBM and SARC (Figure 7C). MSH2 had a positive correlation with Siglec-9 in KIRC, LIHC, and PRAD, and had a negative correlation in (Figure 7C). MSH6 showed a positive association with Siglec-9 in HNSC, KIRC, LGG, LIHC, PAAD, PRAD, and STAD, and a negative association in GBM, LAML, LUSC, SARC, TGCT, and UCS (Figure 7C). PMS2 was positively correlated with Siglec-9 in KIRC, LIHC, PAAD, PRAD, and STAD, and negatively correlated in GBM, LUSC, SARC, SKCM, and THCA (Figure 7C). Siglec-9 was negatively correlated to EpCAM in most tumors including BLCA, BRCA, CESC, CHOL, COAD, KIRC, KIRP, LAML, LGG, LUAD, LUSC, OV, PAAD, PRAD, STAD, TGCT, and THCA (Figure 7C). Collectively, Siglec-9 had a significant association with MMRs in KIRC, PAAD, and PRAD (Figure 7C).

DNA methyltransferase1 (DNMT1), DNA methyltransferase2 (DNMT2), DNA methyltransferase3 (DNMT3), and DNA methyltransferase3B (DNMT3B) were included for the analysis of the correlation between Siglec-9 expression and DNA methyltransferases. It showed that Siglec-9 expression was correlated with DNA methylation in various cancers, especially in BLCA, GBM, HNSC, KICH, KIRC, LGG, PAAD, PRAD, TGCT, and THCA (Figure 7D). In particular, Siglec-9 was positively correlated with DNMT1 (Pearson’s rho = 0.11, p = 0.029), DNMT2 (Pearson’s rho = 0.2, p = 8.5 e-6), and DNMT3 (Pearson’s rho = 0.21, p = 5.4 e-6). There was no significant correlation or less correlation between Siglec-9 expression and DNA methyltransferases in ACC, LUSC, and THYM.

The analysis above showed that Siglec9 expression was correlated with the prognosis of patients, especially of ACC and LGG patients (Figure 3). Due to the larger sample size, similar survival pattern to other tumors, and closer association with the tumor immune microenvironment, TMB, MSI, MMRs, and DNMTs, we selected LGG for further analysis. Data from the Chinese Glioma Genome Atlas (CGGA) was employed for further verification. Expression of Siglec-9 in different grades of glioma (WHO II, WHO III, and WHO IV) was analyzed. The results showed that Siglec-9 expression was higher in WHO IV (GBM) than that in WHO II and WHO III (LGG) (Figure 8A). OS analysis showed that high expression of Siglec-9 was correlated with poor prognosis in primary WHO grade II and recurrent WHO grade III glioma (Figure 8B). There was no significant correlation between OS and Siglec-9 expression in recurrent WHO grade II (Supplementary Figure S9A) and primary WHO grade III glioma (Supplementary Figure S9B). In addition, Siglec-9 expression in IDH mutation and 1p/19q co-deletion status subtypes was analyzed, because IDH mutation status and 1p/19q co-deletion status were beneficial to diagnosis, prognosis prediction, and treatment of glioma (Aoki et al., 2018). No difference between IDH mutant and wild type in WHO grade II was observed, and higher Siglec-9 expression was observed in wild type compared to IDH mutant in WHO grade III (Figure 8C). Besides, Siglec-9 expression was significantly higher in the none 1p/19q co-deletion than in the 1p/19q co-deletion subtype in both WHO grade II and WHO grade III (Figure 8D), which coincided with the former result of low expression of Siglec-9 in Codel subtype in the TCGA-LGG (Figure 2B). Subsequently, Siglec-9 expression in subtypes that combined IDH mutation with 1p/19q was also analyzed in LGG and GBM. The result showed that Siglec-9 was significantly differentially expressed among these subtypes (Figure 8E). Then, TIMER was used to assess the correlation between Siglec-9 expression and immune cell infiltration in CGGA-LGG. Similar to TCGA-LGG, Siglec-9 expression was positively correlated with B cells, CD4⁺ T cells, neutrophils, macrophages, and myeloid dendritic cells except for CD8⁺ T cells according to TIMER immune score in CGGA-LGG (Figure 9A). A strong correlation was observed between Siglec-9 expression and immune score, stromal score, and ESTIMATE score (Figure 9B). Analysis of correlation between Siglec-9 expression and specific immune cells indicated that high Siglec-9 expression was accompanied by high infiltration level of M1/2 macrophages and neutrophils, and low infiltration level of activated NK cells, Tfh, and naïve CD4⁺ T cells (Figure 9C).

To clarify the specific mechanism by which siglec9 affected LGG survival, GSEA analysis was used to identify key biological functions of Siglec-9, between the low and high Siglec-9 subgroups in the TCGA-LGG and CGGA-LGG (WHO grade II and III). Results indicated that various pathways involved in adaptive immune response, cytokine production, inflammatory response, antigen processing, and presentation were enriched in the high Siglec-9 expression subgroup in Gene Ontology (GO) enrichment analysis (Figures 10A,B; Supplementary Figures S10A–F) and KEGG analysis (Figures 10C,D; Supplementary Figures S10G,H). Immune process-related cellular components, such as MHC protein complex, immunological synapse, and secretory granule membrane, were significantly enriched in GO enrichment analysis (Supplementary Figures S10E,F). Interestingly, neural cell function-related pathways were significantly enriched in low Siglec-9 expression in all analyses mentioned above.