Standardized RNA-seq and clinical data of patients with glioma, including 174 patients with GBM and 529 with lower grade glioma [World Health Organization (WHO) grade 2–3; lower grade glioma LGG] were obtained from The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/). The gene expression level of SOAT1 in 33 cancer and 31 normal tissues was obtained from the Genotype-Tissue Expression (GTEx) database (https://www.gtexportal.org/home/-index.html). The RNA-seq data was downloaded from UCSC XENA (https://xenabrowser.net/datapages/). The correlation between the SOAT1 gene expression level and the clinicopathological characteristics of glioma patients in TCGA datasets were investigated. Furthermore, the gene expression levels of SOAT1 between cancer and normal tissues were analyzed in the integrated datasets (combined TCGA with GTEx databases).

The patients were divided into high- and low-expression groups based on the median value of SOAT1 expression. Kaplan-Meier curves (K-M curves) and prognostic risk score were used to visualize the significant difference in overall survival (OS) between the high- and low-expression groups. Receiver operating characteristic (ROC) curves were used to show the prognostic value of SOAT1 in glioma. P-values and hazard ratio (HR) with 95% confidence interval (CI) were generated using log-rank tests and univariate or multivariate Cox regression analyses.

The single-sample GSEA (ssGSEA) algorithm was used to perform Spearman’s correlation analysis between SOAT1 expression and immune cell infiltration (18) using the R package “GSVA” (19). Cell type-level proportion analysis was also performed using gene expression profiling interactive analysis (GEPIA) 2021 (http://gepia2021.cancer-pku.cn/). Spearman’s correlation analysis was performed to assess the correlation between SOAT1 expression and immune infiltration through tumor immune estimation resource version 2 (TIMER2.0; http://timer.comp-genomics.org/).

To further elucidate the correlation between immune cell infiltration and SOAT1 expression in glioma, the TISIDB (an integrated repository portal for tumor-immune system interactions. http://cis.hku.hk/TISIDB/index.php) tool was used to analyze the correlation between SOAT1 expression and immune cell infiltration. The landscape of correlation between SOAT1 expression and 28 tumor-infiltrating lymphocytes across multiple cancer types was demonstrated using the “lymphocyte” module. The ESTIMATE algorithm was also used to analyze the immune and stromal scores, to evaluate the tumor immune microenvironment of SOAT1 expression in glioma. In addition, Pearson’s correlation coefficient was performed to explore the correlation between immune checkpoint gene expression and SOAT1 expression. To further investigate the association between SOAT1 expression and immune cell migration, Pearson’s correlation coefficient was performed between SOAT1 expression and chemokines/chemokine receptor genes based on TCGA cohort.

The clinicopathological characteristics of 55 patients with gliomas and 6 patients with refractory epilepsy treated by surgery in the Fifth Affiliated Hospital of Zhengzhou University (Zhengzhou Henan China) between January 2016 and December 2020 were retrospectively analyzed, and corresponding paraffin-embedded specimens were collected for immunofluorescence. Among the 55 cases of glioma, 14 cases were WHO grade 2, 21 cases were grade 3, and 20 cases were GBM. The 6 cases of refractory epilepsy were all pathologically diagnosed as gliosis.

Paraffin-embedded specimens were cut into 3-μm sections.

HE staining: Paraffin tissue sections were deparaffinized in xylene and rehydrated in degraded ethanol, respectively. After washing with dH₂O, slides were stained with hematoxylin and eosin solution in sequence followed by being washed with dH₂O. Then slides were dehydrated in degraded ethanol and immersed in xylene.

Immunofluorescence assay: Following deparaffinization with xylene and rehydration, antigen retrieval was performed by microwave treatment in ethylenediaminetetraacetic acid (pH 8.0) antigen retrieval buffer for 23 min. The endogenous peroxidase was blocked with 3% H₂O₂. Non-specific binding was blocked for 30 min using bovine serum albumin. The sections were washed in phosphate-buffered saline (PBS; pH 7.4). Diluted primary antibodies against SOAT1 (sc-69836; Santa Cruz Biotechnology, Inc.) and FoxP3 (2A11G9l; Santa Cruz Biotechnology, Inc.) were added to the tissue and incubated overnight at 4°C in a humidified chamber. Samples were then incubated with horseradish peroxidase-labelled secondary antibody for 50 min at room temperature in the dark. Sulfo-cyanine 3 tyramide signal amplification (CY3-TSA; G1223, Wuhan Servicebio Technology Co., Ltd.) and fluorescein isothiocyanate (FITC-TSA; G1222, Wuhan Servicebio Technology Co., Ltd.) were used for fluorescence signal amplification and autofluorescence quenchers (G1221; Wuhan Servicebio Technology Co., Ltd.) were used to remove autofluorescence. DAPI (G1012; Wuhan Servicebio Technology Co., Ltd.) was used for nuclear staining. Results were visualized and photographed under a fluorescence microscope (Nikon ECLIPSE C1; Nikon Corporation). Quantitative evaluation was performed by examining each section using at least three different high-power fields to calculate the percentage of positive cells in each section. Positive rate = number of positive cells/total number of cells. The median value of SOAT1 expression was 4.9%, and the percentage of positive cells in the high-expression group was ≥4.9%, while that in the low-expression group was <4.9%. The statistical summary is presented in Supplementary Table S3B . The present study was reviewed and approved by the Ethics Committee of the Fifth Affiliated Hospital of Zhengzhou University. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.

Patients were divided into the SOAT1 high- and low-expression groups based on the median SOAT1 mRNA expression value, and DEG analysis was conducted between these two groups, with adjusted P<0.05 and |log2 (FoldChange)|>2 set as the screening thresholds. To further explore the biological processes (BPs) and signaling pathways involved in SOAT1, R package “clusterProfiler” (20) was employed to perform GO and KEGG analyses based on DEG analysis. Adjusted P<0.05 was considered to indicate significant enrichment.

A PPI network was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (https://string-db.org/) with the following qualifications: “evidence”, “experiments”, “low confidence (0.150)”, “no more than 50 interactors”. Co-expression analysis was also performed using “R” software with the screening thresholds: |Pearson’s correlation coefficient|>0.6, P-value<0.05. A Venn diagram was used to perform an intersection analysis of the PPI network and co-expression analysis. To further investigate the correlation between SOAT1 and intersection genes, Pearson’s correlation coefficient was adopted.

All statistical analysis was performed using R software (version 3.6.3). According to whether the data was normally distributed and whether the variance was uniform, Wilcoxon rank sum test, Welch t-test, or t-test were used to compare the difference in SOAT1 expression between two different groups, respectively. One-way ANOVA, Welch one-way ANOVA or Kruskal-Wallis test were performed for multi-group comparisons. Kaplan-Meier curves and log-rank tests were performed to compare the difference in survival between the high- and low-expression groups. ROC curves, prognostic risk score, and univariate or multivariate Cox proportional hazards regression analysis were performed to evaluate the prognostic value of SOAT1 in glioma. P<0.05 was considered to indicate a statistically significant difference.

A pan-cancer analysis was performed to compare the mRNA level of SOAT1 between tumor and corresponding normal tissues in the integrated datasets (combined TCGA with GTEx databases), which revealed that SOAT1 was upregulated in 22/33 tumor types and downregulated in 3 tumor types ( Figure 1A ). In addition, the mRNA level of SOAT1 was significantly upregulated in both LGG and GBM compared to normal brain tissues ( Figures 1B, C ). Furthermore, SOAT1 expression was elevated as the tumor grade increased ( Figure 1D ). The SOAT1 expression of isocitrate dehydrogenase (IDH) wild type was higher than that of IDH mutant type in WHO grade 3 and overall gliomas ( Figure 1F ; Supplementary Table S1 ). SOAT1 expression in 1p19q non-co-deleted gliomas was higher than that of 1p19q-co-deleted gliomas ( Figures 1E , H ; Supplementary Table S1 ). Since IDH wild type and 1p19q non-co-deletion are malignant biomarkers of glioma (21), the high SOAT1 expression may predict the advanced malignancy of glioma. Nevertheless, no difference in SOAT1 expression was identified between different histological types and IDH types in grade 2 gliomas ( Figures 1G, I, J ). In combination, SOAT1 was widely overexpressed in a variety of tumors, and the high SOAT1 expression was found to be correlated with malignant clinicopathological characteristics in glioma.

Using K-M curves and log-rank tests, the high SOAT1 expression was found to be associated with poor prognosis in overall and grade 3 gliomas, but not in grade 2 gliomas and GBM ( Figures 2A–D ). The prognostic risk score demonstrated that the fatality rate in the SOAT1 high-expression group was significantly higher than that in the low-expression group. The risk score was positively correlated with the SOAT1 expression level ( Figure 2E ). To further observe the predictive effect of SOAT1 expression on the OS of glioma patients, ROC curves were employed. It was found that SOAT1 had a marked predictive ability for the 1-, 3-, and 5-year OS of overall and grade 3 gliomas. The survival area under the curve (AUC) of SOAT1 expression was 0.766, 0.744, and 0.708 for the 1-, 3-, and 5-year OS, respectively, in overall glioma. The AUC of SOAT1 expression was 0.755, 0.645, and 0.603 for the 1-, 3-, and 5-year OS, respectively, in WHO grade 3 gliomas ( Figures 2F, G ).

In addition, univariate Cox regression analysis indicated that SOAT1 was a poor prognostic factor for glioma patients; the HR (95% CI) was 2.863 (2.358-3.477). However, multivariate Cox regression analysis lost statistical significance ( Supplementary Table S2 ). Furthermore, according to WHO 2021 classification of brain tumors, we classified gliomas into astrocytoma IDH-mutant grade 2, astrocytoma IDH-mutant grade 3, oligodendroglioma grade 2, oligodendroglioma grade 3, astrocytoma IDH wild type, and GBM in TCGA cohort. We divided patients into high/low groups by the median value of SOAT1 expression and performed the survival analysis in each type of glioma. Unfortunately, no statistical significance of survival analysis was identified between SOAT1 -high and -low groups in all types of gliomas ( Supplementary Figure S5 ). Only the astrocytoma IDH mutant grade 3 showed a little tendency toward poor prognosis in the high-expression group (HR=1.69, 0.82-3.51, P=0.16. Supplementary Figure 5A ). Therefore, high expression of SOAT1 was correlated with poor prognosis in glioma, but SOAT1 was not an independent prognostic factor in glioma.

Immunotherapy, including therapeutic vaccines and engineered T-cells based on tumor-immune cell interactions and blockers, has become an important strategy in the field of cancer research, including glioma (22–26). The relationship between SOAT1 expression and immune cell infiltration in the TME was determined using the ssGSEA algorithm based on the TCGA cohort. The expression of SOAT1 was found to be positively correlated with macrophages, neutrophils, T effector memory (Tem), T helper 17 cells (Th17), and activated dendritic cells (aDC), and negatively correlated with plasmacytoid DC (pDC), natural killer (NK), and CD56 bright cells in the GBM ( Figure 3A ). SOAT1 expression was positively correlated with T helper cells, macrophages, aDC, eosinophils, and T helper 2 cells (Th2), and negatively correlated with NK CD56 bright and pDC cells in the LGG ( Figure 3A ). Next, the relationship between the proportion of immune cell infiltration and the expression level of SOAT1 was demonstrated by searching the GEPIA2021 database ( Figure 3B ).

To further identify and quantify the immune and matrix components in GBM and LGG, immune, stromal, and ESTIMATE scores were employed, and the results indicated that SOAT1 expression was positively correlated with all three scores in GBM and LGG ( Figure 3C ). In addition, correlation analysis results for SOAT1 expression and several tumor infiltrating lymphocytes (TILs) from the TIMER2.0 database are displayed in Figures 3D, E . The expression level of SOAT1 was negatively correlated with purity of tumor in GBM and LGG (r=-0.27, P=1.36e-03, r=-0.03, P=-5.10e-01) ( Figures 3D, E ). Of note, SOAT1 expression was positively correlated with multiple TILs in GBM, including macrophages (r=0.522, P=5.85e-11), neutrophils (r=0.436, P=9.76e-08), and myeloid dendritic cells (r=0.641, P=3.15e-17). In addition, SOAT1 expression was positively correlated with multiple TILs in LGG, including CD4+ T cells (r=0.409, P=9.44e-21), myeloid dendritic cells (r=0.569, P=2.23e-42), and neutrophils (r=0.742, P=7.63e-85).

To further validate the correlation between SOAT1 expression and immune infiltration, the correlation between SOAT1 expression and 28 TILs from the TISIDB database was examined. The pan-cancer analysis of the correlation between SOAT1 expression and 28 TILs is displayed via a heatmap ( Figure 4A ). Specifically, SOAT1 expression was positively correlated with multiple TILs in GBM and LGG ( Figures 4B–T ; Supplementary Figures S1 and S2 ). The above results indicated that the expression level of SOAT1 was positively correlated with various tumor infiltrating immune cells in the TME.

Treg is a subgroup of CD4+ cells. In previous bioinformatics analysis, SOAT1 expression was positively correlated with multiple CD4+ T cells and Tregs based on different dataset ( Figures 3A, E , Figures 4B, J, K, S ; Supplementary Figures S1Q, S , and Supplementary Figures S2J, R ). Forkhead Box P3 (FoxP3) is a protein-coding gene that is crucial for the development and inhibitory function of Treg (27). FoxP3 appears to be distinctive for Tregs as it is expressed in CD4+ CD25+ T cells and CD4+ CD25+ CD8- thymocytes whereas it is not found in other thymic cells, T cells, B cells, natural killer cells, or natural killer T cells (28, 29). Foxp3 also was found to control the expression of gene programs, which define and maintain Treg cell identity and function (30, 31). Therefore, FoxP3 is reliable and a constant marker that has served to isolate and characterize Tregs. Furthermore, Foxp3+ Tregs are impacted by different environmental conditions and metabolic differences associated with diverse transcriptional patterns (32). Foxp3+ Tregs have been extensively studied in gliomas and frequently infiltrate high-grade malignant gliomas (33). Since FoxP3+ Tregs play a crucial role in glioma-mediated immunosuppression (34), we wonder whether the SOAT1 expression was associated with immunosuppression in glioma.

To further validate this correlation, immunofluorescence was performed by double-labeling SOAT1 and FoxP3 on paraffin-embedded specimens from 55 cases of glioma and 6 of gliosis. Pathologists have made diagnosis based on histology and molecular characteristics. Histologically, among the 55 gliomas, 20 were GBM, 21 were grade 3 gliomas, and 14 were grade 2 gliomas. Of the 21 grade 3 gliomas, 13 were oligodendroglioma grade 3 and 8 were astrocytoma grade 3; among the 14 grade 2 gliomas, 5 were oligodendrogliomas grade 2, and 9 were astrocytoma grade 2. Consistent with our previous study, SOAT1 expression was increased as the tumor grade increased ( Figure 5A ; Supplementary Table S3A ), while there was no difference in clinicopathologic features such as age, gender, epilepsy, histological type, preoperative Karnofsky performance status, and tumor diameter between the high- and low-expression groups. The difference of SOAT1 expression between the astrocytoma grade 2–3 and the oligodendroglioma grade 2–3 did not reach statistical significance (P=0.489) probably due to the small sample size. FoxP3+ Treg infiltration was observed in 17 cases of glioma, mostly in high-grade gliomas (9 GBM, 7 WHO grade 3, and 1 WHO grade 2). Although the difference in FoxP3+ Treg infiltration between SOAT1 high- and low-expression groups did not reach statistical significance (χ2 = 0.24, P=0.622), the positive rate of FoxP3 in the SOAT1 high-expression group (35.7%) was higher than that in the SOAT1 low-expression group (25.9%) ( Supplementary Table S3A ). In addition, it was found that some fields of SOAT1 high-expression group co-localized with FoxP3+ Treg infiltration ( Figure 5B ). The above results suggest that a positive trend of the correlation between SOAT1 expression and Treg infiltration was observed in this cohort and some of the Treg cells or tumor cells in their immediate vicinity may be expressing SOAT1.

Immune checkpoints have been proven to be involved in tumor immunosurveillance escape (35). Immunotherapy based on immune checkpoint blockade is increasingly considered as the most promising therapy for GBM in addition to operative treatment (36). Pearson’s correlation coefficient was performed between the SOAT1 expression and 35 immune checkpoint genes in LGG and GBM based on the TCGA dataset. The results demonstrated that SOAT1 expression was positively correlated with various checkpoint genes in GBM and LGG ( Figures 6A, B ). In particular, SOAT1 expression had a significant positive correlation with PD-L1 (CD274, Aliases for PD-L1 gene) (r=0.52 P=4.01e-38 in LGG, r=0.39 P=2.06e-07 in GBM), monoclonal antibodies against which have been used in clinical trials of GBM (37). These results indicated that SOAT1 may be involved in tumor immune escape.

Next, the relationship between SOAT1 expression and multiple chemokine/chemokine receptor gene levels in GBM and LGG was comprehensively analyzed. As shown in Figures 6C–F , SOAT1 expression was positively correlated with multiple chemokine gene levels such as CCL13 (r=0.35, P=2.68e-06), CCL20 (r=0.37, P=6.30e-07), and CCL18 (r=0.35, P=2.82e-06) in GBM, and CCL5 (r=0.46, P=1.70e-29), CCL2 (r=0.39, P=1.85e-20), and CCL20 (r=0.35, P=3.81e-17) in LGG. SOAT1 expression was also positively correlated with multiple chemokine receptor gene levels, including those of CCR5 (r=0.58, P=9.09e-17), CCR1 (r=0.55, P=5.17e-15), and CCR2 (r=0.47, P=1.26e-10) in GBM, and CCR5 (r=0.665, P=2.97e-67), CCR1 (r=0.63P=2.95e-60), and CCR2 (r=0.53, P=1.71e-40) in LGG. These results indicated that SOAT1 overexpression may increase the secretion of chemokines to attract immune infiltration into the TME.

To further understand the role of SOAT1 in glioma pathogenesis, DEG analysis was performed based on the TCGA dataset, and the |logFC| and adjusted P-value were visualized using volcano plots ( Figure 7A ). The results showed that 662 genes were upregulated, and 289 genes were downregulated under the screening threshold (adjusted P<0.05 and |log2 (FoldChange)| >2). GO and KEGG analyses based on DEGs analysis demonstrated that a variety of pathways and BPs were significantly enriched, including antigen binding, immunoglobulin complex, humoral immune response (GO analysis; Supplementary Figure S3A ), neuroactive ligand-receptor interaction, PI3K/AKT signaling pathway, viral protein interaction with cytokine and cytokine receptor, chemokine signaling pathway, and cAMP signaling pathway (KEGG analysis; Figure 7B ). To further visualize the genes involved in the significantly enriched pathway, a chordal graph was used ( Figure 7C ; Supplementary Figures S3B and S4A ). These results demonstrated that SOAT1 was involved in a variety of BPs and signaling pathways in the pathogenesis of glioma.

A PPI network was constructed using the STRING website and the 50 identified SOAT1-binding proteins ( Figure 8A ). Subsequently, a co-expression analysis of SOAT1 in glioma was performed using Pearson’s correlation coefficient, and the results showed that 1287 genes were positively correlated, and 24 genes were negatively correlated with SOAT1 expression in glioma under the screening threshold (adjusted P<0.05 and |Pearson’s correlation coefficient|>0.6). In addition, the top 5 positively correlated genes were MOB kinase activator 1A (MOB1A), CKLF like MARVEL transmembrane domain containing 6 (CMTM6), major facilitator superfamily domain containing 1 (MFSD1), LHFPL tetraspan subfamily member 2 (LHFPL2), and LIM and senescent cell antigen-like-containing (LIMS1), and the top 5 negatively correlated genes were mitochondrially encoded 12S RRNA (MT-RNR1), mitochondrially encoded 16S RRNA (MT-RNR2), mitochondrially encoded cytochrome C oxidase II (MT-CO2), Rho GDP dissociation inhibitor gamma (ARHGDIG), and RUN domain containing 3A (RUNDC3A) ( Figure 8B ).

The Venn diagram showed that 4 genes were shared between the PPI network and the co-expression analysis: Calnexin (CANX), embryonic ectoderm development (EED), katanin catalytic subunit A1 (KATNA1), and phospholipase A2 group XV (PLA2G15) ( Figure 8C ). The correlations between SOAT1 expression in glioma were demonstrated in Figures 8D–G . CANX is a calcium-binding, endoplasmic reticulum-associated protein that interacts transiently with newly synthesized N-linked glycoproteins, facilitating protein folding and assembly. The interaction between SOAT1 and CANX was detected by co-sedimentation via a density gradient assay (38). Embryonic ectoderm development encodes a member of the polycomb-group (PcG) family. PcG family members form multimeric protein complexes, which are involved in maintaining the transcriptional repressive state of genes over successive cell generations, and the interaction between SOAT1 and CANX was detected by affinity chromatography technology assay (39). KATNA1 is the catalytic subunit of a complex which severs microtubules in an ATP-dependent manner. The interaction between SOAT1 and KATNA1 was detected using an anti-tag coimmunoprecipitation assay (https://www.ebi.ac.uk/). PLA2G15 encodes the protein that hydrolyzes lysophosphatidylcholine to glycerophosphorylcholine and a free fatty acid. The interaction between SOAT1 and PLA2G15 was detected by genetic interference assay (40). In combination, these genes may play a significant role in the regulatory mechanism of SOAT1 in gliomas.