Spatial transcriptomics reveals unique metabolic profile and key oncogenic regulators of cervical squamous cell carcinoma

Abstract

Background

As a prevalent and deadly malignant tumor, the treatment outcomes for late-stage patients with cervical squamous cell carcinoma (CSCC) are often suboptimal. Previous studies have shown that tumor progression is closely related with tumor metabolism and microenvironment reshaping, with disruptions in energy metabolism playing a critical role in this process. To delve deeper into the understanding of CSCC development, our research focused on analyzing the tumor microenvironment and metabolic characteristics across different regions of tumor tissue.

Methods

Utilizing spatial transcriptomics (ST) sequencing technology, we conducted a study on FFPE (formalin-fixed paraffin-embedded) tumor samples from CSCC patients. Coupled with single-cell RNA sequencing (scRNA-seq) data after deconvolution, we described spatial distribution maps of tumor leading edge and core regions in detail. Tumor tissues were classified into hypermetabolic and hypometabolic regions to analyze the metabolism profiles and tumor differentiation degree across different spatial areas. We also employed The Cancer Genome Atlas (TCGA) database to examine the analysis results of ST data.

Results

Our findings indicated a more complex tumor microenvironment in hypermetabolic regions. Cell-cell communication analysis showed that various cells in tumor microenvironment were influenced by the signalling molecule APP released by cancer cells and higher expression of APP was observed in hypermetabolic regions. Furthermore, our results revealed the correlation between APP and the transcription factor TRPS1. Both APP and TRPS1 demonstrated significant effects on cancer cell proliferation, migration, and invasion, potentially contributing to tumor progression.

Conclusions

Utilizing ST, scRNA-seq, and TCGA database, we examined the spatial metabolic profiles of CSCC tissues, including metabolism distribution, metabolic variations, and the relationship between metabolism and tumor differentiation degree. Additionally, potential cancer-promoting factors were proposed, offering a valuable foundation for the development of more effective treatment strategies for CSCC.

Introduction

Cervical cancer exhibits high incidence and mortality rates, ranking fourth among all cancers in women and approximately 342,000 women succumb to this cancer annually [1]. Cervical squamous cell carcinoma (CSCC) has emerged as the most prevalent histological type, constituting about 75-80% of cervical cancers. Currently, the main treatment modalities for cervical cancer include surgery and radiation therapy. While localized CSCC carries a favorable prognosis, the five-year survival rate stands at a mere 16.5% for patients with late stage, metastatic or recurrent cancers [2]. Since Human papillomavirus (HPV) infection is the primary cause, HPV vaccination is currently the most direct and efficacious strategy for preventing cervical cancer [3,4,5]. However, despite its effectiveness, the protective efficacy of vaccination remains limited, and ensuring high vaccination rates poses a challenge.

Metabolic reprogramming in cancer cells, which adjusts metabolism to promote cell growth and proliferation, is emerging as a key feature of cancer [6, 7]. Within the tumor microenvironment, different cell types engage in metabolic crosstalk, allowing cancer cells to thrive in challenging conditions [8, 9]. These metabolic adaptations not only enable cancer cells to sustain rapid growth and survival under nutrient-deprived and hypoxic conditions but also influence the surrounding microenvironment by shaping immune responses, angiogenesis, and extracellular matrix remodeling [10,11,12]. Consequnently, understanding tumor metabolism and microenvironment reshaping is critical for elucidating the mechanisms underlying CSCC progression. Chen et al. conducted RNA sequencing (RNA-Seq) of HeLa cells and demonstrated that abnormal glucose metabolism can accelerate the progression of cervical cancer [13]. Kim et al. analyzed The Cancer Genome Atlas (TCGA) database and demonstrated the significant role of various molecules in disrupting glucose metabolism, regulating the tumor immune microenvironment, and triggering immune evasion in cancer cells [14]. Utilizing TCGA data and clinical information, Dai et al. demonstrated that aerobic glycolysis was regulated in cervical cancer cells, promoting cancer cell growth and proliferation by accelerating cell cycle progression [15]. Growing evidence suggests that disruptions in energy metabolism are crucial for cervical cancer progression [16, 17]. However, most of these studies rely on cell experiments and bulk RNA-seq data, which offer an overall gene expression profile of cell populations across different tissue types. Given the complexity of metabolic networks, the heterogeneity of the tumor microenvironment, and the diversity of intercellular communication, investigating the relationship between visualizing metabolic reprogramming and cell-cell interactions remains a challenging task.

Single-cell RNA sequencing (scRNA-seq) enables characterization of transcriptomic profiles in human cancers with single-cell resolution and helps to identify cell populations with distinct biological functions [18, 19]. Spatial transcriptomics (ST), is a novel biotechnology that provides a comprehensive analysis of gene expression patterns and spatial distribution of various cell types in different tissues, addressing the spatial information deficiency in scRNA-seq data [20,21,22]. Ou et al. utilized Stereo-seq and spatial transcriptome technology to investigate cervical cancer cases, revealing higher transcriptional and translational activity, cell proliferation, oxidative phosphorylation, and immune response in tumor regions compared to other areas [23]. Fan et al. conducted a multi-omics analysis of CSCC, establishing a correlation between changes in the tumor immune microenvironment and different cell states [24]. Beyond cervical cancer, ST technology has significantly advanced research in various types of cancers. Wu et al. utilized Stereo-seq to identify a new invasion area in liver cancer patients, providing insights for the development of treatments for advanced liver cancer and other solid tumors [25]. The human body is a complex and highly organized organism, with intricate metabolic patterns and immune responses. Malignant tumors, while still challenging to understand fully, have a significant impact on human health. With the development of technology, we are now able to gain a better understanding of tumor heterogeneity, metabolic profiles, and tumor environments through bioinformatics analysis [26]. Furthermore, we can investigate their influence and mechanisms on pathophysiological functions of tumors through cellular and molecular experiments. In this study, we adopted a comprehensive approach that combined bioinformatics analysis and experimental techniques, including scRNA-seq, ST and in vitro cellular and molecular experiments. Our analysis delved deeply into the metabolism of tumor regions in CSCC cases, investigating the diverse metabolic patterns across different spatial locations and their impact on the tumor microenvironment. Our goal is to offer valuable insights for the development of novel anti-tumor treatment strategies that target metabolic reprogramming.

The APP gene, widely recognized for its significant association with Alzheimer’s disease, has predominantly been studied within the realm of neurodegenerative disorders [27]. However, emerging research has highlighted its role in various signalling pathways that affect the survival and proliferation of cancer cells [28,29,30]. TRPS1 has been identified as a crucial regulator of transcriptional networks involved in tumor progression and differentiation [31, 32]. Despite these findings, there is currently no research explicitly investigating whether APP and TRPS1 serve as valuable targets for tumor progression in CSCC or their potential therapeutic implications. This knowledge gap underscores the need for further exploration of these factors in CSCC.

In this study, ST analysis was conducted on two CSCC tissues and one normal tissue, along with scRNA-seq data from cervical cancer samples. The study elucidated the spatial distribution and immune microenvironment characteristics of the tumor leading-edge (LE) region and tumor core (TC) region, as well as investigated the correlation between tumor differentiation and metabolism. The findings revealed that endothelial cells, CD74 + macrophages, and other immune cells in the tumor microenvironment might be influenced by signalling molecules released by cancer cells, such as APP. Moreover, we observed a positive and significant association between APP expression and metabolic activity in cancer cells, and both APP and the transcription factor TRPS1 significantly promoted cancer cell proliferation, migration, and invasion. In summary, our analysis data demonstrated the close association between the metabolic spatial profiles of CSCC tissues and the molecules in the tumor microenvironment. Through in vitro cellular experiments, we demonstrated that these molecules play an essential role in the development of CSCC. Our research suggests that the metabolic spatial features of CSCC are closely related to its pathogenesis and progression, offering valuable insights for the prognosis and treatment of CSCC.

Materials and methods

Single-cell RNA-seq data processing

The sequencing data were processed using Cell Ranger (v.7.0.1) against the GRCh38 human reference genome. Subsequent quality control was applied to the obtained cells based on unique molecular identifier counts, the number of genes detected, and the proportions of mitochondrial, red blood cell-associated, and ribosomal gene counts per cell using Seurat (v.5.0.1) [33].

Specifically, cells meeting the following criteria were retained for further analysis: more than 200 and fewer than 10,000 genes, unique molecular identifier counts less than 60,000, a proportion of mitochondrial gene counts less than 25%, a proportion of red blood cell-associated gene counts less than 25%, and a proportion of ribosomal gene counts less than 0.05%. Library size normalization was performed in Seurat on the filtered matrix to obtain the normalized count, which was then used for all downstream analysis.

Differential expression analysis was performed using the FindAllMarkers function to identify significantly differentially expressed genes across clusters, with a threshold of p-value < 0.05 (Wilcoxon rank-sum test). Additionally, the top_n function was utilized to select the top 20 and top 5 marker genes for each cluster based on average log2 fold change.

Dimension reduction and cell-type annotation

Dimension reduction and unsupervised clustering were applied to the single-cell data following the workflow in Seurat. Briefly, 3,000 highly variable genes were selected for downstream analysis by using the Find Variable Features function with the parameter ‘nFeatures = 3000.’. The harmony package (v.1.1.0) [34] was used for scRNA-seq data integration. Principal component analysis (PCA) was performed on an integrated data matrix to reduce the dimensionality of the scRNA-seq dataset. Using the Elbow plot function of Seurat, the top 30 PCs were subjected to the downstream analysis. The main cell clusters were identified using the Find Clusters function in Seurat, with the resolution set as 0.8. Subsequently, the clusters were visualized by constructing UMAP plots. Conventional markers described in a previous study were used to categorize each cell into a known biological cell type.

Stereo-seq library preparation and sequencing

We collected one normal sample and two CSCC samples from two patients. The two CSCC samples were respectively located at the tumor core and tumor border. Spatially resolved transcriptomes were generated from these samples using Stereo-seq following established protocols. For each FFPE sample, it was cut into two consecutive sections with a thickness of 10 μm. One section was attached to a glass slide and stained with H&E following the previous protocol. The second section was mounted on a Stereo-seq N transcriptomics chip (Cat# 211KN114-EA, STOmics). RNAs were captured and reverse-transcribed using the Stereo-seq procedure following the vendor’s manual and prior publications [21, 35, 36]. Firstly, after drying the tissue-mounted Stereo-seq chip slide, paraffin was melted at 60 °C for 1 h and removed through deparaffinization using a xylene substitute and ethanol. Following the complete removal of paraffin, de-crosslinking of the tissue on the slide-chip was conducted using the STOmics reagent kit (Cat# 211KN114, STOmics). Next, the chip was immersed in pre-cooled methanol for tissue fixation at -20 °C for 17 min. After drying the chip, the tissue section was incubated in permeabilization buffer at 37 °C for 30 min. Subsequently, the FFPE Dimer Mix was added, and the chip was incubated at 25 °C for 1 h. Captured RNAs were then reverse-transcribed and ligated to probes adhered to the surface of the Stereo-seq chip through overnight incubation at 42 °C. The cDNAs were subsequently released from the chip, followed by size selection, amplification, and purification. The cDNA concentration was quantified, and finally, sequencing libraries were constructed using 20ng of cDNA from each sample with a library preparation kit (Cat# 111KL114, STOmics) and sequenced on the DNBSEQ-G400 sequencing platform.

Stereo-seq raw data processing

Fastq files generated by DNBSEQ-G400 were paired-end. Read 1 contained CID (1–25 bp) and MID (26–31 bp), while read 2 consisted of the cDNA sequences. The BGI STOmics analytical pipeline SAW (https://github.com/STOmics/SAW) was used to process the raw FASTQ files. Briefly, CID sequences from the first reads were mapped to the designed coordinates of the in situ captured chip using the sequence-to-coordinate correspondence stored in the mask file, allowing for one base mismatch. Reads with low-quality MID (containing N bases or more than 2 bases with a quality score lower than 10) or adapter sequences were filtered out. Retained reads were aligned to the reference genome hg38 and annotated to their corresponding genes. Duplicated reads, which contain the same CID, MID (allowing 1 mismatch to correct for sequencing and PCR errors), and the same gene locus, were collapsed. Finally, the gene expression profile matrix with coordinate information was generated. The expression profile matrix was divided into nonoverlapping bins covering 50 × 50 DNB area (bin50) for subsequent analysis. All bins are remained without any filter.

Annotation of spots in stereo-seq slides with deconvolution analysis

To comprehensively characterize the spatial transcriptomic landscape within a tumor region, a bin size of 50 (50 × 50, 250 × 250 μm) was utilized as a fundamental unit (spot) for subsequent analysis (Supplementary Table S2). To spatially map cell types as defined by scRNA-seq analysis within the Stereo-seq data, we employed cell2location model (v.0.5-alpha) [37] with hyperparameters N_cells_per_location = 6 and detection_alpha = 20. Specifically, we used cell2location to decompose multi-cell spatial transcriptomics data into spatially resolved estimates of cell-type abundance. First, this model derives cell-type expression signatures by calculating the average expression of each gene across cell types from the scRNA-seq data. Next, it performs hierarchical non-negative decomposition of the gene expression profiles at spatial locations (spots with multiple cells) to match them with the reference signatures, thereby identifying cell-type distributions. This tool effectively decomposes multi-cell ST data, providing spatially resolved estimates of cell-type abundance. Each bin was annotated with the most abundant cell type within it.

Division of regions into high and low metabolic regions

We assessed the activity of chosen pathways through two distinct approaches. Firstly, we employed the scanpy.tl.score_genes() function in Scanpy (v.1.9.3) [38] to calculate gene set scores within ST data, visualizing the results using the scanpy.pl.spatial() function. Additionally, we computed scores for metabolic gene sets using the GSVA (v.1.50.0) [39], and presented the outcomes using the ComplexHeatmap (2.18.0) [40]. We calculated the scores of metabolic pathways for each spatial region from the ST data and utilized these scores as biological features for the analysis of metabolic activity in tumor samples. To better understand the metabolic heterogeneity within the tumor, we divided the LE and TC regions into 2 and 3 parts, respectively, with equal area for each part. The gene sets of the six metabolic pathways mentioned above were merged, and metabolic scores of representative cancer cells in each part were calculated. Initially, quartiles and medians of the metabolic pathway scores were computed for each part to analyze the distribution of metabolic activity within each part. Based on the comparison between the median score of each part and the overall median score, regions with scores higher than the overall median were classified as hypermetabolic regions, while those with scores lower than the overall median were classified as hypometabolic regions. To address concerns about insensitivity around the median, regions within ± 10% of the overall median were categorized as intermediate metabolic regions.

Chromosomal CNV and transcription factor analysis

Inferred copy number variation (CNV) analysis was carried out using the InferCNV (v.1.18.0) (https://github.com/broadinstitute/inferCNV). The expression intensity of genes across chromosomal positions for each area cancer cells were compared against the mean value calculated for the endothelial cells and fibroblasts, and displays the relative gene expression on each chromosome in the form of a heatmap. Analysis of transcription factor regulatory network in each region of cancer cells was performed using SCENIC R package (v.1.3.1) [41].

Constructing cellular trajectories using RNA velocity

The RNA velocity analysis was performed using scVelo (v.0.2.5) [42] in Python. After gene selection and normalization, the first- and second-order moments were calculated using the scv.pp.moments() function. The velocities were then obtained using the scv.tl.velocity() function in dynamical mode. These velocities were projected onto diffusion maps and visualized as streamlines using the scv.pl.velocity_embedding_stream() function.

Differential expression and gene set enrichment analysis

The FindMarkers function in Seurat, taking the Wilcoxon rank-sum test as the statistical method and applying a minimum expression fraction of 10% (min.pct = 0.1), was used to analyze the differences between samples or cell types. Differentially expressed genes were identified based on an adjusted p-value < 0.05 and an absolute log2 fold change ≥ 1.5. The genes for difference analysis were analyzed using the clusterProfiler (v.4.10.0) [43], R package for Gene Ontology (GO), Kyoto Encyclopedia of Genes (KEGG), and GSEA enrichment [44]. The results of differential expression analysis and enrichment analysis were visualized using ggplot2 (v.3.4.4) (https://ggplot2.tidyverse.org/).

Ligand–receptor interactome

The cell–cell interaction analysis was performed by CellPhoneDB (v.4.0.0) [45] based on receptor-ligand interactions between two cell types/subtypes. Permutation test repeated 1000 times was used to evaluate the significance for ligand–receptor pairs across each cell type. Individual ligand or receptor expression was threshold with a cutoff value based on the average log gene expression distribution for all genes across all the cell types, ensuring that only interactions where at least 5% of the cells expressed the involved ligand or receptor were considered. The threshold for significance was set at P < 0.05. Statistical evaluations of the interactions were conducted using the cpdb_statistical_analysis_method from CellPhoneDB, which generated significant means and deconvoluted results. The significant interactions were then extracted using the search_analysis_results function, enabling the identification of meaningful communication networks within the cellular microenvironment.

Bulk RNA-seq analysis

The bulk RNA-seq data of CSCC samples were downloaded from the TCGA CESA cohort (https://portal.gdc.cancer.gov/projects/TCGA-CESC) using the gdc-client tool. A total of 257 expression files and 255 clinical files were available for analysis. Low-expressed genes (those with counts equal to zero in more than 25% of the samples) were filtered out to ensure robustness in our analysis. After filtering, data from 233 CSCC patients were included in the study. Using the np.corrcoef() function from the numpy module in Python, we calculated the correlation of different genes in the expression matrix. Additionally, we employed the stats.pearsonr() function from the Stats module to determine the significance of the calculated correlations. We utilized the survival (v.3.5.7) (https://github.com/therneau/survival) and survminer (v.0.4.9) (https://github.com/kassambara/survminer) R packages. Through the surv_cutpoint() function, we identified the optimal cut-points for the expression of the ‘APP’ gene. Subsequently, the surv_categorize() function was employed to categorize the gene based on these cut-points. Following this, the survfit() function was utilized to fit a survival analysis model and calculate survival curves for each group. Finally, the ggsurvplot() function was used for visualization, generating survival curves with accompanying p-values and a risk table. Protein interactions were analysed using STRING database (v12.0) (https://string-db.org [46]. The interactions with a combined score > 0.7 were selected to construct the protein-protein interaction (PPI) networks using the Cytoscape (v3.8.2).

Multiplex immunofluorescence assay

Detection was based on a 4-µm thick slide cutting from tumor tissues. After deparaffinisation and rehydration, the slides were subjected to epitope retrieval through boiling for 2 min in citric acid antigen repair buffer (pH6.0). Subsequently, endogenous peroxidase was blocked by incubation for 20 min in 3% hydrogen peroxide and later the protein was blocked in 10% goat serum albumin at 37 °C for 30 min. Then the four antigens were labelled by cyclic staining, including incubation of primary and secondary antibodies, tyramine signal amplification (TSA) visualization and removing the TSA–antibody complex in citric acid solution by microwave treatment for 5 min at 97 °C. In each round, antibody labelling was followed after epitope retrieval and protein blocking as mentioned above. After cycle staining, each slide was counterstained with DAPI for 5 min and images were scanned using a Pannoramic MIDI (3DHISTECH; Budapest, Hungary). Furthermore, primary antibodies are as follows: (1) CD74 (HUABIO; cat #HA601117; 1:1000), (2) APP (Abclonal; cat #A17911; 1:100), (3) TRPS1 (Abclonal; cat #A7743; 1:100), (4) KRT3 (Abclonal; cat #A0411; 1:100). For secondary antibodies, we utilized goat anti-rabbit/mouse secondary antibody labeled with HRP (KPL; cat #074-1506; cat # 074-1806) and incubated them for 1 h at 37 °C.

Immunohistochemical staining

The tumor sections (4 μm-thick) underwent deparaffinization, hydration, antigen retrieval, and were then sealed with 3% hydrogen peroxide. Subsequently, tissue sections were blocked with 3% bovine serum albumin for 30 min at room temperature, followed by overnight incubation at 4 °C with the primary antibody. The primary antibodies are as follows: (1) IDO1 (HUABIO; cat #HA721331; 1:200), (2) CD274/PD-L1 (Proteintech; cat #66248-1-Ig; 1:5000), (3) SPP1 (OriGene; cat #TA806784; 1:150). After washing with PBS, the slices were incubated with species-specific HRP labeled secondary antibodies for 50 min at room temperature. The secondary antibodies are goat anti-rabbit/mouse secondary antibody labeled with HRP (Pinuofei Biological; cat #PN0046; #PN0080). After incubation, the slices were washed by PBS and freshly configured DAB color developing solution was added. Then ddH₂O was used to terminate the color development. Harris hematoxylin was used to restain the nucleus. After dehydration and sealing with neutral resin, the slices were observed with the microscope and the images were collected and analyzed.

RNA isolation and quantitative reverse-transcription polymerase chain reaction

Total RNA was extracted from the HeLa cell lines using TRlzol reagent(Ambion, Austin, TX, USA), and reverse transcribed using NovoScript^®Plus All-in-one 1st Strand cDNA Synthesis SuperMix (gDNA Purge) (Novoprotein, Suzhou, China)according to the manufacturer’s protocols. The complementary DNA was then amplifed to detect APP, TRPS1, ITGA2, PRKCE, PKM, PSME4 and NF1 by qRT-PCR using a SYBR Green Master mix (Yeasen, Shanghai, China) according to the manufacturer’s protocols. All experiments were repeated independently three times. The mRNA expression was normalized to the expression of GAPDH mRNA and calculated using the 2-ΔΔCt method. The PCR primers used are listed in Supplementary Table S11.

CCK-8 assay

The proliferation rate of HeLa cells was analyzed using the Cell Counting Kit-8 (CCK-8, Vazyme, Nanjing, China) following the manufacturer’s protocols. The cells were seeded and cultured in medium containing 10% FBS in 96-well microplates. The working solution of CCK-8 reagent was prepared by adding 10 µL of the reagent to 90 µL of medium. Next, 100 µL of the working solution was added to each well and incubated for 45 min. We performed this assay at 0 h, 24 h and 48 h after transfection of siRNAs. All experiments were performed in triplicate. The absorbance at 450 nm was measured using a microplate reader (Bio-Rad, Hercules, CA, USA). Cell proliferation was determined based on the relative absorbance.

Transwell migration and invasion assay

The migration and invasion assays were performed using using an 8-mm pore size transwell system (Corning, NY, USA). For migration assay, HeLa cells were seeded into the upper chamber with serum-free medium (5 × 104 cells). For the invasion assay, the upper chamber was coated with Matrigel (Corning, NY, USA), and the subsequent steps were similar to the migration assay. After 12 h, cells in the upper chamber were transfected with siRNAs and medium with 30% FBS was added to the bottom of the chamber. After the cells migrated or invaded for 36 h, they were fixed with fixative solution (4% formaldehyde, methanol-free) (Biosharp, Anhui, China) and stained with 0.1% crystal violet (Macklin, Shanghai, China). For each insert, 5 fields were randomly selected to be captured under a microscope at 20×, and cell numbers were quantified using ImageJ2. The cell number of each treated group was normalized to the cell number of the control group to calculate the migration and invasion index.

Small interference RNA

Small interference RNAs (siRNA) against human APP and TRPS1 and the control siRNA were purchased from RIBOBIO (Guangzhou, China). Transfections were performed with Lipofectamine RNAiMAX reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s recommended protocol. HeLa cells were seeded to be 60-80% confluent at transfection. Lipofectamine RNAiMAX reagent were diluted in Opti-MEM medium (Gibco, Grand Island, NY, USA) and mixed with diluted siRNAs. After incubation for 10 min, siRNA-lipid complex was added to cells. The siRNAs were used at 60nM and cells were analyzed after 36 h.

Results

ST analysis reveals spatial expression patterns of cells in the LE and TC regions of CSCC

Through Stereo-seq ST sequencing, CSCC tissues and normal tissue (NOR) were analyzed (Supplementary Table S1). The normal cervical tissue was obtained from an area adjacent to the pelvic mass. The CSCC samples were obtained from the LE and TC regions, respectively, and these samples were HPV positive (Fig. 1A). We integrated published scRNA-seq data from cervical cancer patients to assist in the interpretation and optimization of the Stereo-seq data for further investigation [47]. The scRNA-seq data we analyzed include both CA and CSCC cases as references for cell annotation in ST data analysis (Fig. 1C). Despite their differences in cell composition, this approach aids in characterizing cell types and identifying the distribution of cell types within CSCC tissues (Supplementary Figure S1A and Table S3). Following quality control, we obtained data from a total of 15,318 cells and 29,663 genes (Fig. 1B and Supplementary Table S3). We identified 11 major cell types based on classic cell markers, including Fibroblasts (4,237), Cancer cells (3,705), Smooth muscle cells (420), Endothelial cells (179), T cells (2,676), Macrophages (1,210), B cells (943), NK cells (895), Neutrophils (860), Mast cells (127), and DC cells [66]. We found that cancer cells mainly expressed CDH1, EPCAM, CDKN2A, and TP63 in cervical cancer tissues (Fig. 1B, Supplementary Figure S1B, C and Supplementary Table S4). Endothelial cells showed high expression of CDH5, PECAM1 and VWF, while smooth muscle cells exhibited high expression of MYH11, MYLK, and RGS5. Fibroblasts displayed high expression of APOD and LAMA2 (Supplementary Figure S1B, C and Supplementary Table S4) [48]. Apart from structural cells, we also identified immune cells including T cells (CD3D+, CD3E+) and Macrophages (FCGR2A+, CD163+) [49]. These findings provided insights into the cellular composition of cervical cancer tissue, serving as a cell annotation reference for follow-up ST analysis.

Fig. 1

Spatial expression patterns of cells in the leading edge and tumor core regions of CSCC. (A) Overview of the schematic representation of the study design. Tumor tissues from CSCC patients were processed into single-cell suspensions, and unsorted cells were subjected to scRNA-seq using 10x Genomics. ST was performed on normal tissue and tumor slices using Stereo-seq, generating ST data. Created with BioRender. (B) UMAP plots of cells identified from scRNA-seq data of three CSCC tissues. Different colors represent distinct cell types. (C) The integrated analysis of single-cell transcriptome and spatial transcriptome data using cell2location is depicted. Created with BioRender. (D-F) Cellular spatial distribution and proportions of reads count for each cell type in normal tissue (D), leading edge region (E), and tumor core region (F)

To systematically localize the cell types identified by scRNA-seq in CSCC tissues, we utilized cell2location to integrate scRNA-seq with ST, mapping cells to CSCC tissues [37]. In the NOR, fibroblasts, smooth muscle cells, endothelial cells, and a small proportion of immune cells were observed (Fig. 1D). Fibroblasts, had the highest proportion among all the cell types in the LE. Particularly, immune cells, particularly macrophages, showed higher proportions in the LE, compared to those in the NOR and TC. In addition, those cells were mainly located in the distal region of the tumor (Fig. 1E). The TC region was mainly composed of cancer cells, with a few scattered endothelial cells and fibroblasts (Fig. 1F). These findings revealed the spatial distribution of different cell types within CSCC tissue, enhancing our understanding of the tumor microenvironment of cervical cancer.

Spatial transcriptomic variations and functional implications in LE and TC regions of CSCC

We conducted a differential gene analysis between the LE and TC, revealing differences in the expression of multiple immune-related genes. Notably, genes associated with macrophage function, such as NRP1, NRP2, MSR1, and SPP1, were significantly upregulated [50, 51] (Fig. 2A and Supplementary Table S5). Spatial transcriptomics data revealed a lower proportion of immune cells in the NOR and TC compared to the LE. In contrast, an aggregation of macrophages was observed in the area with fewer cancer cells in the LE region (Fig. 1D-F).

Fig. 2

Differential analysis of immune characteristics in different tumor regions. (A) Volcano plots (P-value and fold change) comparing gene expression between the leading edge and tumor core. Upregulated (orange) and downregulated genes (cyan) are highlighted, and differentially expressed immune-related genes (dark red) or metabolic genes (red) are labelled. Horizontal dashed line indicates adjusted P-value of 0.05, while the vertical dashed lines represent log2FC = -1 and 1. (B) Volcano plots (P-value and fold change) comparing gene expression of macrophages between the leading edge and tumor core(left) or normal tissue(right). Upregulated (orange) and downregulated genes (cyan) are highlighted. Horizontal dashed line indicates adjusted P-value of 0.05, while the vertical dashed lines represent log2FC = -1 and 1. (C) Macrophages spatial location(left) and spatial expression of SPP1 in leading edge region(right). (D) Immunohistochemical images depicting endogenous protein expression levels (brown) of SPP1 in leading edge region of CSCC sample (left) and normal sample (right), with black arrows highlighting macrophages. (E) Spatial expression of selected immune checkpoint genes (PD-L1 and IDO1) in normal tissue, leading edge region, and tumor core region. (F) Immunohistochemical images depicting endogenous protein expression levels (brown) of PD-L1 and IDO1 immune checkpoint genes in normal tissue, leading edge region, and tumor core region samples. (G) GSEA showing enrichment of Oxidative phosphorylation signalling pathwa in cancer cells from the leading edge and tumor core regions

In the LE, macrophages were aggregated in spots with very few cancer cells. These macrophages expressed SPP1 at significantly higher levels in the LE compared to the TC and NOR (Fig. 2B and Supplementary Table S6). Consistently, the expression of SPP1 was found to correlate with the distribution of these macrophages (Fig. 2C). Immunohistochemistry assays further confirmed higher expression of SPP1 in macrophages in the LE compared to the NOR region (Fig. 2D). Previous research has shown that SPP1 + macrophages exhibit immunosuppressive properties [52]. Additionally, these macrophages also exhibited high expression of CTSB (Fig. 2B and Supplementary Table S6), which is a prognostic biomarker for cancer progression and therapeutic target associated with immune cell infiltration and immunosuppression [53]. During cancer treatment, immune checkpoint genes play crucial roles in enabling cancer cells to evade immune system attacks, leading to tumor immune escape [54]. Our spatial analysis revealed substantially elevated expression of two immune checkpoint genes, PD-L1 and IDO1, in the TC compared to both the LE and NOR regions (Fig. 2E). The immunohistochemical staining of FFPE tissue samples further validated our spatial transcriptomics analysis, demonstrating higher expression levels of PD-L1 and IDO1 in the TC region (Fig. 2F).

In addition to the immune-related genes, a variety of metabolic genes showed significantly differential expression between the LE and TC (Fig. 2A and Supplementary Table S5). To investigate the functional implications of these metabolic gene changes and better understand the features of cancer cells in these two regions, we performed Gene Set Enrichment Analysis (GSEA) which revealed that oxidative phosphorylation in cancer cells was activated in the LE compared to TC (Fig. 2J), indicating spatial differences in energy requirements [55]. There might be metabolic differences between regions of the CSCC tumor architecture, which could significantly impact tumor growth and spread. In summary, the analysis of ST data revealed differences in immune features and metabolism among different tumor regions.

Dynamic resolution of tumor heterogeneous metabolic landscapes and cancer cell differentiation

Given that metabolic activity is a critical component of intra-tumor heterogeneity [56,57,58], we delved deeper into the metabolic features of different tumor regions and their correlation with the cancer cell differentiation status. To comprehensively evaluate metabolic activity in different regions of tumor tissue, we focused on six metabolic pathways: oxidative phosphorylation, glycolysis, hypoxia, lactic acid metabolism, lipid metabolism, and pentose phosphate [59, 60] (Supplementary Table S7). We calculated metabolic scores for these pathways in cancer cells and mapped onto tissue spatial locations (Fig. 3A, B). We observed that cancer cells in LE region showed a certain degree of metabolic activity, and their metabolic levels showed a relatively even distribution (Fig. 3A). In contrast, there were differences in metabolic levels among different areas in the TC region (Fig. 3B). These differences could be attributed to the diverse metabolic strategies adopted by cancer cells in response to varying microenvironmental conditions, ultimately impacting tumor progression and metastasis.

Fig. 3

Metabolic score and copy number variation analysis in different tumor regions of CSCC. (A-B) Spatial score of oxidative phosphorylation, hypoxia, glycolysis, lactate metabolism, lipid metabolism, and pentose phosphate pathway in the leading edge (A) and tumor core region (B) of Stereo-seq samples. (C-D) Grid division of the leading edge and tumor core regions, with yellow representing high metabolic scores and blue representing low metabolic scores. Line plot displaying the distribution of Median metabolic pathway scores across the five regions (C). The heatmap illustrating GSVA scores for various metabolic pathways in different grid regions (D). (E) Comparison of copy number variation (CNV) between cancer cells in different regions, including leading edge, hypermetabolic region in the tumor core region, and hypometabolic region in the tumor core region. The CNV status of each cell is displayed for chromosomes 1 to 22

To better understand the metabolic heterogeneity within the tumor, the LE and TC regions were divided into equal parts, which showed varying metabolic pathway scores. Regions classified as hypermetabolic had scores significantly above the overall median, whereas hypometabolic regions had scores significantly below it. Intermediate metabolic regions, defined as those within ± 10% of the overall median, represented areas with metabolic activity close to the median value. These classifications provide a nuanced view of the metabolic landscape within the tumor, highlighting areas of high, low, and intermediate metabolic activity (Fig. 3C). Consistent with the method of ranking the metabolic activity based on the median score of six representative pathways, the results from GSVA demonstrated similar trends in metabolic activity within the tumor regions (Fig. 3D). Interestingly, all regions in the LE were hypermetabolic regions, whereas there was a mix of hypermetabolic and hypometabolic spots in the TC, revealing the complex and dynamic nature of tumor metabolism.

To assess copy number variation (CNV) of cancer cells across different regions, we compared chromosomal CNV between the LE and TC regions. Hypermetabolic regions Leading edge1 and Leading edge2 from the same sample, were combined for a more comprehensive analysis. This consolidation allowed for a more robust comparison of CNV profiles between the LE and TC regions. Unsupervised hierarchical cluster analysis of CNV values revealed that hypermetabolic regions were clustered together, suggesting their similarity and monoclonal origin. Notable differences in CNV patterns of cancer cells between hypermetabolic and hypometabolic regions were observed, suggesting potential associations with metabolic pathway activation (Fig. 3E). Examination of the expression of metabolic genes in different tumor regions showed a significant upregulation of many metabolism-related genes, including PRKCE, GNGT1, and ABHD2, in cancer cells within hypermetabolic regions (Fig. 4A). Additionally, genes such as MDM4, EPAS1, and ATM also showed higher expression levels in endothelial cells and fibroblasts in the hypermetabolic regions (Fig. 4A). This finding suggests that not only the metabolic activity of cancer cells themselves but also surrounding cell types might be involved in regulating the metabolic process within tumor tissues.

Fig. 4

Spatial and functional characteristics of different metabolic tumor regions in CSCC. (A) Heatmap displaying the expression profiles of 50 metabolic genes in different tumor regions. Colors represent the expression levels of the corresponding genes; the “blue-white-red” pattern maps to the minimum-average-maximum expression values. (B) Heatmap showing transcription factors and their expression levels in different tumor regions. Color intensity represents z-score values, and the size of the dots indicates the regulon specificity score (RSS) of specific transcription factors in different tumor regions. (C) Spatial expression of the TRPS1 in cancer cells of leading edge and tumor core region. (D) PPI network graph depicting the interaction and enrichment analysis results of TRPS1 target genes. Each pathway is presented as circular nodes in the network graph, with node size proportional to the number of input genes, and node color indicating the cluster it belongs to. Pathways with a similarity score greater than 0.3 are connected by edges, and edge thickness represents the degree of similarity between them. (E) RNA velocity streams of leading edge and tumor core region, colored by cell annotations for leading edge and tumor core regions

Transcriptional regulatory networks were distinct across regions at different levels of metabolic activities (Fig. 4B). The analysis showed that the transcription factor TRPS1 was highly expressed in the TC region, which was previously demonstrated involved in the occurrence and development of breast cancer through aberrant transcriptional regulation of oncogenes [31]. Mapping the expression of TRPS1 onto spatial locations revealed higher expression levels in cancer cells within hypermetabolic regions, with relatively high expression in the LE region as well (Fig. 4C). Enrichment analysis of target genes for TRPS1 revealed they were mainly enriched in cell-cell junction, cell adhesion molecule binding, and cell morphogenesis pathways (Fig. 4D), suggesting increased proliferation and differentiation potential in cancer cells within hypermetabolic regions [61, 62]. ZEB1, a transcription factor, is recognized as an important regulator for tumor metastasis and invasion [63]. It exhibited markedly higher expression in the LE compared to the TC (Fig. 4B). Mapping the expression of ZEB1 onto spatial locations revealed its abundant expression in both LE and TC regions, with its target genes most highly enriched in the GMP metabolic process (Supplementary Figure S2A, B). Furthermore, RNA Velocity analysis indicated a differentiation trajectory of cancer cells from hypermetabolic to hypometabolic regions in both LE and TC regions (Fig. 4E). We also found that the cancer stem cell-associated gene CD44 was highly expressed in hypermetabolic regions, while its expression was reduced in hypometabolic regions (Supplementary Figure S2C) [64,65,66]. This indicated a close correlation between metabolism and cell differentiation, suggesting that important role of metabolism in regulating cell differentiation. These findings improve our understanding of intra-tumor heterogeneity by revealing a distinct transcriptional regulatory profile across hypo- and hyper-metabolic regions.

The tumor microenvironment varies across different metabolic regions

Cancer cells were surrounded by clustered fibroblasts and endothelial cells within the hypermetabolic spots (Fig. 5A). Further analysis of cellular composition proportions revealed that in the TC, over 80% of the cells are cancer cells. However, the proportions of other cell types, such as endothelial cells and fibroblasts, remained higher in the hypermetabolic region (Tumor core1) compared to the hypometabolic regions (Tumor core2 and Tumor core3) (Fig. 5A). These findings suggests that the hypermetabolic state of cancer cells may be regulated by surrounding cells, possibly through cell-cell communication. We thus assessed the receptor-ligan interactions between cancer cells and other cells within the tumor microenvironment. Multiple pairs including SEMA3C-NRP2 and Dehydroepiandrosterone-bySTS-ESR1 were interacted between cancer cells and fibroblasts. Previous studies have shown that the binding of SEMA3C to NRP2 can promote angiogenesis [67], and SEMA3C can regulate tumor progression [68]. In addition, we found that communication between cancer cells and DC cells, endothelial cells, and macrophages primarily relies on the APP-CD74 signalling pathway (Fig. 5B). Differential gene expression analysis of macrophages, DC cells, and endothelial cells between hypermetabolic and hypometabolic regions revealed distinct differences in the expression patterns (Figure S3A). For instance, in hypermetabolic regions, endothelial cells show increased mRNA expression levels of ZEB1 (Supplementary Figure S3A and Supplementary Table S8-S10), suggesting that metabolism regulates the tumor microenvironment through the EMT-inducing transcription factor ZEB1. Enrichment analysis showed varied signalling pathways of these cell types between hyper- and hypo-metabolic regions (Figure S3B). Notably, we observed that some immune response-related pathways were activated in DC cells within the hypermetabolic region, including ‘positive regulation of phagocytosis’.

Fig. 5

Spatial and functional characteristics of different metabolic tumor regions in CSCC. (A) Spatial distribution of cancer cells (red dashed line indicates hypermetabolic region, green dashed line indicates hypometabolic region) (left), and other cell types in the different regions of the leading edge and tumor core (middle). Proportion of each cell type in each region (right). (B) Heatmap of cell communication analysis between cancer cells and other cell types. Dot size reflects the significance level of cell communication, and color indicates the average expression level of interacting molecules between cancer cell clusters and clusters of other cell types. (C) Correlation analysis between the expression level of APP and the expression level of genes related to tumor metabolism (ITGA2 and PRKCE) in different cell types. (D) Correlation analysis between the expression level of APP and the tumor metabolic score in different gridded regions. (E) Correlation analysis between the expression level of APP and the expression level of TRPS1 in different cell types. (F-G) Spatial distribution of gene expression of APP (F), PRKCE and ITGA2 (G) in cancer cells. (H) Immunofluorescence microscopy observation of cervical cancer slices, with DAPI (blue), KRT (purple), CD74 (green), APP (red), and TRPS1 (orange) antibody staining

APP-CD74 has been implicated in cellular communication within the tumor microenvironment in several cancer-related studies [69]. However, previous research has not elucidated its specific functional relevance to tumor progression [70]. By confirming the importance of APP-CD74 in tumor progression, we aim to address the current knowledge gap regarding its role in cancer biology. The correlation analysis results suggest that the expression of APP within the tumor microenvironment is positively associated with the expression levels of various metabolic genes, such as PRKCE, ITGA2, PKM, PSME4, and NF1, genes associated with cancer progression [71,72,73,74,75] (Fig. 5C and Supplementary Figure S4A). Moreover, the expression levels of APP and TRPS1 were significantly positively correlated in the hypermetabolic regions (Fig. 5E). However, there is no significant correlation between the expression of CD74 and metabolic genes (Supplementary Figure S4B). To further investigate the impact of the tumor microenvironment on cellular metabolism, our focus will be on exploring the relationship between APP and metabolic-related genes. Upon closer examination, a positive correlation was found between the metabolic scores of cancer cells in five gridded regions and the expression of APP (Fig. 5D). We observed higher levels of APP expression in cancer cells located in regions with hypermetabolic activity (Fig. 5F). Interestingly, other metabolic-related genes like PRKCE, ITGA2, PKM, PSME4, and NF1 showed a similar spatial distribution pattern (Fig. 5G and Supplementary Figure S4D). These findings suggest a potential link between the expression levels of APP and multiple metabolic genes, indicating that APP might play a role in influencing tumor metabolic activity. Immunofluorescence staining of human CSCC tissues revealed that CD74 protein was aggregated around cancer cells stained with KRT3. These cells exhibited abundant expression of APP, which was co-localized with TRPS1 (Fig. 5H). These results suggest that CD74 is expressed in the extracellular matrix surrounding the cancer cells, while APP is prominently expressed within the cancer cells, co-localizing with the transcription factor TRPS1. This indicates a potential reciprocal regulatory relationship between APP and TRPS1. The results of correlation analysis and immunofluorescence experiments suggest that APP likely plays a crucial role in regulating the metabolic activity of cancer cells. Moreover, its mechanism of action may involve metabolic-related genes and the transcription factor TRPS1.

APP and TRPS1 are associated with the progression of CSCC

Our research has identified the overlap of high expression of APP regions and hypermetabolic regions where the transcription factor TRPS1 is also expressed (Figs. 4C and 5F). Prior studies have shown that the cleavage product of APP, known as Aβ, has the potential to impact mitochondrial oxidative stress and energy metabolism [76]. Furthermore, previous research has found that abnormal expression of TRPS1 could drive the evolution of the tumor genome, potentially enhancing the adaptability of cancer cells [31]. Based on our previous analysis, APP and TRPS1 both are closely associated with the metabolism of CSCC cancer cells. Therefore, we further explored their functional relationship and their impact on the metabolism and progression of CSCC. By analysing transcriptomic expression profiles of 233 CSCC patients from TCGA database, we found that not only is the expression of APP positively correlated with TRPS1, but it is also positively correlated with the expression levels of multiple metabolic-related genes (such as ITGA2, PRKCE, PKM, PSME4, and NF1) (Fig. 6A and Supplementary Figure S4C). These findings indicate that APP and TRPS1 may play a key role in regulating the metabolic status of CSCC.

Fig. 6

Spatial and functional characteristics of different metabolic tumor regions in CSCC. (A) Correlation analysis between the expression levels of APP in RNA-seq data of TCGA CSCC tissues and the expression levels of the transcription factor TRPS1 and genes related to tumor metabolism (ITGA2 and PRKCE). (B) Bar graph showing siRNA-mediated knockdown of APP (left) and TRPS1 (right) RNA expression in HeLa cells. (C) Bar graph showing TRPS1 RNA expression in control and knockdown of APP HeLa cells (left) and APP RNA expression in control and knockdown of TRPS1 HeLa cells. (D) Bar graph showing the RNA expression of tumor metabolism-related genes (ITGA2 and PRKCE) in control and knockdown of APP HeLa cells. (E) Cell proliferation measured by CCK8 assay in HeLa cells transfected with si-APP, si-TRPS1, or si-NC. (F-G) The effect of APP or TRPS1 knockdown on HeLa cells’ migration (F) and invasion (G) determined by transwell assays. (H-I) Estimation of the survival probability of CSCC patients based on APP (H) or TRPS1 (I) gene expression levels. The analysis was performed using TCGA’s dataset of 255 CSCC samples. (J) Schematic summary illustrating immune characteristics and associations of APP and TRPS1 with CSCC progression, generated using BioRender (https://biorender.com/). Significant differences were determined using the two-tailed unpaired Student’s t-test (****p < 0.0001; ***p < 0.001; **p < 0.01; *p < 0.05; ns, not significant). P values < 0.05 were considered statistically significant

Cancer cell metabolism abnormalities are closely related to cell proliferation [7]. To explore the effects of APP and TRPS1 on cancer cell metabolism and proliferation in the development of CSCC, we first performed PPI analysis, which confirmed no direct interaction between APP and TRPS1 (Supplementary Figure S4E). Next, we individually knocked down APP and TRPS1 in the HeLa cell line (Fig. 6B). Initially, we observed a decrease in TRPS1 expression when APP expression was reduced. However, knocking down TRPS1 did not lead to a significant alternation in APP expression (Fig. 6C), indicating that APP may exert its function partially through the regulation of TRPS1. Additionally, we evaluated the expression of metabolic-related genes after knocking down APP and observed a significant decrease in mRNA expression levels of genes such as ITGA2 and PRKCE (Fig. 6D and Supplementary Figure S4F), consistent with the transcriptomic analysis of CSCC. Subsequently, we performed CCK8 assays and Transwell experiments to analyse the impact of APP and TRPS1 on cancer cell proliferation, migration, and invasion.

CCK8 assays showed that knocking down either APP or TRPS1 markedly attenuated the proliferation of HeLa cells within 48 h (Fig. 6E). The migration and invasion of cells were significantly inhibited in response to either APP or TRPS1 knockdown (Fig. 6F). These results demonstrate the important functional roles of APP and TRPS1 in regulating cervical cancer cell growth, further supporting our ST data analysis findings. To validate the pro-tumor effects of APP and TRPS1, we also conducted survival analysis using the data from TCGA CSCC cohort, including 233 patients. Higher expression of APP was found to be significantly correlated with poorer survival outcomes (p = 0.003, Fig. 6H). Similarly, higher levels of TRPS1 expression were also associated with lower survival rates among CSCC patients (Fig. 6I, p = 0.034).

Overall, our experimental results have uncovered the regulatory roles of APP and TRPS1 on multiple metabolic-related genes in cancer cell metabolism, aligning with the survival analysis using the genomics data from the TCGA CSCC cohort. Furthermore, both transcriptional regulators are crucial in cancer cell proliferation, migration, and invasion. Although the precise regulatory interplay between APP and TRPS1 warrants further investigation, our study provides novel functional evidence, laying a solid foundation and identifying these factors as promising targets for CSCC treatment.

Discussion

The relationship between changes in tumor metabolism and the occurrence and development of tumors entails intricate biological processes and multiple molecular mechanisms [7]. Cancer cells can adjust their own energy sensing pathways to enhance adaptation to metabolic stress, which increases their survival rate under low-nutrient conditions and thus contributes to tumor drug resistance [77]. Currently, there are several studies on the metabolism and tumor progression of CSCC, but most of these studies lack a deep understanding of the interactions between different cell types and the impact of their spatial positions on the tumor microenvironment [78, 79]. Although ST technology and scRNA-seq analysis of cellular spatial distribution of CSCC have contributed to understanding of the tumor microenvironment, comprehensive knowledge of its metabolic characteristics is still lacking [23, 24, 80]. This study utilized Stereo-seq ST technology and scRNA-seq datasets to explore the heterogeneity of metabolic regions within CSCC tissues and correlate them with cancer cell differentiation status. The analysis of the tumor microenvironment unveiled the crucial roles of APP and TRPS1 in regulating CSCC cancer cell metabolism and proliferation that are associated with patient prognosis (Fig. 6J). These findings provide a new theoretical foundation and potential targets for the treatment and personalized therapy of CSCC.

Various types of cells and microenvironments exist within tumor tissues. Typically, the LE region refers to the part adjacent to normal tissues. It may be significantly influenced by the tumor microenvironment and immune cell infiltration, showing signs of metastasis, infiltration, or transformation [25]. The TC region, on the other hand, represents the central part of the tumor where cancer cells have higher genomic instability, making them more prone to metastasis and dissemination [81]. We utilized stereo-seq and scRNA-seq to further analyze the spatial molecular characteristics of the LE and TC regions of CSCC by deconvolution using the cell2location method [37]. We observed varying distributions of immune cells and metabolic activities within the CSCC tumor microenvironment across the LE and TC regions. Specifically, in the LE region, there was an accumulation of macrophages expressing high levels of SPP1. Furthermore, previous research has also demonstrated that the changes in macrophage polarity defined by the expression of the CXCL9 and SPP1 genes are a crucial characteristic of the tumor microenvironment [82]. Macrophages in the LE region, characterized by high SPP1 and CTSB expression, may have immunosuppressive properties. Further investigation is needed to understand their role in immune evasion and their potential impact on treatment strategies and patient prognosis. In the LE region of CSCC tissues, apart from the accumulation of SPP1 + macrophages, there are indications of heightened inflammation-related pathways, reflecting the complex immune regulatory mechanisms within the tumor microenvironment [83, 84]. The TC region exhibits a higher immune-suppressive environment, characterized by elevated levels of PD-L1 and IDO1 expression, along with lower activity in inflammation-related pathways. These findings further underscore the crucial roles of immune checkpoint genes and inflammatory pathways in tumor development [54, 85]. In terms of metabolism, the LE region displays elevated metabolic activity, such as cGMP metabolism, oxidative phosphorylation, and purine metabolism. These findings not only highlight the heterogeneity within CSCC tumor but also provide crucial clues for the discovery of personalized therapeutic targets. In line with our research findings, recent studies have also underscored the heterogeneity of various regions within tumors. For instance, a recent study identified distinct transcriptional characteristics in the TC and LE regions of oral squamous cell carcinoma, corroborating our observations [86].

Further delving into the metabolic characteristics within the CSCC tumor, we found that distinct metabolic signatures exist even within the LE and TC regions. Additionally, we observed differential patterns of CNV in high and low metabolic areas of cancer cells, which may be associated with alterations in metabolic-related genes, suggesting a relationship between metabolic status and patterns of genetic variation. John et al.‘s study clarified that p53 promotes the differentiation of cancer cells by regulating cellular metabolic pathways, particularly the accumulation of α-ketoglutarate [87]. In our research, we found a close correlation between metabolic activity and the differentiation status of cancer cells, suggesting that cancer cells may differentiate from regions of hypermetabolic activity to those of hypometabolic activity. These findings are consistent with previous studies showing that cancer cells often transition from a hypermetabolic, stem cell-like state to a hypometabolic, differentiated state [7, 88]. This is observed in cancers such as breast cancer and glioblastoma, where high metabolic activity not only sustains stem cell properties but also enables tumor heterogeneity and promotes tumor growth [89,90,91]. Additionally, our research has revealed close connections between metabolic activity, cancer cell differentiation, and transcriptional regulatory networks. Transcription factor analysis revealed higher expression of TRPS1 and ZEB1 in hypermetabolic regions. TRPS1 and ZEB1 have been extensively studied in tumor development and progression [92,93,94,95]. Our research found that their transcriptional regulation of cancer cells may be closely associated with tumor metabolism.

We have discovered specific cell-cell communication between cancer cells and other cell types, such as SEMA3C-NRP2 and APP-CD74, which may play a role in tumor growth and metastasis [68]. In Bian et al.‘s study, it was found that the heightened uptake of methionine by cancer cells reduces methionine availability for cytotoxic T cells, resulting in diminished cytotoxic T cell function and suppression of effector T cells [96]. Additionally, in our study, we also found that in different regions of varying metabolic activity in cancer cells, other cell types exhibit different transcriptional patterns, reflecting how cancer cells can regulate the status and function of immune and stromal cells surrounding them through metabolism. Furthermore, we have identified key regulatory factors APP and TRPS1, which play crucial roles in regulating the metabolism, proliferation, migration, and invasion of CSCC cancer cells. Kleffman et al. found that the secretion of APP by melanoma cells promotes cancer cell metastasis and activates an anti-inflammatory phenotype [97]. Yang et al.‘s study on breast cancer revealed the critical role of TRPS1 in driving heterochromatin activation and cancer cell genomic evolution [31]. These conclusions corroborate our research findings and complement each other, elucidating their significant roles in the onset and progression of cancer. Targeting these factors therapeutically could present new avenues for treatment. For instance, the development of APP or TRPS1 inhibitors could potentially disrupt the pro-tumorigenic signalling pathways they influence [98,99,100,101].

In our study, there are some limitations that need to be considered. Firstly, our sample size is relatively small. Although we employed rigorous statistical methods during the analysis to ensure the reliability of the results, the small sample size may still lead to some degree of bias in the outcomes. However, we have addressed this limitation by validating our findings through experimental results. Through immunohistochemistry and immunofluorescence experiments, we have obtained results consistent with those from the ST data analysis. Through in vitro experimentation, we conducted functional validation of APP and TRPS1 by manipulating their expression in the HeLa cell line. This indicates that even with a limited sample size, our research findings still hold a certain degree of credibility and significance. Additionally, the spatial resolution of ST may limit the precise localization of cell subpopulations, thereby affecting the interpretation of results. Furthermore, sample preservation and processing may impact the transcriptional characteristics of the cells, particularly regarding RNA degradation and expression changes. Finally, future studies should consider integrating other technologies and improved algorithms to overcome these limitations and provide a more comprehensive perspective on metabolic analysis [102,103,104,105,106,107]. To improve the reliability of future studies, larger sample sizes, especially from multicenter studies that include diverse populations, should be considered to validate our findings.

Conclusions

Overall, our research findings bridge the gap in knowledge regarding tumor metabolism and CSCC development. By elucidating the pivotal roles of APP and TRPS1 in regulating tumor metabolism and proliferation, we offer new theoretical foundations and potential targets for CSCC treatment. Moreover, our study underscores the significance of recognizing metabolic heterogeneity within tumors and its implications for cancer treatment. In conclusion, our discoveries lay the groundwork for the advancement of novel therapeutic approaches targeting CSCC tumor metabolism, offering promising avenues for enhancing patient prognosis and personalized treatment modalities.