Abstract Glutaminase (GLS), a crucial gene regulating glutaminolysis, has received much attention as it was found to regulate tumor metabolism and copper-induced cell death. However, its biological roles and mechanisms in human cancers remain obscure. Consequently, the integrated pan-cancer analyses and biological experiments were conducted to elucidate its oncological functions. We found GLS was differentially expressed in human cancers and upregulated GLS predicted poor survival, clinicopathological progression, and tumor heterogeneity. Single-cell analysis found GLS was closely related to various biological functions and pathways. Spatial transcriptomic analysis found GLS expression was mainly derived from tumor cells, which implies tumor cells may have a stronger ability to utilize glutamine than antitumor immune cells in the tumor microenvironment (TME). Meanwhile, we noticed GLS expression was strongly related to the infiltration of various immune cells and stromal cells, the expression of immunomodulatory genes, the activity of some conventional antitumor agents, and the therapeutic response of immunotherapy. Moreover, enrichment analyses suggested GLS was related to various metabolic reprogramming, innate and adaptive immunity suppression, and extracellular matrix remodeling. Finally, we observed GLS was highly expressed in our gastric cancer (GC) cohort. As an independent risk factor for GC prognosis, high-GLS was closely related to pathological progression. Inhibiting GLS expression in GC cells effectively prevented proliferation, migration, and invasion and triggered apoptosis. In conclusion, GLS is an underlying biomarker for oncological progression, prognosis, TME, antitumor drug sensitivity, and immunotherapy response. Targeting GLS can facilitate the implementation of individualized and combined treatment strategies. Supplementary Information The online version contains supplementary material available at 10.1038/s41598-024-84916-w. Keywords: GLS, Pan-cancer analysis, Gastric cancer, Biomarker, Tumor microenvironment, Immunotherapy Subject terms: Tumour biomarkers, Bioinformatics, Immunotherapy, Cancer microenvironment, Gastrointestinal cancer, Tumour biomarkers, Tumour immunology, Gene expression Introduction Globally, cancer has been recognized as a leading cause of death, as its incidence and mortality are rising quickly^[38]1. The aberrant alternations of cancer-associated genes drive the complex process of tumorigenesis, progression, and evolution, resulting in significant heterogeneity and poor prognosis, which remains considerable challenges for tumor treatment^[39]2–[40]4. Immunotherapy has made remarkable progress and has become a key therapeutic strategy for multiple tumors in recent years, whereas only subsets of tumor patients benefit from it^[41]5. Consequently, it is extremely indispensable to seek potential genetic biomarkers for predicting which patients may be suitable for immunotherapy. In addition to glucose metabolism, uncontrolled glutamine metabolism is another crucial prerequisite for tumor biosynthesis and energy supply^[42]6. Specifically, glutamine is a pivotal carbon pool and nitrogen pool for malignant cells, and its catabolism provides precursors for the biosynthesis of nucleotides and non-essential amino acids, thereby facilitating the proliferation and division of tumor cells^[43]7–[44]9. Besides, glutamine participates in the generation of glutathione and nicotinamide adenine dinucleotide phosphate to maintain cellular redox balance, thereby creating a favorable microenvironment for cellular malignant transformation and survival^[45]10,[46]11. Moreover, due to the decreased flux of oxidative phosphorylation in tumor cells resulted from the Warburg effect, glutamine endlessly supplies tricarboxylic acid (TCA) cycle intermediates α-ketoglutarate to replenish the energy expenditure needed for sustained proliferation^[47]12,[48]13. Briefly, tumor cells can reprogram aforementioned glutamine metabolism pathways to support their energy requirement and biosynthesis. Glutaminase (GLS), a crucial regulator of glutamine metabolism, is responsible for hydrolyzing glutamine to glutamate and ammonia^[49]14. Physiologically, GLS expression is primarily concentrated in the kidney and brain tissues, involved in cellular energy metabolism, cerebral neurotransmitter synthesis, and renal acid-base maintenance^[50]15,[51]16. However, GLS expression is often dysregulated in malignant cells as metabolic reprogramming occurs. Its dysregulation continuously drives glutaminolysis for the complement of TCA cycle and biosynthetic consumption, assisting malignant tumors to spread wildly. Therefore, targeting GLS-mediated glutamine metabolism to deprive the glutamine utilization of tumor cells is promising as a therapeutic strategy for cancer, which has received extensive attention. Remarkably, a recent study has reported GLS, as a cuproptosis-related gene, negatively regulates copper-triggered cell death^[52]17. Since this novel cell death pattern may be involved in the tumorigenesis and progression, it once again arouses academic interest in GLS research. Although GLS is recognized to drive tumorigenesis, the detailed roles of GLS in prognostic outcome, tumor microenvironment (TME), immunoregulation, and immunotherapy still need to be systematically investigated based on multi-omics analyses. In this research, a systematic pan-cancer analysis was performed to mainly elucidate the association of GLS with tumor progression, prognosis, single-cell function, genetic alteration, antitumor drug sensitivity, immune microenvironment, and immunotherapy response, aiming to clarify its biological functions in human tumors and provide a reference for decision-making of individualized and combined antitumor therapies. In addition, we validated its expressed differences in gastric cancer (GC) as well as its potential roles in tumor progression and prognosis using a clinical GC cohort and in vitro biological experiments. Briefly, we found GLS is an underlying prognostic, therapeutic, and immunological biomarker in human cancers. Materials and methods Data collection Transcriptomic datasets of 33 cancer types and their matching clinical data were derived from The Cancer Genome Atlas (TCGA) ([53]https://portal.gdc.cancer.gov), Genotype-Tissue Expression (GTEx) ([54]https://gtexportal.org), and Therapeutically Applicable Research to Generate Effective Treatments (TARGET) ([55]https://ocg.cancer.gov/programs/target) databases. The transcript of GLS in each sample was extracted and log[2](TPM + 1) transformation was performed in these transcriptomic profiles. The abbreviations and expanded names of 33 tumor types are displayed in Table [56]1. The samples derived from tumor tissues and corresponding normal tissues were included for further pan-cancer analysis but the samples with missing gene expression or lacking required clinical information were eliminated. The proteomics datasets of GLS were derived from the Clinical Proteomic Tumor Analysis Consortium (CPTAC) ([57]https://pdc.cancer.gov/pdc/) and all proteomic expression values were converted into unitless Z-scores. Furthermore, the single-cell spatial transcriptome datasets were acquired from the Gene Expression Omnibus (GEO) ([58]https://www.ncbi.nlm.nih.gov/geo/) and 10x Genomics ([59]https://www.10xgenomics.com/) databases. The workflow in this study is demonstrated in Fig. [60]1. Table 1. The expanded names and abbreviations of 33 tumor types. Number Tumor types Abbreviations 1 Adrenocortical carcinoma ACC 2 Bladder urothelial carcinoma BLCA 3 Breast invasive carcinoma BRCA 4 Cervical squamous cell carcinoma and endocervical adenocarcinoma CESC 5 Cholangiocarcinoma CHOL 6 Colon adenocarcinoma COAD 7 Lymphoid neoplasm diffuse large B-cell lymphoma DLBC 8 Esophageal carcinoma ESCA 9 Glioblastoma multiforme GBM 10 Head and neck squamous cell carcinoma HNSC 11 Kidney chromophobe KICH 12 Kidney renal clear cell carcinoma KIRC 13 Kidney renal papillary cell carcinoma KIRP 14 Acute myeloid leukemia LAML 15 Brain lower grade glioma LGG 16 Liver hepatocellular carcinoma LIHC 17 Lung adenocarcinoma LUAD 18 Lung squamous cell carcinoma LUSC 19 Mesothelioma MESO 20 Ovarian serous cystadenocarcinoma OV 21 Pancreatic adenocarcinoma PAAD 22 Pheochromocytoma and paraganglioma PCPG 23 Prostate adenocarcinoma PRAD 24 Rectum adenocarcinoma READ 25 Sarcoma SARC 26 Skin cutaneous melanoma SKCM 27 Stomach adenocarcinoma STAD 28 Testicular germ cell tumors TGCT 29 Thyroid carcinoma THCA 30 Thymoma THYM 31 Uterine corpus endometrial carcinoma UCEC 32 Uterine carcinosarcoma UCS 33 Uveal melanoma UVM [61]Open in a new tab Fig. 1. [62]Fig. 1 [63]Open in a new tab Diagram of main workflow in this study. Differential expression analysis The mRNA expression of GLS between tumor and normal tissues was compared in TCGA samples with or without matched GTEx data across 33 tumor types. The co-expression correlation analysis across TCGA tumors was performed to clarify the correlation between GLS and the remaining 9 confirmed cuproptosis-related genes (CRGs). Furthermore, we investigated the total protein levels of GLS utilizing the CPTAC samples, and Wilcoxon rank sum tests were selected to compare protein expression differences between different groups. Lastly, the immunohistochemical (IHC) staining images were acquired to supplement and validate the GLS protein expression in different tissues through inquiring about the Human Protein Atlas (HPA) ([64]https://www.proteinatlas.org/) database, in which the quantity of stained cells (including not detected, < 25%, 25–75% and > 75%) and staining intensity (including negative, weak, moderate, and strong) were used to assess the expression levels of GLS protein. Association between GLS expression and clinical characteristics The relevance of GLS expression to distinct clinical features was investigated. We further detected the expression difference of GLS across the different demographic and oncological characteristics, including gender, age, pathological TNM stage, and histological grade across TCGA tumor samples. Furthermore, the “Gene expression comparison” module in TNMplot ([65]https://tnmplot.com/analysis/) was manipulated to analyze whether GLS expression was associated with tumor metastasis. The differential levels of GLS among normal, tumor, and metastatic samples were compared based on the RNA-Seq data and gene chip data. Differential expression analyses were measured via the Kruskal-Wallis test with a post-hoc Dunn’s test or One-way ANOVA with a post-hoc Tukey Honestly Significant Difference test when appropriate. Survival analysis We selected the overall survival (OS), disease-free interval (DFI), and progression-free interval (PFI) to assess the potential impact of GLS expression on survival outcomes. After excluding the patients with follow-up time or survival time less than thirty days, the potential prognostic role of GLS expression in TCGA and TARGET cohorts was analyzed. The “maxstat (v.0.7.25)” R package was utilized to calculate the optimal cutoff value of the GLS expression. The “survival (v.3.3.1)” and “survminer (v.0.4.9)” R packages were applied to analyze prognostic differences. The Kaplan-Meier (K-M) survival curves were delineated in diverse GLS expression subgroups and Cox regression was performed to determine whether these survival curves were significant. Additionally, the “K-M plotter” module in Kaplan-Meier Plotter ([66]https://kmplot.com/analysis/) was manipulated as a supplement to illuminate its prognostic effects on tumor patients. Single-cell analysis and spatial transcriptomic atlas Single-cell sequencing technology has become a crucial tool to explore tumor progression, heterogeneity, and therapeutic response. We dissected the tumor single-cell sequencing data to examine the relationship between GLS expression and 14 single-cell functional states across multiple tumor types via the CancerSEA ([67]http://biocc.hrbmu.edu.cn/CancerSEA/) platform, where the gene set variation analysis (GSVA) and Spearman’s correlation were employed separately to calculate the tumor single-cell functions and the correlation between these functions and GLS expression. The functional states with the most significance (p < 0.001) and corresponding expression distribution diagrams were also simultaneously obtained in specific single-cell samples in the CancerSEA database. With no STAD data available in the CancerSEA database, it was necessary to further elaborate on the biological roles of GLS in STAD. Hence, we acquired a single-cell transcriptome dataset ([68]GSE167297), which included 41,554 cells, downloaded from the TISCH ([69]http://tisch.comp-genomics.org/) database. To identify clusters and patterns in the STAD dataset, we applied the Uniform Manifold Approximation and Projection (UMAP) method to visualize high-dimensional data as a two-dimensional heatmap, and to visualize GLS expression at single-cell resolution. Meanwhile, the Kruskal-Wallis test was utilized to compare its biological expression among different cell types. According to whether these cells expressed GLS or not, we divided them into the GLS-positive expression group and GLS-negative expression group and calculated the proportion of each cell type in the two groups respectively. To evaluate the biological effects of GLS expression on tumor cells, various biological pathways related to immunity, metabolism, mitochondria, cell death, proliferation, and signaling were enriched in malignant tumor cells using the “AUCell (v.1.26.0)” R package and the referenced Hallmark and MitoCarta3.0 gene sets^[70]18,[71]19. The differences in pathway scores between the two groups of malignant tumor cells were compared via the “limma (v.3.60.0)” R package. Single-cell spatial transcriptome samples were obtained from the GEO ([72]GSE203612, [73]GSE175540, and [74]GSE179572) and 10x Genomics databases. Detailed sample sources are provided in Table [75]S1. To accurately identify the cell composition in each spot on the slides, deconvolution analysis was performed based on the referenced scRNA-seq data. Firstly, we constructed a signature matrix by calculating the average expression of specific expression genes for each cell type in each spot. Next, based on the “get_enrichment_matrix” and “enrichment_analysis” functions in the “Cottrazm (v.1.0.0)” R package, we generated the enrichment score matrix for cell composition identification. Then, the most abundant cell type and GLS expression landscape in each spot were calculated and visualized utilizing the “SpatialDimPlot” and SpatialFeaturePlot” functions in the “Seurat (v.5.0.3)” R package. Furthermore, a malignant spot was defined if its cells were exclusively malignant, and a normal spot was defined if its cells were exclusively non-malignant. Otherwise, the mixed malignant spot was defined. Wilcoxon rank sum tests were utilized to assess statistical differences in GLS expression between the two groups. Landscapes of genetic alteration and genomic heterogeneity The cBioPortal ([76]https://www.cbioportal.org/) was manipulated to elucidate the genetic alteration across the TCGA tumors, involved in genetic mutation, amplification, structural variant, and deep deletion, based on the “Query” searching module and “TCGA PanCancer Atlas Studies” datasets. Then the “Immune-Mutation” module in TIMER2.0 ([77]http://timer.cistrome.org/) was utilized to clarify the association of the GLS mutation with multiple immune cell infiltration. Lastly, to delineate whether genomic heterogeneity was relevant to GLS expression, we analyzed chromosome ploidy, homologous recombination deficiency (HRD), and mutant-allele tumor heterogeneity (MATH) based on the Sangerbox3.0 platform ([78]http://vip.sangerbox.com/), in which the ploidy and HRD data were derived from a previous study^[79]20, and MATH scores were calculated using the simple nucleotide variation datasets across TCGA tumor samples from TCGA portal based on the “inferHeterogeneity” function of “maftools (v.2.8.5)” R package. Analysis of antitumor drug sensitivity To uncover the association of GLS expression with antitumor drug sensitivity, we retrieved the pharmacogenomic datasets by searching the CellMiner database ([80]https://discover.nci.nih.gov/cellminer/home.do). The classic NCI-60 cancer cell line was selected for cytotoxic screening of antitumor compounds in vitro^[81]21. The drug activities for chemotherapies, targeted therapies, and endocrine therapies with Food and Drug Administration (FDA) approval were examined to search for potential antitumor therapeutic agents. Correlation analysis was performed between GLS mRNA expression (Z-value) and drug activity. Analysis of tumor immune microenvironment To investigate the association of GLS expression with tumor immune microenvironment, we performed a Single-sample Gene Set Enrichment Analysis (ssGSEA) algorithm for various immune cell infiltration analyses in TCGA tumor samples. Based on markers for various immune cells reported in previous literature, the “GSVA (v.1.46.0)” R package was applied to perform the ssGSEA analysis^[82]22. B cells, T cells, dendritic cells (DCs), eosinophils, mast cells, natural killer (NK) cells, macrophages, neutrophils, and cytotoxic cells were selected to analyze their infiltration levels. Additionally, the “Immune-Gene” module in TIMER2.0 was manipulated to explore the association between GLS expression and various infiltrated cells in the TME. CD4^+ T lymphocytes, CD8^+ T lymphocytes, cancer-associated fibroblasts (CAFs), endothelial cells (ECs), regulatory T cells (Tregs), and myeloid-derived suppressor cells (MDSCs) were filtrated for further correlation analysis. The default immune deconvolution methods in TIMER2.0 were utilized to calculate the infiltration abundances of various cells. Analysis of immunogenomics and immunotherapeutic response Immune checkpoints (ICPs), as immune regulators, act a pivotal part in maintaining tumor immune microenvironment. Targeting ICPs remains an especially considerable immunotherapy modality today. We acquired an immunomodulatory gene list from a previous study that included 24 inhibitory checkpoints^[83]20. The co-expression connections between GLS expression and these immunomodulatory checkpoints were analyzed across TGCA tumors through correlation analysis. Furthermore, in malignant tumors, the dysregulated chemokines as well as their receptors can change immune cell recruitment and activation, often towards a tumorigenic state, which has become a potential target for tumor immunotherapy^[84]23–[85]25. Thus, correlation analyses were performed to examine whether GLS expression was relevant to chemokines and their receptors across various TGGA cancer types. To illustrate the roles of GLS in immunotherapy, several tumor immunogenicity predictors, including microsatellite instability (MSI), tumor mutational burden (TMB), and tumor neoantigen that were acquired from Sangerbox3.0 platform, were selected to estimate the immunotherapeutic response in TCGA tumors. The radar map visualized the correlations between GLS expression and MSI, TMB, and tumor neoantigen. Since immunophenoscore (IPS) is a well-established and remarkable predictor of immunotherapy response, we categorized the TCGA tumor patients into different subgroups based on the GLS expression median and compared the IPS levels between the different subgroups. The IPS scores in TCGA samples were acquired from the TCIA ([86]https://tcia.at/home) database. A higher IPS score implies a more sensitive response to immunotherapy. Additionally, two independent immunotherapy cohorts (IMvigor210 and [87]GSE91061) in the GEO database were retrieved to validate the potential effect of GLS on tumor immunotherapy. The patients who developed complete or partial responses after immune checkpoint inhibitor (ICI) treatment were identified as immunotherapeutic responders and those who developed stable or progressive diseases were identified as non-responders. The differences in GLS expression between the responders and non-responders were explored, and relevant immunotherapeutic survival curves were obtained via the “Start KM Plotter for immunotherapy” module in the Kaplan-Meier Plotter database to infer immunotherapeutic prognosis in different GLS expression subgroups. GLS-related gene enrichment analysis The STRING ([88]https://cn.string-db.org/) was manipulated to acquire the top 20 interactors relevant to GLS based on the protein-protein interaction (PPI) network. “Full STRING network” was chosen in the settings options and active interaction sources were limited to “experiments” and “databases”. Meanwhile, the “Correlation Analysis” module in the GEPIA2 ([89]http://gepia2.cancer-pku.cn/) platform was utilized to obtain the top 100 similar genes of GLS through Pearson correlation analysis. Whereafter, a correlation heatmap of the selected genes was delineated across TCGA tumors to supply the correlativity between GLS and GLS-related genes utilizing the “Exploration- Gene_Corr” module of the TIMER2.0 platform. To discover the influence of GLS on pathways and biological processes, the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analysis were performed using the aforementioned 120 genes based on the “clusterProfiler (v.4.4.4)” R package. After that, we classified the TCGA tumor samples into different expression subgroups based on the GLS expression median. With |LogFC|>1.0 and adjusted p < 0.05 as the threshold value, the differential expression genes between the two subgroups were calculated via the “DESeq2 (v.1.36.0)” R package. The Gene Set Enrichment Analysis (GSEA) analysis subsequently was conducted via the “clusterProfiler” R package utilizing the differentially expressed genes in different TCGA tumor types. Gastric cancer population We collected the clinical information of 295 patients who were pathologically diagnosed with GC at the Department of Gastrointestinal Surgery of the First Hospital of Jilin University from January 2011 to December 2016. All enrolled patients underwent gastric surgery but did not undergo preoperative neoadjuvant chemoradiotherapy. We extracted clinical information of these GC patients, such as age, sex, body mass index, smoking history, drinking history, pathological TNM stages (AJCC 8th), histological grades, lympho-vascular invasion, perineurium invasion, and MSI, from the hospital’s medical record system. The clinical information regarding these GC samples is summarized in Table [90]S2. The data collection and further study were approved by the Ethics Committee of the First Hospital of Jilin University (No.2021 − 493) and relevant informed consent was obtained. Regular follow-up was performed in accordance with our previous study^[91]26. All follow-up works were completed until August 2022 and OS for each patient was calculated separately during follow-up. Immunohistochemical staining and scoring In the above GC population, 295 tumor tissues and 277 paired para-tumor tissues were acquired for IHC staining. Briefly, the paraffin sections were baked at 65 °C for 30 min, followed by dewaxing and rehydration. After repairing the antigen with sodium citrate (MVS-0101, Maixim, China), the sections were placed in a wet box and an endogenous peroxidase blocker (SP KIT-A3, Maixim, China) was added. Subsequently, normal non-immune goat serum (SP KIT-B2 Maixim, China) was added for antigen blocking. Then, the tissue sections were incubated with GLS monoclonal antibody (#56750, CST, USA) overnight at 4℃, followed by a secondary antibody (KIT-5020, Maixim, China) for 30 min. After color detection with 3,3-diaminobenzidine tetrahydrochloride (DAB-1031, Maixim, China), re-staining with hematoxylin and differentiation with hydrochloric alcohol were carried out. Finally, the tissue sections were dehydrated and mounted. The ASI system (Applied Spectral Imaging Ltd., Israel) was applied to estimate the staining intensity and the percentage of stained tumor cells in each tissue section. The staining intensity was classified as 0 (no staining), 1 (light yellow), 2 (light brown), or 3 (brown). According to the percentage of stained tumor cells (pi) and the staining intensity (i), the H-score (H-score = Σ(pi×i)) was calculated to assess GLS protein expression levels in different tissues. Subsequently, differences in GLS expression between tumor tissues and para-tumor tissues were tested based on H-score. Based on the optimal cutoff value of the H-score, these patients were categorized into different expression subgroups, and the OS was simultaneously visualized utilizing “survival” and “survminer” R packages. Cox regression was performed to compute the hazard ratio (HR) for OS. The distribution of pathological features between the different subgroups was subsequently assessed using the Chi-squared test. After univariate and multivariate Cox regression analyses were performed using “survival” R package, the clinical nomograms based on or without GLS H-score were constructed via “survival” and “rms (v.6.3.0)” R packages to conveniently predict prognostic survival in GC patients. The 3-, 5-, and 7-year time-dependent ROC curves and calibration curves were plotted to evaluate the nomogram’s predictive utility via “timeROC (v.0.4)”, “survival”, and “rms” R packages. Similarly, we plotted ROC curves for the clinical and pathological variables shown in the nomogram model to compare the difference in predictive accuracy between GLS-based risk model and these variables. Cell culture and GLS targeted siRNA transfection Immortalized gastric normal epithelial cell line GES-1 and GC cell lines, including AGS, MKN-1, and SNU-638, were provided by ZQXZbio (Shanghai, China) and Procell Life Science & Technology (Wuhan, China). These cells were cultured in RPMI-1640 medium with 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin solution and placed in a 37℃-incubator containing 5% CO[2]. Once the cell confluence reached 60–80% in culture plates, siNC or siGLS (GenePharma, China) was transfected into GC cells utilizing lipofectamine RNAiMAX (13778030, Invitrogen, USA) according to the manufacturer’s instruction. The oligonucleotide sequences for siRNAs are displayed in Table [92]S3. Briefly, Opti-MEM medium was used to dilute lipofectamine and siRNA (10µM), respectively. Subsequently, they were mixed and incubated at room temperature for 5 min to prepare the siRNA-liposomes that were then added to each well. Quantitative real-time PCR (qPCR) assays When GC cells were transfected for 48 h, we extracted cellular RNA using Molpure® Cell/Tissue Total RNA Kit (19221ES50, YEASEN, China). The concentrations of the extracted RNA were quantified utilizing a NANODROP 2000c spectrophotometer (Thermo Scientific, USA). Subsequently, the reverse transcription process was performed in C1000™ Thermal Cycler (BIO-RAD, Singapore) device utilizing the ToloScript RT EasyMix for qPCR Kit (#22106, Tolo Biotech, China). All primers used for qPCR were synthesized following the sequence in Table [93]S4. The qPCR reaction was performed utilizing the 2×Q3 SYBR qPCR Master Mix (#22204, Tolo Biotech, China) premixed solution and LightCycler® 480II Instrument (Roche, Switzerland). Additionally, differential expression of GLS mRNA between normal GES-1 and GC cells was also performed following the qPCR procedure described above. Western blotting assays Total protein of GC cells was extracted using RIPA Lysis Buffer (SL1010, CooLaber, China) and quantified using BCA Protein Assay Kit (ZJ101, Epizyme, China). Protein electrophoresis was performed using SDS-PAGE Gel Rapid Preparation Kit (PG112, Epizyme, China) and Tris-Gly buffer system. Subsequently, the proteins in the gel were transferred to the 0.45 μm pore size polyvinylidene fluoride (PVDF) membranes (IPVH00010, MILLIPORE, USA). After sealing with 5% skim milk for 2 h, the PVDF membranes were incubated overnight at 4℃ in anti-GLS antibody (#56750, CST, USA) and anti-β-Actin antibody (abs171598, Absin, China). The next day, the secondary antibody labeled with HRP (H6162, UElandy, China) was incubated at room temperature for binding to primary antibody. The ECL Femto Light Chemiluminescence Kit (SQ201, Epizyme, China) was used for chemiluminescence detection of protein bands. Cell counting Kit-8 (CCK-8) assays CCK-8 assays were applied to explore the role of GLS knockdown in GC cell proliferation. GC cells were evenly spread into 96-well plates with 1 × 10^4 cells per well and cultivated to the desired confluence. The cell viability was detected at 24 h, 48 h, and 72 h after siRNA transfection. 10 µl of CCK-8 (C6005, UElandy, China) solution was added to each well and incubated in the incubator sheltered from light for 1 h. Ultimately, the absorbance of each well was detected at 450 nm wavelength utilizing a microplate reader (Thermo Scientific, USA). Wound healing assays 1 × 10^6 GC cells were spread into each well of 6-well plates after transfection. Once their confluence was close to 100%, the cell layer was scratched vertically utilizing the 200 µl pipette tip and photographed under a microscope immediately. After these GC cells were cultivated in an FBS-free medium for 48 h, the scratched sites were photographed again. The ImageJ (v.1.46r) software was utilized to calculate the cell migration distance between the corresponding photographs. Transwell migration assays After transfection for 24 h with serum-free culture, we seeded 1 × 10^5 GC cells into the upper chamber (3422, Corning, USA) with 200 µl of serum-free RPMI-1640 medium. Accordingly, we added 600 µl of medium containing 20% FBS into the lower transwell chamber. While GC cells were cultured for 24 h, we lightly wiped the GC cells in the upper chamber using cotton swabs, and then 4% paraformaldehyde was used to immobilize GC cells entering the lower chamber. Finally, the migrating cells were stained using 0.1% crystal violet. Five random microscopic fields were photographed and ImageJ software was used to estimate the cell number in each field. Matrigel invasion assays A 10-fold dilution of the Matrigel (356234, Corning, USA) was performed using a cold FBS-free medium. Next, we added 100 µl of diluted Matrigel into the upper transwell chamber and incubated in the incubator for at least 3 h to form a matrix film in the upper chamber. Subsequently, we seeded 1 × 10^5 transfected GC cells with 200 µl of FBS-free medium into the upper chamber. The subsequent operation was the same as the description in the transwell migration assay. Cell apoptosis assays GC cells were seeded into 6-well plates with 2 × 10^5 cells per well and cultivated for 24 h. Subsequently, these cells were transfected by siRNA and continued to be cultured for 48 h. The YF^®488-Annexin V and PI Apoptosis Kit (Y6002, UElandy, China) was utilized to assess the cell apoptosis levels based on the manufacturer’s instruction. In brief, GC cells were digested using EDTA-free trypsin and collected in the centrifuge tubes, then these cells were recovered in optimal cell culture conditions and medium for 30 min. Then they were washed twice using cold PBS and were re-suspended with Annexin V binding buffer. After being filtered by 40 μm cell screens, they were stained using YF^®488-Annexin V and PI and incubated at room temperature for 15 min away from light. Ultimately, 200 µl of Annexin V binding buffer was supplemented before detecting the apoptosis via flow cytometry (BD, USA). Statistical analysis The normality of distribution and homogeneity of variance of continuous data were examined via the Shapiro-Wilk test and Levene’s test, respectively, since they determined which statistical method was the optimal option. Wilcoxon rank-sum test or independent-sample T test was selected to test the statistical differences between two different subgroups. Similarly, Kruskal-Wallis test or One-way ANOVA was applied to analyze the significant differences among multiple subgroups. Subsequently, a post-hoc test was performed for pairwise comparisons by Dunn’s test or Tukey Honestly Significant Difference test when appropriate. Wilcoxon signed-rank test was performed to test the statistical significance in paired samples. Correlation analyses were measured through the Spearman or Pearson test when appropriate. The K-M curves were delineated based on diverse GLS expression subgroups categorized by the optimal cutoff value. Cox regression was performed to test the statistical differences in survival curves. Univariate and multivariate Cox regression analyses were applied to identify prognostic risk factors. Two-tailed p < 0.05 was considered statistically significant. The R software (v.4.2.0) was utilized for statistical analysis and the “ggplot2 (v.3.3.6)” R package was utilized to visualize analyzed data. Results GLS is abnormally expressed in human cancers Through comparing transcription data between tumor and normal samples in the TCGA database, we discovered GLS expression was notably upregulated in various malignancies, including CHOL, COAD, ESCA, HNSC, LIHC, and STAD, whereas that was significantly downregulated in BRCA, GBM, KICH, KIRC, KIRP, LUAD, LUSC, PCPG, and UCEC (Fig. [94]2A). After matching TCGA samples with GTEx data, upregulated GLS were found additionally in DLBC, LAML, PAAD, READ, and THYM, but downregulated GLS were observed additionally in ACC, BLCA, CESC, LGG, OV, PRAD, SKCM, THCA, and UCS (Fig. [95]2B). Additionally, differential expression of paired samples between tumor tissues and adjacent tissues was observed in BRCA, CHOL, ESCA, HNSC, KICH, KIRC, LIHC, LUSC, PRAD, STAD, THCA, and UCEC (Fig. [96]2C). Overall, we were surprised to notice GLS was upregulated in the majority of digestive malignancies, while that was downregulated in almost all urogenital malignancies. Notably, we unveiled a significant positive co-expression relationship between GLS and currently known CRGs except the CDKN2A gene (Fig. [97]2D). In the HPA database, compared with normal tissues, higher IHC staining was observed in colon cancer, cutaneous melanoma, thyroid cancer, pancreatic carcinoid, lung cancer, and lymphoma (Fig. [98]2E). In terms of the total protein levels, the analysis results based on CPTAC samples uncovered the significantly upregulated protein expression of GLS in COAD, HNSC, and LIHC, but downregulated expression in GBM, KIRC, and UCEC (Fig. [99]2F). These differences were consistent with the trend of GLS mRNA expression in the corresponding tumors. Fig. 2. [100]Fig. 2 [101]Open in a new tab Differential expression of GLS in pan-cancer. (A) GLS mRNA expression between tumor tissues and unpaired normal tissues in TCGA samples. (B) GLS mRNA expression between tumor tissues and unpaired normal tissues in the matched TCGA and GTEx samples. (C) GLS mRNA expression between tumor tissues and paired normal tissues in TCGA samples. (D) Co-expression relationship between GLS and other 9 cuproptosis-related genes in TCGA cancers. (E) GLS protein expression in tumor and normal tissues based on the immunohistochemical images of HPA database. (F) Total protein levels of GLS between tumor and normal tissues in CPTAC samples. *p < 0.05; **p < 0.01; ***p < 0.001. GLS expression is relevant to clinical characteristics in human cancers With the assistance of the TNMplot tools, we attempted to explore the connection between GLS expression and tumor-metastatic progression. We noticed a trend toward higher GLS expression in metastatic tumor tissues compared with normal and primary tumor tissues in esophageal, hepatocellular, and skin cutaneous melanoma, as well as thyroid cancers (Fig. [102]3A). After that, we made a thorough inquiry about the expression differences of GLS in further detailed clinical characteristics. Interestingly, we realized GLS expression levels were upregulated in KIRC, KIRP, LIHC, and LUAD in women compared with that in men (Fig. [103]3B). We also found GLS expression levels in BRCA, KIRP, LIHC, and PAAD were increased in the patients aged 60 years or less compared with that in the older (Fig. [104]3C). Furthermore, we uncovered that GLS expression was associated with the pathological features in many cancers. Specifically, the differential expression of GLS was detected in KIRP, LIHC, PRAD, READ, and THCA based on the categories of T stage (Fig. [105]3D). Likewise, that was found in KIRP, LIHC, PRAD, and THCA based on the categories of N stage (Fig. [106]3E). Also, we indicated the significant connection between GLS expression and histological grades in LIHC and UCEC (Fig. [107]3F). Additionally, we clearly observed GLS expression was associated with TNM stages in KIRP, LIHC, THCA, and UCEC (Fig. [108]3G). Taken together, these findings indicated a strong association between GLS expression and advanced clinicopathology in multiple tumors, especially in LIHC, KIRP, and THCA. Fig. 3. [109]Fig. 3 [110]Open in a new tab Relationship between GLS expression and clinical characteristics in pan-cancer. (A) Association between GLS expression and tumor metastasis in TNMplot database. (B-G) Differential expression of GLS in different genders (B), ages (C), T stages (D), N stages (E), pathological grades (F), and TNM stages (G) based on TCGA data. GLS is a prognostic biomarker in human cancers In TCGA and TARGET cohorts, we discovered high GLS expression was correlated with worse OS in STAD, BRCA, ESCA, HNSC, KIRP, LAML, LGG, LIHC, MESO, and UCEC, worse PFI in BLCA, CESC, COAD, LUSC, PRAD, SKCM-primary (SKCM-P), and UVM, and worse DFI in KIPAN (KICH, KIRC, and KIRP), PAAD, and THCA (Fig. [111]4). In the Kaplan-Meier Plotter database, we noticed that the patients with higher GLS expression in gastric, liver, colon, and ovarian cancer had significantly worse OS (Fig. [112]S1). Collectively, high GLS expression was well established to predict poor prognosis in these 20 tumor types, including BLCA, BRCA, CESC, COAD, ESCA, HNSC, KIRP, LAML, LGG, LIHC, LUSC, MESO, OV, PAAD, PRAD, SKCM, STAD, THCA, UCEC and UVM. Fig. 4. [113]Fig. 4 [114]Open in a new tab Survival analysis of GLS expression in pan-cancer, including OS, PFI, and DFI. Single-cell functional analysis and spatial transcriptomic atlas of GLS in human cancers Single-cell sequencing analysis was conducted via the CancerSEA database to examine whether GLS was relevant to different functional states of tumor cells at single-cell resolution. As displayed in Fig. [115]5A, GLS expression was positively relevant to tumor stemness, epithelial-mesenchymal transition (EMT), and inflammation in multiple tumor types, but these correlations were relatively weak. Detailed relevance between GLS expression and 14 functional states across distinct cancers is displayed in Table [116]S5. Then, we investigated the correlation between GLS and the functional states in specific single-cell tumor samples. The results indicated GLS expression levels were strongly positively relevant to metastasis, inflammation, EMT, and proliferation in AML, positively relevant to angiogenesis and inflammation in retinoblastoma (RB), and negatively relevant to DNA repair, DNA damage, and apoptosis in UVM (Fig. [117]5B). Meanwhile, the T-SNE diagrams showed GLS expression distribution at single-cell resolution in LAML, RB, and UVM. Fig. 5. [118]Fig. 5 [119]Open in a new tab Analysis of tumor single-cell sequencing in pan-cancer. (A) Correlation between GLS expression and 14 single-cell functional states. (B) T-SNE expression diagram and the most significantly correlated single-cell functional states in AML, RB, and UM samples. (C) Cell type annotation and GLS expression identification by UMAP in STAD. (D) Differential expression of GLS among different lineages in STAD. (E) Proportion of each cell type in GLS-positive and GLS-negative expression groups in STAD. (F) Biological pathways with significant differences between GLS-positive and GLS-negative groups in malignant cells of STAD. ***p < 0.001. Since the above findings didn’t address the biological roles of GLS in STAD, single-cell analyses were performed specifically to elucidate that. As demonstrated in Fig. [120]5C, different cell populations in STAD were clearly distinguished in the UMAP diagram, and differential GLS expression among different cells was shown. We found GLS expression was relatively higher in malignant cells among different cell lineages (Fig. [121]5D). Meanwhile, we noticed the proportions of CD8^+ T lymphocytes, plasma cells, and fibroblasts were reduced by about half in the GLS-positive group compared with the GLS-negative group in STAD (Fig. [122]5E). Furthermore, some notorious cancer-related pathways were significantly enriched in STAD. The immunological pathways such as “IL6 JAK STAT3 Signaling” and “Inflammatory Response”, the metabolic and mitochondrial pathways related to glutamine metabolism and oxidative phosphorylation, and the proliferation and cell death pathways such as “Myc Targets”, “Cuproptosis” and “Apoptosis”, as well as multiple key cancer-related signaling such as “PI3K Akt MTOR Signaling”, “MTORC1 Signaling”, “Hypoxia”, “EMT” and “KRAS Signaling”, were significantly observed in GLS-positive malignant cells (Fig. [123]5F). These findings imply GLS may have the potential to regulate tumor immunity, metabolism, proliferation, progression, and cell death, which warrants further investigation and experimental validation. Meanwhile, we elucidated the spatial expression patterns of GLS in human tumor samples based on single-cell transcriptomic and spatial transcriptomic data. As demonstrated in Fig. [124]6, we found GLS expression had a highly similar spatial localization to that of tumor cells in BRCA, colorectal cancer (CRC), KIRC, LUAD, LUSC, OV, SKCM, and gastrointestinal stromal tumors (GIST). In these tumor tissues, higher GLS expression levels were observed in malignant and mixed-malignant spots compared to the spots composed of non-malignant cells, in which GLS expression levels were the lowest. These findings imply that the spatial expression of GLS in TME is mainly derived from tumor cells, and therefore tumor cells may have a stronger ability to utilize glutamine than non-tumor cells such as immune cells. Fig. 6. [125]Fig. 6 [126]Open in a new tab Spatial transcriptome expression atlas of GLS in human cancers. GLS expression is associated with genetic alteration and genomic heterogeneity We queried the information on GLS alteration across TCGA pan-cancer cohorts through retrieving the cBioPortal database. As indicated in Fig. [127]7A, its alteration frequency in UCEC was the highest (5.86%), and mutation and amplification were the two dominating types in all TCGA tumors. As displayed in Figs. [128]7B and 103 mutation sites were identified in its protein structure where the genetic missense was the most common mutation type. The mutation sites seemed to be mainly concentrated in the glutaminase component segment (244–530). For example, we observed the occurrence of missense or nonsense mutations (amino acid change: R387Q) with underlying clinical significance in three UCEC patients. Their detailed mutation sites were shown in the stereoscopic protein structure of GLS (Fig. [129]7C). To further investigate the significance of GLS mutation for UCEC patients, we explored the association between GLS mutation and multiple infiltrated immune cells via the TIMER2.0 database. Obviously, the cohorts with mutated GLS in UCEC samples showed more abundances of CD8^+ T cells, B cells, myeloid dendritic cells, and M1 macrophages than the cohorts with wild-type GLS (Fig. [130]S2). This could mean GLS mutation affects its enzymatic metabolic activity, resulting in changes in the immune cell infiltration and tumor immune microenvironment. Unfortunately, there is currently a lack of data to elucidate the roles of GLS mutations in UCEC prognosis. Fig. 7. [131]Fig. 7 [132]Open in a new tab Landscapes of genetic alteration and genomic heterogeneity in pan-cancer. (A) GLS alteration frequency across TCGA pan-cancer cohorts of cBioPortal database. (B) 103 mutation sites in the protein structure of GLS based on the cBioPortal database. (C) The stereoscopic protein structure of GLS and mutation site (AA change: R387Q) in UCEC. (D-F) Roles of GLS alteration in the genomic heterogeneity, including chromosome ploidy (D), HRD (E), and MATH (F). Moreover, the relevance between GLS and genomic heterogeneity was analyzed. We witnessed GLS expression had a significant positive relevance to chromosome ploidy in COAD, STAD, and UCEC (Fig. [133]7D), suggesting that high-GLS expression may cause abnormal chromosome numbers in tumor cells. We also discovered GLS expression had a significant positive relevance to HRD in eight tumors, including BRCA, COAD, HNSC, KIRP, LIHC, PRAD, STAD, and UCEC, implying that the homologous recombination repair of these tumor cells is dysfunctional and patients with high-GLS expression may be highly sensitive to poly (ADP-ribose) polymerase (PARP) inhibitors and platinum-based chemotherapy agents (Fig. [134]7E). Simultaneously, we found higher GLS expression predicted higher MATH, that is, the greater tumor heterogeneity, in COAD, LIHC, LUAD, STAD, and UCEC (Fig. [135]7F). Overall, these results demonstrated GLS expression levels were dramatically relevant to increased genomic instability and tumor heterogeneity, especially in STAD, COAD, LIHC, and UCEC, which may be one of the reasons for their worse prognosis. GLS expression is relevant to the therapeutic sensitivity of antitumor drugs We analyzed whether GLS expression levels were relevant to the sensitivity of antitumor drugs utilizing the pharmacogenomic datasets in the CellMiner database where we selected the classical NCI-60 cancer cells to detect the activities of drugs in vitro. The findings indicated GLS expression was relevant to enhanced activities of various targeted drugs, such as afatinib, dacomitinib, erdafitinib, erlotinib, gefitinib, ibrutinib, lapatinib, midostaurin, neratinib, and vandetanib. However, GLS expression was associated with therapeutic resistance to several common chemotherapy and endocrine therapy drugs, such as actinomycin D, docetaxel, fulvestrant, vinorelbine, and paclitaxel (Fig. [136]8). These results provide substantial implications and references for the selection of clinical treatment