Integrated spatial multi-omics delineates fatty acid degradation fuels malignant evolution at the tumour periphery in cervical squamous cell carcinoma.

The spatial transcriptomics data for normal cervical and CSCC samples have been deposited in the Genome Sequence Archive (GSA) of the National Genomics Data Centre (Access Link: https://ngdc.cncb.ac.cn/gsa-human/browse/ ; ID: HRA006531 ). Additionally, the metabolic data for these samples are available in the Metaspace database (Access Link: https://metaspace2020.eu/api_auth/review?prj=baf8d214-b6a0-11ee-adab-673ac9ff27e7&token=0q4-Q9kPOMpK ).

This study was conducted in compliance with the Declaration of Helsinki and Good Clinical Practice guidelines, with approval granted by the Medical Ethics Committee of Tongji Medical College, Huazhong University of Science and Technology (approval number: TJ-IRB20221247). Samples for ST and spatial-MSI were collected from six treatment-naïve patients diagnosed with CSCC, and two patients with benign uterine adenomyosis. All patients were of East Asian (Han Chinese) ethnicity, and their detailed clinical characteristics were provided in Table S1 . After isolating the tumour tissue, the specimen was immediately placed on ice, and embedded in optimal cutting temperature compound (OCT, Sakura, USA) at −30 °C. Subsequently, it was stored at −80 °C for further sequencing. The sequencing of ST was conducted by Ouyi Biomedical Technology Co., Ltd (Shanghai, China) following the protocols of 10x Genomics Visium. 21 Briefly, cervical tissues embedded in OCT were cryo-sectioned into 10 μm-thick slices. These slices were then mounted on an ST microarray, obtained from the Spatial Transcriptome Team ( https://www.10xgenomics.com/ ). Subsequently, the sections were immersed in isopropanol for 1 min and stained with haematoxylin and eosin (H&E) to generate bright-field images. Following staining, the tissues underwent permeabilization, and cellular mRNA was captured by primers at the gene delivery site. All cDNA synthesised from the captured mRNA contained a uniform spatial barcode. A sequencing library was constructed from this cDNA and sequenced using the Illumina NovaSeq 6000. Ultimately, the spatial barcode enabled the alignment of sequencing data with the corresponding tissue slice image, facilitating the construction of a spatial gene expression map. Adjacent tissue sections used for ST were subjected to metabolic sequencing via AFADESI-MSI analysis, performed by Luming Biotech (Shanghai, China), following the manufacturer's instructions. 21 In brief, 10 μm-thick tissue slices were mounted on a positively charged desorption plate (Thermo Scientific, U.S.A) and dried sequentially at −20 °C for 1 h, followed by room temperature for 2 h prior to MSI analysis. The AFADESI-MSI system from Beijing Victor Technology Co., LTD (Beijing, China) was interfaced with a Q-Orbitrap mass spectrometer (Q Exactive, Thermo Scientific, U.S.A). For analysis, the negative ion mode utilised acetonitrile (ACN)/H2O (8:2) as the spray solvent, while the positive ion mode utilised ACN/H2O (8:2, 0.1% FA). The solvent flow rate was maintained at 5 μL/min, with a transfer gas flow rate of 45 L/min. The spray voltage was set to 7 kV, with 3 mm distance between the sample surface and the spray, as well as between the spray and the ion transfer tube. Mass spectrometry parameters included a resolution of 60,000, a mass range of 70–1000 Da, an automatic gain control (AGC) target of 2E6, the maximum injection time of 200 ms, an S-lens voltage of 55 V, and a capillary temperature of 350 °C. MSI scanning was conducted at a speed of 0.2 mm/s in the X direction and a step size of 100 μm in the Y direction. Fastq reads were aligned to the human reference genome GRCh38 using SpaceRanger (10x Genomics). Gene expression matrices, spatial barcodes, coordinates, and H&E images were imported into Seurat (v4.3.0). For integration of eight 10x Visium samples, we employed Seurat's standard workflow to mitigate batch effects. Individual samples were normalised using SCTransform, followed by identification of 3000 highly variable features. Integration anchors were established via FindIntegrationAnchors, and datasets were harmonised using IntegrateData, yielding a batch-corrected combined matrix. PCA was performed on the integrated data. A shared nearest-neighbour (SNN) graph was constructed from the first 30 principal components (PCs) using FindNeighbors, and clustering executed via FindClusters at resolutions spanning 0.1–1. A resolution of 0.5 optimally delineated distinct spot clusters. Uniform manifold approximation and projection (UMAP) was subsequently applied to the top 30 PCs for visualisation. 15 publicly available Stereo-seq datasets (CNP0002543) were processed identically in Seurat for normalisation and clustering, though without integration. Dimensionality reduction and clustering parameters matched those of the 10x Visium analysis. scRNA-seq data from Zhang et al. served as the reference for deconvolution using Robust Cell Type Decomposition (RCTD), a method robust to reference-technique mismatches. 22 , 23 This reference was annotated into seven major lineages: epithelial cells ( EPCAM , CDH1 , KRT5 , TP63 ); myeloid cells ( ITGAX , CSF1R , CD14 , FCGR3A ); endothelial cells ( CLDN5 , VWF , CDH5 , KDR ); fibroblasts ( COL1A1 , COL1A2 , LUM ); B cells ( MZB1 , CD79A , MS4A1 ); T cells ( CD3D , CD3E , CD2 ); mast cells ( KIT , IL1RL1 , MS4A2 ). All 10x Visium data were deconvolved using run. RCTD (RCTD package; doublet_model = full), quantifying per-spot cell type proportions. 23 Cluster-specific cell type abundances were compared via Wilcoxon tests, with median proportions informing hierarchical clustering. Epithelial-enriched clusters (tumour/normal) were classified based on epithelial cell dominance. Stereo-seq datasets were analogously deconvolved. Epithelial and non-epithelial clusters were manually annotated per sample using cell type abundances and H&E histology. Spatial continuity was computed as described using the following formula: 7 Spatial continuity = ∑ i ∑ j ∈ Neigh ( i ) 1 ( cluster j = Tumour ) ∑ i ∑ j ∈ Neigh ( i ) 1 where i indicated a spot from tumour cluster and Neigh(i) was the seven spots including the surrounding six neighbour spots and the spot i itself. By this method, the higher the spatial continuity, the more spatially concentrated the tumour spots are. To characterise relationships between epidermal gene functions and tumour purity (epithelial percentage), we applied segmented regression (segmented R package). 24 Linear models initially regressed gene-set scores against epithelial percentages per spot. Breakpoints in these relationships were identified using the segmented function, with fitted curves visualised via ggplot2. We utilised inferCNV to assess copy number variant in six tumour samples. 25 Initially, we extracted normal epithelial clusters from samples N1 and N5, which were then used as reference in inferCNV, individually (N1 epithelial clusters, N5 epithelial clusters) and in combination. The copy number profiles exhibited high similarity across all three reference groups. Based on these findings, we proceeded with the combined reference cluster approach for the final analysis. Tumour layers were delineated using an advanced clustering algorithm in Cottrazm, which employed ST and H&E stained images to segregate spots into distinct clusters. 26 For subsequent inter-layer comparisons, however, we utilised the SCT assay generated by Seurat instead of relying on morphologically adjusted gene expression. Differentially expressed genes (DEGs) between layers were identified using the FindAllMarkers function in Seurat, with parameters set as |Log 2 FC| > 0.25, P -value <0.05. These DEGs played a pivotal role in enriching Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, cancer state functions, and cancer hallmarks. 27 Pathway activity in CancerSEA and hallmarks was evaluated using the AddModuleScore function in Seurat. Additionally, cervical intraepithelial neoplasia (CIN) signatures from Liu et al., ( Table S2 ) were leveraged to investigate the initiation of CSCC. 28 We utilised Monocle 2 to conduct trajectory analysis on tumour focal layers. 29 Initially, clusters from the tumour layers in samples SM27 and SM35 were isolated, retaining only genes expressed in at least five spots. Subsequent analysis involved calculating DEGs between layers using the differentialGeneTest function. The top 1500 significant genes, ranked by q-value, were selected for spot ordering. Dimensionality reduction was carried out using the DDRTree method, and trajectory visualisation was generated with the plot_cell_trajectory function. Pseudo-time values for each spot were integrated into the Seurat object. Ultimately, the spatial distribution of pseudo-time was visualised using Seurat's Spatial plot. To investigate the spatial characteristics of distinct tumour layers, the tumour boundary was delineated for each tumour focus. Initially, spots were selected across all layer, with each spot surrounded by six adjacent spots on the spatial slide. If all surrounding spots resided within the tumour focus, the spot was classified as an inner tumour spot. Conversely, if fewer than six surrounding spots were within the tumour focal, the spot was designated as a boundary spot. Next, the boundary distance of focal spots was computed based on the shortest path to boundary spots. Finally, the distances of spots from each layer were aggregated. In the cell-cycle analysis, we employed Seurat to assign the cell-cycle scores to spatial transcriptomic data. Each cell was classified into a specific cycle phase (G1, S, and G2M) based on established marker genes. Subsequently, cells were grouped into layers according to their spatial clusters. We visualised the cellular trajectory, colour-coded by cell-cycle phase, and mapped the spatial distribution of these phases across the tissue layers. We conducted a TF regulatory analysis using the R package SCENIC. 30 Briefly, the gene expression matrix for selected spatial spots was estimated using GENIE3. Subsequently, the RcisTarget package was utilised to identify TF motifs. The regulon activity score for each spot was then computed using the AUCell package. Finally, regulon activities exceeding 0.1 were selected for visualisation. We utilised mistyR to estimate the dependence coefficient between cell types and functional pathways. 31 Different models were employed in the mistyR analysis. The intrinsic model assessed cell type co-localisation based on the percentage of each cell type within individual spots. In the para view, we evaluated the impact of neighbouring cell type proportions on hallmark pathways of tumour foci across all spots. This para view elucidated the spatial dependence of tumour development mechanisms on the abundance of surrounding cell types. To determine whether spatial stratification also exists in additional specimens, we utilised external data from Stereo-seq to evaluate the estimated scores of spatial layers and certain cancer hallmarks across distinct tumour foci. 10 Two samples (TJH 37/90) with obvious tumour boundaries were selected for further analysis. Using the FindAllMarkers function, we identified gene sets associated with different layers in samples SM27 and SM35 ( Table S3 ). Subsequently, we employed the AddModuleScore function with default parameters to score the ST data. Tumour foci were then isolated using the subset function, partitioned into distinct regions, and subjected to detailed visualisation to more precisely characterise their spatial features. Raw mass spectrometry data were converted into imzML format using imzMLConverter software. 32 Subsequent preprocessing steps, including data normalisation, spectral smoothing, baseline reduction, peak picking, alignment and filtering, were performed with the Cardinal tool. 33 Positive and negative ions in the spatial metabolomics dataset were identified by processing imzML files using the pySM framework. 34 Specifically, metabolite annotation was performed against two reference databases: the in-house library SmetDB (developed by Shanghai Luming Biotechnology Co., Ltd.) and the public Human Metabolome Database (HMDB) ( www.hmdb.ca ). Metabolite extraction from the imaging data (resolution: 60,000) was conducted using Metabolite Signal Match (MSM) scoring, with a mass error tolerance of 5 ppm and isocalc_sig set to 0.01. To eliminate false positives, a decoy adduct strategy was employed (fdr.n_im = 20), with the target false discovery rate (FDR) set at 0.1 (10%) and 20 replicates for threshold calibration. Only annotations with MSM scores that passed FDR validation were retained for subsequent analysis. Multivariate statistical analysis was performed to cluster the samples using orthogonal partial least squares discriminant analysis (OPLS-DA). Differential metabolites between tumour and normal samples were identified based on ion intensity, which was randomly sampled at the same times from two group samples. A T-test with a P- value <0.05 were used to achieve the significant differential metabolites. In the regional comparison of metabolite intensity, we firstly manually conducted image fusion and spatial matching between the MS data and the H&E image through the software MsiReader. Then, we selected regions of interest according to the H&E image, the location of which would be mapped to the MS data. The metabolic profiles in different regions were used in OPLS-DA. The variable importance in projection (VIP) value was calculated for each ion to indicate their contribution to the OPLS-DA model. Besides, the Student's T-test was also used in this comparison, significant metabolites between regions were with a VIP value > 1.0 and P -value <0.05 in independent T-tests. To investigate the dynamic intensity of metabolites in the tumour focal foci, we selected the tumour focal foci from samples SM27 and SM35 using MsiReader. Within the focal foci, boundary spots were labelled as distance 0. The distance of each spot from the boundary was then determined by calculating the minimal Euclidean distance between the queried spot and all boundary spots. Pearson correlation analysis was performed to assess the relationships between metabolite intensity and distance from the boundary across all spots in the focal foci. Metabolites with a P- value <0.05 were classified as dynamic, indicating significant variation from the local centre to the boundary. The detailed clinical data of patients with CC, including gene expression profiles, survival status, and survival times, were retrieved from the TCGA database ( https://www.cancer.gov/tcga ) using the TCGAbiolinks R packages, as summarised in Table S4 . The TJH cohorts’ clinical characteristics were comprehensively outlined in Table S5 . The layer 3- specific gene signature was generated by intersecting the marker genes of cluster 0 (sample SM27) and cluster 9 (sample SM35), with an adjusted P -value 0.6 as the screening thresholds. For each patient in the TCGA cohort, the signature score was subsequently computed using gene set variation analysis (GSVA). Patients were stratified into high- and low-score groups based on the median estimated score of the layer 3 signature or the targeted protein staining, and univariate Cox regression analysis was employed to evaluate the prognostic significance of specific factors in the CC cohorts. We identified DEGs using the Wilcoxon rank-sum test implemented in the “FindAllMarkers” function of Seurat (Version 4.3.0). DEGs were screened based on the following thresholds: |Log 2 FC| > 0.25 and adjusted P- adj <0.05. Gene Ontology (GO) enrichment analysis for DEGs was conducted using the DAVID database ( https://david.ncifcrf.gov/home.jsp ). Gene set signature scores were computed for each spot utilising the “AddModuleScore” function in Seurat (Version 4.3.0). Additional details on the gene set signatures were displayed in Table S2 . The CC cell lines (HeLa and SiHa) were obtained from the American Type Culture Collection (ATCC), and maintained in Dulbecco's modified Eagle's medium (DMEM, Gibco, #11965118) supplemented with 10% foetal bovine serum (FBS, EVERY GREEN, #11011-8611) and 1% penicillin-streptomycin (Gibco, #15140163). Cells were cultured at 37 °C in a humidified atmosphere containing 5% CO 2 . The used cells were confirmed to be free of mycoplasma contamination, and the testing was conducted by Wuhan Servicebio Technology Co., Ltd. Fresh CSCC tissue specimens were collected under sterile conditions and transported in cold, serum-supplemented DMEM. Upon arrival, tissues were washed, minced into small fragments (approximately 0.5–1.0 mm 2 ), and enzymatically digested using collagenase I (Solarbio, #D8071, 2 mg/mL) at 37 °C for 20–40 min under gentle agitation. The resulting cell suspension was filtered through a 70 μm strainer, centrifuged at 400 × g for 5 min, and resuspended in growth factor-reduced basement membrane extract (Bio-Techne China Co. Ltd., #3533-010-02). Droplets containing approximately 5000–10,000 cells were seeded in pre-warmed plates and allowed to polymerise. Organoids were maintained in a defined, serum-free medium based on Advanced DMEM/F12, supplemented with HEPES (Boster, #PYG0019, 1X), GlutaMAX (Thermo Fisher Scientific, #35050061, 1X), antibiotic-antimycotic solution (Solarbio, #P7630, 1X), mycoplasma elimination reagent (Yeasen, #40607ES03, 1X), recombinant human Noggin (PeproTech, #120-10C, 100 ng/mL), fibroblast growth factor 7 (PeproTech, #100-19, 200 ng/mL), SB202190 (MCE, #HY-10295, 1 μM), nicotinamide (MCE, #HY-B0150, 2.5 mM), N-acetylcysteine (MCE, #HY-B0215, 1.25 mM), forskolin (MCE, #HY-15371, 10 μM), B-27 supplement (Gibico, #17504044, 1X), Y-27632 (MCE, #HY-10583, 10 μM), and the TGF-β receptor inhibitor A83-01 (Selleck, #S7692, 500 nM). Organoids were passaged every 10–16 days by mechanical and enzymatic dissociation (TrypLE) and were used experimentally between passages 5–10 to ensure phenotypic consistency. Patient information utilised for organoid construction was provided in Table S6 . The detailed methodologies for calculating IC 50 in CC-derived organoids have been thoroughly documented in our previous study. 35 Briefly, CC-derived organoids were dissociated into single-cell suspensions, and 2000 cells per well were seeded into a 96-well plate. Following 48 h of culture, varying concentrations of mitotane (MCE, HY-13690) or perhexiline (MCE, HY-B1334A) were administrated for 120 h. Adenosine triphosphate (ATP) levels were subsequently quantified using the CellTiter-Glo 3D Viability Assay (Promega, #G9683) in accordance with the manufacturers’ recommendations. IC 50 values were determined using GraphPad Prism 6. Validation of bioinformatic results was experimentally confirmed using multiplex mIF and IHC. For these analyses, 4 μm sections, either paraffin-embedded or frozen, were prepared from collected CC samples. Both methodologies have been extensively delineated in prior literature. 35 Patient-specific details for mIF and IHC were provided in Table S5 . The following primary antibodies were employed: EPCAM (abcam, #ab223582, 1:500), SPRR3 (Atlas, #HPA044467, 1:5000), TNC (Atlas, #HPA004823, 1:100), Snail (abcam, #ab85936, 1 μg/mL), and Ki-67 (abcam, #ab16667, 1:200). Additional reagent specifications were displayed in Table S7 . The siRNA targeting specific genes was designed by Beijing Qingke Biotechnology Co., Ltd. CC cells were seeded in 6-well plates 24 h before transfection to reach approximately 70% confluency. For transfection, 5 μL of Lipofectamine 3000 (Thermo Fisher Scientific, #L3000015) and 10 μL of siRNAs were separately diluted in 250 μL of Opti-MEM (Thermo Fisher Scientific, #31985070), incubated for 5 min at room temperature, and then gently mixed. After a 15 min incubation, 500 μL of the transfection complex was added to each well along with 500 μL of complete DMEM. The medium was refreshed 16 h post-transfection. Cells were harvested after 72 h for mRNA and protein detection of the target genes. The siRNA sequences of target genes were as follows: (1) HADH #1-sense: GACUGGAUACUACGAAGUU; HADH #1-antisense: AACUUCGUAGUAUCCAGUC. (2) HADH #2-sense: CAGAUUACAAGCAUAGCUA; HADH #2-antisense: UAGCUAUGCUUGUAAUCUG. (3) HADHA #1-sense: GGACUAGCUGAUAAGAAGA; HADHA #1-antisense: UCUUCUUAUCAGCUAGUCC. (4) HADHA #2-sense: GGUUAUAUAUGCCGCAAUU; HADHA #2-antisense: AAUUGCGGCAUAUAUAACC. (5) ECHS1 #1-sense: CCGGAGAUCUUAAUAGGAA; ECHS1 #1-antisense: UUCCUAUUAAGAUCUCCGG. (6) ECHS1 #2-sense: GGUGCUAACUUUGAGUACA; ECHS1 #2-antisense: UGUACUCAAAGUUAGCACC. HADH #1-sense: GACUGGAUACUACGAAGUU; HADH #1-antisense: AACUUCGUAGUAUCCAGUC. HADH #2-sense: CAGAUUACAAGCAUAGCUA; HADH #2-antisense: UAGCUAUGCUUGUAAUCUG. HADHA #1-sense: GGACUAGCUGAUAAGAAGA; HADHA #1-antisense: UCUUCUUAUCAGCUAGUCC. HADHA #2-sense: GGUUAUAUAUGCCGCAAUU; HADHA #2-antisense: AAUUGCGGCAUAUAUAACC. ECHS1 #1-sense: CCGGAGAUCUUAAUAGGAA; ECHS1 #1-antisense: UUCCUAUUAAGAUCUCCGG. ECHS1 #2-sense: GGUGCUAACUUUGAGUACA; ECHS1 #2-antisense: UGUACUCAAAGUUAGCACC. CC cell pellets were lysed in RIPA buffer (Servicebio, #G2002) supplemented with protease/phosphatase inhibitors (30 min on ice), followed by sonication and centrifugation (12,000 rpm, 20 min, 4 °C). Protein concentrations were quantified and denatured. Proteins were separated by SDS-PAGE gels (60 V for 30 min, then 120 V for 1.5 h) and transferred to PVDF membranes (280 mA, 2 h, 4 °C). The membranes were blocked with 5% bovine serum albumin (BSA) (30 min, room temperature), incubated with primary antibodies (overnight, 4 °C), washed with tris-buffered saline/Tween 20 (TBST), and probed with HRP-conjugated secondary antibodies (1 h, room temperature). After TBST washes, signals were detected using Western Bright ECL substrate (Advansta, #K-12045-C20). Detailed information of the antibodies used for WB was as follows: E-cadherin (Cell Signalling Technology, #3195, 1:1000), N-cadherin (Cell Signalling Technology, #13116, 1:1000), Snail (Cell Signalling Technology, #3879, 1:1000), HADH (Proteintech, #19828-1-AP, 1:5000), HADHA (Proteintech, #10758-1-AP, 1:5000), ECHS1 (Proteintech, #11305-1-AP, 1:5000), β-Actin (Cell Signalling Technology, #4967, 1:1000) (see Table S7 ). The experimental workflow involved RNA extraction, cDNA synthesis, and qRT-PCR analysis. CC cells were lysed in TRIzol reagent (1 mL), and incubated at room temperature for 5 min followed by chloroform phase separation, isopropanol precipitation, and 75% ethanol washes. RNA concentration was quantified prior to cDNA synthesis using the HiScript II Q RT Super Mix Kit (Vazyme, #R223-01). Target gene expression was analysed via the Bio-Rad CFX96 Real-Time System Manager (C1000 Thermal Cycler) with the qPCR SYBR Green Master Mix (Vazyme, #Q111-02/03), with GAPDH serving as the internal control. The primer sequences used were as follows: (1) HADH -F: TCCGTTGTCCACAGCACAGACT; HADH -R: GGAGGAAGTGTTGCTGGCAAAG. (2) HADHA -F: GCCGACATGGTGATTGAAGCTG; HADHA -R: GGAGAGCAGATGTGTTACTGGC. (3) ECHS1 -F: CCAGGACTGTTACTCCAGCAAG; ECHS1 -R: CACATCATGGCAAGCTCACAGC. (4) GAPDH -F: GTCTCCTCTGACTTCAACAGCG; GAPDH -R: ACCACCCTGTTGCTGTAGCCAA. HADH -F: TCCGTTGTCCACAGCACAGACT; HADH -R: GGAGGAAGTGTTGCTGGCAAAG. HADHA -F: GCCGACATGGTGATTGAAGCTG; HADHA -R: GGAGAGCAGATGTGTTACTGGC. ECHS1 -F: CCAGGACTGTTACTCCAGCAAG; ECHS1 -R: CACATCATGGCAAGCTCACAGC. GAPDH -F: GTCTCCTCTGACTTCAACAGCG; GAPDH -R: ACCACCCTGTTGCTGTAGCCAA. CC cells in logarithmic growth phase were trypsinised, resuspended in phosphate buffer solution (PBS), and counted. Subsequently, cells (1 × 10 3 cells/well) were seeded into 12-well plates and cultured for 14 days (37 °C, 5% CO 2 ), with the medium renewal every three days. Prior to fixation with 4% paraformaldehyde (2 mL/well, 30 min, room temperature), colonies were microscopically examined. After washing with PBS, the colonies were stained with 0.1% crystal violet (20 min, room temperature), rinsed under running water, and air-dried. Finally, the plates were scanned using a colony scanning system, and colony numbers and area were quantified using ImageJ software. CC cells were digested and seeded at a density of 2 × 10 3 cells/well in 96-well plates for overnight culture (37 °C). The following day, the medium was replaced with fresh complete DMEM, and cells were incubated for 24 h. A CCK-8 (Vazyme, #A311-01/02) working solution was prepared by mixing serum-free DMEM with CCK-8 reagent at a 9:1 ratio. After removing the medium, 100 μL of the CCK-8 solution was added to each well. The plates were then incubated at 37 °C in the dark, and OD 450 measurements were taken at 1.5, 2, and 3-h intervals over five consecutive days. Proliferation rates were analysed using GraphPad Prism 6. Experimental conditions included siRNA targeting specific genes ( HADH , HADHA , and ECHS1 ) and concentration gradients of fatty acid oxidation inhibitors (mitotane or perhexiline). Drug concentrations were determined based on manufacturer-recommended protocols established for use in other disease models. CC cells in the logarithmic phase were seeded (4 × 10 4 cells/insert) into 8.0 μm pore size Transwell inserts (Corning, #3422), with 500 μL of complete DMEM added to the lower chambers. After a 16-h incubation at 37 °C, the cells were fixed with 4% paraformaldehyde at room temperature for 15 min and subsequently washed three times with PBS. The lower chambers were then treated with 1 mL of crystal violet stain for 30 min, followed by PBS rinsing and gentle swabbing to remove non-migrated cells. The membranes were air-dried at 37 °C for 10 h, carefully excised with a blade, and mounted (migrated-side up) on slides using neutral resin. Finally, the samples were scanned for ImageJ based quantification of cell penetration. Fatty acid degradation activity was assessed by measuring the NADP + /NADPH ratio, free fatty acid levels, and ATP content using commercial kits (Enhanced NADP + /NADPH Assay Kit, Beyotime #S0180S; Amplex Red Free Fatty Acid Assay Kit, Beyotime #S0215S; ATP Assay Kit, Beyotime #S0026). CC cells (SiHa and HeLa), transfected with targeting siRNAs or treated with pharmacological inhibitors (mitotane or perhexiline), were processed according to kit protocols. Following standard curve generation, working solution preparation, and spectrophotometric analysis, elevated NADP + /NADPH ratios and free fatty acid levels were observed, concurrent with decreased ATP content, indicating impaired fatty acid degradation activity. For xenograft models: female BALB/c nude mice (obtained from Beijing Huafukang Biotechnology Co., Ltd, 4–5 weeks old) were housed under specific pathogen-free (SPF) conditions at the Laboratory Animal Centre, Huazhong University of Science and Technology (Wuhan, China). After seven days of adaptive feeding, nude mice were subcutaneously injected into the right dorsal axilla with 100 μL of PBS suspension containing 5 × 10 6 SiHa cells or 6 × 10 6 HeLa cells. Two weeks post-injection, the mice were randomly allocated into the following groups by a single researcher (uninvolved in subsequent stages): (1) vehicle control for mitotane (n = 5); (2) mitotane treatment (n = 5, 100 mg/kg, intraperitoneally injected once daily); (3) vehicle control for perhexiline (n = 5); (4) perhexiline treatment (n = 5, 80 mg/kg, administered by oral gavage every other day). Tumour volume was measured every two days using a calliper and calculated with the formula: (length × width × width)/2. After three weeks of treatment, the mice were euthanised. Tumour tissues were harvested, weighed, fixed in formalin, and paraffin-embedded for subsequent analysis. For PDX models: female NCG (NOD/ShiLtJGpt-Prkdcem26Cd52Il2rgem26Cd22/Gpt) mice (obtained from Nanjing GemPharmatech Co., Ltd, 4 weeks old) were also housed under SPF conditions. Fresh tumour tissues collected from respective patient donors were mechanically minced into fragments (approximately 3–4 mm 3 ) and subcutaneously implanted into the lower dorsal flank or axilla of the NCG mice. Once stable tumour growth was confirmed, serial transplantation was performed by implanting tumour fragments into new recipient NCG mice. Two weeks post-implantation, mice were randomly assigned to the following groups: (1) vehicle control for mitotane (n = 5); (2) mitotane treatment (n = 5, 100 mg/kg, intraperitoneally injected once daily); (3) vehicle control for perhexiline (n = 5); (4) perhexiline treatment (n = 5, 80 mg/kg, administered by oral gavage every other day). Tumour monitoring and subsequent procedures followed the same protocol as described for the xenograft models above. Statistical analysis was conducted using R (Version 4.1.3). Detailed descriptions of the statistical methods were included in the respective figure legends. A P -value of less than 0.05 was used to determine statistical significance. Asterisks denote significance levels as follows: ns, not significant; ∗ P < 0.05; ∗∗ P < 0.01; ∗∗∗ P < 0.001; ∗∗∗∗ P < 0.0001. This collection of cervical tissues was reviewed and approved by the Medical Ethics Committee of Tongji Medical College, Huazhong University of Science and Technology (Approved number: TJ-IRB20221247) in accordance with the Declaration of Helsinki and relevant policies in China. Written informed consent was obtained from all participants. Animal experiments were conducted by institutional guidelines and were approved by the Institutional Animal Care and Use Committee of Huazhong University of Science and Technology (IACUC number: 5065). All the funders played no roles in study design, data collection, data analysis or writing of the report.

To better understand the spatial architecture and tumour development pattern of CSCC, we conducted ST sequencing on eight samples, including six treatment-naïve patients with CSCC (SM21, SM25, SM27, SM32, SM34, and SM35) and two patients with uterine adenomyosis (N1 and N5), as outlined in Fig. 1 A. The two adenomyosis specimens, obtained from hysterectomy procedures, contained histologically normal cervical epithelium upon pathological examination and were thus designated as normal controls for subsequent comparative analyses. Tissue sections (10 μm thick) from both groups were processed for ST sequencing (10x Genomics Visium) following confirmation of pathological structures via H&E staining. After excluding mitochondrial protein-coding genes, we detected 30,532 spots across the samples, with a median coverage of 2998 genes and 8201 unique molecular identifiers (UMIs) per spot. A professional pathologist (QY) with over fifteen years of expertise in pathology identified tumour foci and normal epithelial regions in the six CSCC samples and two normal samples (bottom panel, Fig. S1A–H ). The distribution of the epithelial marker ( KRT8 ) basically correlated with these manually annotated regions (top panel, Fig. S1A–H ). Notably, sample SM35 contained both normal and tumour foci, as confirmed by pathological assessment and expression of the proliferative marker ( TOP2A ), supporting the reliability of pathologists in delineating tumour foci ( Fig. S1I ). Fig. 1 Integration of ST and scRNA-seq reveals molecular heterogeneity in CC. (A) Schematic overview of the study workflow. (B) UMAP visualisation of 13 transcriptionally distinct clusters, colour-coded by spatial cluster association. (C) Dot plots highlighted cluster-specific marker gene expression. (D) Cluster distribution across eight individual samples. (E) Deconvolution analysis quantified seven cell types within each molecular cluster. Asterisks indicate clusters with enriched proportions of specific cell types relative to others. P- values were calculated using a one-sided Wilcoxon rank-sum test. (F) Spatial cell-type mapping in control (N1 and N5) and tumour (SM21, SM25, SM27, SM32, SM34, and SM35) samples. ∗ P < 0.05. CC: cervical cancer; UMAP: uniform manifold approximation and projection; ST: spatial transcriptomics; scRNA-seq: single-cell RNA-sequencing. Integration of ST and scRNA-seq reveals molecular heterogeneity in CC. (A) Schematic overview of the study workflow. (B) UMAP visualisation of 13 transcriptionally distinct clusters, colour-coded by spatial cluster association. (C) Dot plots highlighted cluster-specific marker gene expression. (D) Cluster distribution across eight individual samples. (E) Deconvolution analysis quantified seven cell types within each molecular cluster. Asterisks indicate clusters with enriched proportions of specific cell types relative to others. P- values were calculated using a one-sided Wilcoxon rank-sum test. (F) Spatial cell-type mapping in control (N1 and N5) and tumour (SM21, SM25, SM27, SM32, SM34, and SM35) samples. ∗ P < 0.05. CC: cervical cancer; UMAP: uniform manifold approximation and projection; ST: spatial transcriptomics; scRNA-seq: single-cell RNA-sequencing. Through dimension reduction and clustering analysis, all spots were divided into 13 molecular clusters, which showed the distinguishing marker genes ( Fig. 1 B and C and Fig. S2A ). Their corresponding DEGs for each cluster were provided in Table S8 (|Log 2 FC| > 0.25, adj- P < 0.05, Wilcoxon rank-sum test). The proportions of these molecular clusters varied significantly across samples ( Fig. 1 D), underscoring substantial intra- and inter-individual heterogeneity in CSCC. To assist the cluster cell type annotation, we deconvoluted the ST data using seven major cell types from a published scRNA-seq data from 20 human cervical tissues (see Methods, Fig. S2B and C ). 22 These cell populations served as references to estimate percentages for spots in ST data using RCTD. 23 Analysis of spatial distribution and marker gene expression of the seven cell types across all spatial spots and individual samples revealed strong concordance between cell type proportions and their UMAP localisation ( Fig. S2D and Fig. S3 ). Notably, epithelial cells exhibited distinct spatial clustering, segregating them from other cell populations. Subsequently, by combining with the clusters’ marker genes ( Fig. 1 C), we annotated the 13 clusters into the following cell types: cluster 1 (B_plasma: CD79A , IGLC2 , IGHG4 and MZB1 ); cluster 2 (stromal_1: CDCA4 , COL16A1 , GPC6 and SWI5 ), cluster 3 (stromal_2: LAMTOR2 and C1GALT1C1 ); cluster 4 (stromal_3: FBXW7 and BRIP1 ), cluster 5 (epithelial cells_1: KRT17 , VEGFA and TM4SF1 ); cluster 6 (epithelial cells_2: KRT15 , MIR205HG and MT1X ); cluster 7 (epithelial cells_3: KRT6A , TRIM29 and KRT14 ); cluster 8 (myoCAFs: COL12A1 , POSTN , MMP11 and FN1 ); cluster 9 (macrophages: C1QB , CTSD and CD68 ); cluster 10 (epithelial cells_4: KRT7 , ECM1 , CNFN and ELF3 ); cluster 11 (endothelial cells: VWF , PECAM1 and PLVAP ); cluster 12 (smooth muscle cells: MYH11 , MUSTN1 and PLN ); and cluster 13 (unknown: NHSL2 and TMUB2 ). By further comparing the abundance of each cell type in each cluster ( Fig. 1 E and Fig. S4A–G ), clusters 5, 6, 7, and 10 were collectively classified as epithelial cells; clusters 2, 3, 8, 12, and 13 as fibroblasts; cluster 1 as immune/fibroblasts; cluster 4 as immune cells; cluster 9 as myeloid cells; and cluster 11 as endothelial cells. Integration of cells by type showed precise alignment between spatially annotated epithelial lesions in ST data and pathological assessment ( Fig. 1 F). Spatial co-localisation analysis revealed extensive interactions between epithelial cells and all other cell types ( Fig. S4H and I ). Together, these findings underscore the remarkable complexity and heterogeneity of epithelial cell populations in CSCC. To further delineate the spatial metabolic landscape of CSCC tumour foci, we performed SM sequencing (AFADESI-MSI) on seven samples (excluding SM21, which failed quality control) using tissue sections adjacent to those employed for ST sequencing ( Fig. 1 A and Fig. S5A–G ). In SM analysis, we identified 1032 negative and 1699 positive metabolites across seven samples. The metabolites exhibited tissue-specific distributions. For example, the metabolite ( m / z = 112.11227, 1-Methyl-1, 3-cyclohexadiene) was predominantly localised in tumour foci (bottom panel, Fig. S5A–G ). Through UMAP, we annotated 15 and 14 clusters for positive and negative ion metabolites, respectively ( Fig. S6A–B and D–E ). The spatial distributions of these clusters closely mirrored the tissue architecture observed in H&E staining ( Fig. S6C and F , and Fig. S5G ), with each cluster exhibiting unique metabolic signatures ( Fig. S6G and H ). These findings demonstrated the utility of spatial-MSI in resolving the spatial regional biological characteristics of CSCC. Next, we asked if there were differences between tumour and normal samples. OPLS-DA revealed clear separation between negative and positive metabolites from normal and tumour samples ( Fig. 2 A and B). 137 differentially expressed metabolites (DEMs) distinguishing the two sample types were identified in Table S9 (VIP> 1, adj- P < 0.05, Wilcoxon rank-sum test), with representative negative and positive DEMs visualised in Fig. 2 C and D, respectively. Functional enrichment analysis indicated that up-regulated DEMs were predominantly involved in unsaturated fatty acids biosynthesis, arginine and proline metabolism, and beta-alanine metabolism ( Fig. 2 E), whereas down-regulated DEMs were primarily linked to ABC transporters, lysine degradation, and galactose metabolism ( Fig. 2 F). These results underscore profound metabolic reprogramming in CSCC development. Fig. 2 ST profiling aligned with transcriptional heterogeneity in CC. (A–B) OPLS-DA of negative (A) and positive (B) ion mode metabolites distinguished tumour and normal tissues. (C–D) Heatmaps of differentially abundant metabolites in negative (C) and positive (D) ion modes. (E–F) Functional enrichment analysis of up-regulated and downregulated metabolites between tumour and normal samples, respectively. ( G) The estimated scores of arginine and proline metabolism pathway across molecular clusters. P- values were calculated using the Wilcoxon signed-rank test. (H) Spatial distribution of arginine and proline metabolism scores in samples N1, N5, SM25, SM27, SM32, SM34, and SM35, respectively. (I) Spatial distribution of certain metabolites linked to arginine and proline metabolism. ∗∗∗∗ P < 0.001. ST: spatial transcriptomics; CC: cervical cancer; OPLS-DA: orthogonal partial least-squares discriminant analysis. ST profiling aligned with transcriptional heterogeneity in CC. (A–B) OPLS-DA of negative (A) and positive (B) ion mode metabolites distinguished tumour and normal tissues. (C–D) Heatmaps of differentially abundant metabolites in negative (C) and positive (D) ion modes. (E–F) Functional enrichment analysis of up-regulated and downregulated metabolites between tumour and normal samples, respectively. ( G) The estimated scores of arginine and proline metabolism pathway across molecular clusters. P- values were calculated using the Wilcoxon signed-rank test. (H) Spatial distribution of arginine and proline metabolism scores in samples N1, N5, SM25, SM27, SM32, SM34, and SM35, respectively. (I) Spatial distribution of certain metabolites linked to arginine and proline metabolism. ∗∗∗∗ P < 0.001. ST: spatial transcriptomics; CC: cervical cancer; OPLS-DA: orthogonal partial least-squares discriminant analysis. To evaluate the consistency between ST and SM data, we focused on the arginine and proline metabolism pathway, which has been reported to be significantly activated in malignancies. 36 , 37 Using the gene set from the KEGG “Arginine and proline metabolism” pathway, we calculated its signature score across all molecular clusters using scMetabolism. Tumour clusters (5, 6, 7 and 10) displayed markedly higher signature scores ( Fig. 2 G). Spatial analysis further demonstrated increased pathway activity within delineated tumour foci (refer to Fig. S1A–H ) ( Fig. 2 H). Additionally, relevant enzymes including SMS ( Fig. S7A and B ) and SAT1 ( Fig. S7C and D ) demonstrated significantly elevated expression in tumour foci, confirming activation of the “Arginine and proline metabolism pathway” in these areas. Correspondingly, downstream metabolites (creatinine, spermine, spermidine, and N-methylhydantoin) were up-regulated in tumour foci, while the initiating metabolite ( l -arginine) exhibited an inverse pattern ( Fig. 2 I). Together, these findings demonstrate strong concordance between ST and spatial-MSI in resolving tissue architecture, particularly within tumour foci, establishing their combined utility for characterising the biological features of CSCC tumour foci. Leveraging the epidermal cell spots identified above, we next performed integrated spatial bio-functional analysis to delineate the growth patterns of tumour foci. Tumour-enriched spots were isolated to quantify spatial continuity and mean malignant cell percentage (see methods, Fig. S8A ). Samples SM27 and SM35 exhibited higher spatial continuity (0.68 and 0.70, respectively) than SM21 and SM34 (0.54 and 0.48; Fig. 3 A). Critically, high-continuity foci demonstrated elevated malignant cell proportions versus low-continuity foci. Specifically, malignant cell percentages in low-continuity samples displayed intra-lesional gradients, whereas high-continuity samples maintained uniformly high proportions ( Fig. 3 B, top), a pattern corroborated by cell density distributions ( Fig. 3 B, bottom). To validate these findings, we analysed 15 independent CSCC samples with Stereo-seq from our prior study. 10 All samples were analogously deconvolved, and further annotated as tumour and non-tumour foci ( Fig. S8B–P ). Consistently, samples TJH37 and TJH90 showed greater spatial continuity and malignant cell proportions than TJH17 and TJH26 ( Fig. 3 C and D). To highlight the molecular and functional differences between samples with high and low continuity, we classified the most extreme samples from both 10x Visium and Stereo-seq datasets into high-purity (n = 4; SM27, SM35, TJH37, TJH90) and low-purity (n = 4, SM21, SM34, TJH17, TJH26) groups based on two key metrics: spatial continuity and cell purity. These grouped samples were then subjected to subsequent analyses. Fig. 3 Samples exhibiting high spatial continuity and epithelial cell percentage revealed epidermal molecular functions. (A) Scatter plot depicted epithelial cell percentage versus spatial continuity for six tumour samples profiled with 10x Visium. (B) Spatial distribution (top) and epithelial cell density (bottom) for representative samples with high (SM27, SM35) and low (SM21, SM34) spatial continuity/epithelial percentage. Vertical dashed lines in density plots indicated the threshold (0.9). (C) Scatter plot of epithelial cell percentage versus spatial continuity for 15 tumour samples profiled with Stereo-seq. (D) Spatial distribution (top) and epithelial cell density (bottom) for representative Stereo-seq samples with high (TJH37, TJH90) and low (TJH17, TJH26) spatial continuity/epithelial percentage. Vertical dashed lines denote the threshold (0.85). (E–F) DEGs analysis between high and low epithelial percentage samples in 10x Visium (E) and Stereo-seq (F) datasets. (G) GO enrichment analysis of genes commonly upregulated in high epithelial percentage samples. (H) Segmented regression analysis correlated epithelial cell percentage with the score of the commonly upregulated gene set. Breakpoints were marked by black vertical dashed lines. DEGs: differentially expressed genes; Stereo-seq: spatial enhanced resolution omics-sequencing; GO: gene ontology. Samples exhibiting high spatial continuity and epithelial cell percentage revealed epidermal molecular functions. (A) Scatter plot depicted epithelial cell percentage versus spatial continuity for six tumour samples profiled with 10x Visium. (B) Spatial distribution (top) and epithelial cell density (bottom) for representative samples with high (SM27, SM35) and low (SM21, SM34) spatial continuity/epithelial percentage. Vertical dashed lines in density plots indicated the threshold (0.9). (C) Scatter plot of epithelial cell percentage versus spatial continuity for 15 tumour samples profiled with Stereo-seq. (D) Spatial distribution (top) and epithelial cell density (bottom) for representative Stereo-seq samples with high (TJH37, TJH90) and low (TJH17, TJH26) spatial continuity/epithelial percentage. Vertical dashed lines denote the threshold (0.85). (E–F) DEGs analysis between high and low epithelial percentage samples in 10x Visium (E) and Stereo-seq (F) datasets. (G) GO enrichment analysis of genes commonly upregulated in high epithelial percentage samples. (H) Segmented regression analysis correlated epithelial cell percentage with the score of the commonly upregulated gene set. Breakpoints were marked by black vertical dashed lines. DEGs: differentially expressed genes; Stereo-seq: spatial enhanced resolution omics-sequencing; GO: gene ontology. Further differential expression analysis between high- and low-purity foci revealed 449 and 470 upregulated genes in 10x Visium and Stereo-seq data, respectively ( Table S10 , Fig. 3 E and F). A 145-gene overlap exhibited high concordance. GO enrichment of these shared genes highlighted “epidermis development” and “epidermal cell differentiation” pathways ( Fig. 3 G). To dissect whether reduced epidermal functions reflected epithelial abundance or intrinsic cellular changes, segmented regression was performed. Above a malignant purity threshold of 0.9, epidermal gene-set scores declined sharply with decreasing purity, transitioning to gradual decay below this breakpoint ( Fig. 3 H). This trend persisted across other CSCC samples ( Fig. S9 ), which indicated that functional erosion arises from altered biological states of malignant cells during progression. Given the enrichment of epidermal functions in high-purity foci, we further probed malignant cell spatial biology in these typical samples. To investigate the spatial structure and developmental dynamics of CSCC tumour foci in detail, we performed high-resolution clustering and functional analysis on samples SM35 and SM27, selected for their higher spatial continuity and purity in tumour foci. Initially, we applied the Cottrazm algorithm to cluster-adjusted gene expression data, identifying 16 distinct clusters in both samples SM35 and SM27 ( Fig. S10A ). 26 Through correlating these clusters with H&E-stained histological features, we delineated tumour-specific clusters (SM35: clusters 5, 6, 9, and 13; SM27: clusters 0, 2, and 9; Fig. 4 A), which exhibited spatially stratified distributions from the core to the periphery of the tumour foci. To assess their developmental progression, we firstly evaluated these clusters using a CIN signature, a biomarker for the translational stage between persistent HPV infection and CSCC carcinogenesis. 38 Notably, inner clusters (cluster 13 in SM35 and cluster 2 in SM27) exhibited significantly higher CIN signature scores compared to outer clusters (cluster 9 in SM35 and cluster 0 in SM27) ( Fig. 4 B). This spatial gradient suggests that the inner foci may represent earlier stages of tumourigenesis, while the outer foci reflect more advanced progression. Fig. 4 Layer-specific transcriptional programs delineated spatial tumour architecture. (A) Spatial segregation of tumour subclusters in samples SM35 (top) and SM27 (bottom). (B) CIN-associated scores across tumour layers in samples SM35 (top) and SM27 (bottom). P -values were calculated using the Kruskal–Wallis test. (C) Heatmap of layer-specific marker genes in samples SM27 and SM35. (D) mIF validated layer-specific SPRR3 and TNC expression in samples SM35 (top) and SM27 (bottom), respectively. (E) GO terms enriched in distinct tumour layers of samples SM27 and SM35. (F) Hypoxia, DNA repair, cell cycle, and EMT scores across varied layers in samples SM27 and SM35, respectively. P -values were calculated using the Kruskal–Wallis test. (G) Spatial annotation (left), tumour layer scores (top right) and functional gene sets (bottom right) in selected tumour foci of TJH37. (H) The correlation between the activity of hypoxia, DNA repair, cell cycle, and EMT and estimated scores of layer 1 (top panel) and layer 3 (bottom panel) in TJH37, respectively. P -values were calculated using the Pearson correlation coefficient. ∗∗∗∗ P < 0.0001. CIN: cervical intraepithelial neoplasia; GO: gene ontology; mIF: multiplex immunofluorescence; EMT: epithelial-to-mesenchymal transition. Layer-specific transcriptional programs delineated spatial tumour architecture. (A) Spatial segregation of tumour subclusters in samples SM35 (top) and SM27 (bottom). (B) CIN-associated scores across tumour layers in samples SM35 (top) and SM27 (bottom). P -values were calculated using the Kruskal–Wallis test. (C) Heatmap of layer-specific marker genes in samples SM27 and SM35. (D) mIF validated layer-specific SPRR3 and TNC expression in samples SM35 (top) and SM27 (bottom), respectively. (E) GO terms enriched in distinct tumour layers of samples SM27 and SM35. (F) Hypoxia, DNA repair, cell cycle, and EMT scores across varied layers in samples SM27 and SM35, respectively. P -values were calculated using the Kruskal–Wallis test. (G) Spatial annotation (left), tumour layer scores (top right) and functional gene sets (bottom right) in selected tumour foci of TJH37. (H) The correlation between the activity of hypoxia, DNA repair, cell cycle, and EMT and estimated scores of layer 1 (top panel) and layer 3 (bottom panel) in TJH37, respectively. P -values were calculated using the Pearson correlation coefficient. ∗∗∗∗ P < 0.0001. CIN: cervical intraepithelial neoplasia; GO: gene ontology; mIF: multiplex immunofluorescence; EMT: epithelial-to-mesenchymal transition. To better characterise the hierarchical structure within both samples, we defined the clusters located at the peripheries of tumour foci as layer 3 (SM35: cluster 9; SM27: cluster 0), those in the central foci as layer 1 (SM35: cluster 13; SM27: cluster 2), and all remaining clusters as layer 2 (SM35: clusters 5 and 6; SM27: cluster 9; Fig. S10B ). We then investigated whether these layers exhibited distinct proximities to the tumour boundary. By analysing the labels of six adjacent spots, we defined boundary spots as those adjacent to at least one non-tumour spot (see Methods, Fig. S10C and F ). The distance of tumour spots was computed based on their spatial proximity to those boundary spots (see Methods , Fig. S10D and G ). Our analysis revealed a progressive decrease in distance to the boundary from layers 1 to 3 ( Fig. S10E and H ), suggesting a clear spatial stratification from the centre to its periphery. To further elucidate the molecular distinctions across different layers, we identified DEGs in layers from samples SM35 and SM27 (|Log 2 FC| > 0.25, adj- P < 0.05, Wilcoxon rank-sum test) ( Table S3 ). Notably, DEGs comparison between samples SM35 and SM27 revealed significant positive correlations in layer 1 (R = 0.55, P < 2.2e-16, Pearson correlation coefficient) and layer 3 (R = 0.39, P = 8.4e-13, Pearson correlation coefficient) ( Fig. S10I and J ), suggesting similar gene expression patterns between these two samples. Specific DEGs common to layers 1, 2, and 3 in both samples were displayed ( Fig. 4 C). The results showed that keratinocyte-related genes, including SPRR3 , SPRR2A , and SFN were predominantly expressed in layer 1, whereas fibroblast-associated and oncogenic genes such as TNC , SERPINB4 , PIK3R1 , and COL7A1 were enriched in layer 3 ( Fig. 4 C and Fig. S11A and B ). mIF staining validated the spatial distribution differences of selected DEGs in SM35 and SM27 samples, further confirming the robustness of our findings ( Fig. 4 D). Functional enrichment analysis demonstrated that layer 1 was primarily involved in pathways related to keratinocyte differentiation, epidermis development, keratinisation, and proteolysis. In contrast, layer 3 was exhibited enrichment in DNA repair, integrin-mediated signalling, response to oxidative stress, and aerobic respiration ( Fig. 4 E). Additionally, assessment of known cancer hallmarks revealed significantly elevated scores for DNA repair, cell cycle, and EMT in layer 3, while layer 1 displayed higher hypoxia scores, suggesting enhanced proliferative and migratory capacities in layer 3 ( Fig. 4 F). To validate spatial stratification, we analysed independent Stereo-seq samples, selecting TJH37 and TJH90 for their high purity and spatial continuity, and annotated them according to their expression profiles ( Fig. 4 G and Fig. S12A ). Layer 1 scores were uniquely elevated in the core tumour foci, whereas layer 3 scores were highest in the outer tumour foci. Functional correlation analysis demonstrated robust positive associations between layer 3 scores and DNA repair, cell cycle and EMT across both samples ( Fig. 4 H and Fig. S12B ). Conversely, hypoxia inversely correlated with layer 3 scores. Additionally, the proportion of proliferating cancer cells (Ki-67 + ) was higher in the outer tumour foci compared to the core foci ( Fig. S12C and D ). Collectively, these findings establish the outer foci (layer 3) as a hub of malignant progression, manifested by decreased epidermal differentiation and enhanced tumour-typical functions such as proliferation, DNA repair and migratory potential. We subsequently investigated the evolutionary trajectories of tumour cells across distinct layers using Monocle2. 29 In samples SM35 and SM27, layers 1, 2, and 3 aligned with the initial, intermediate, and terminal stages of the inferred pseudo-temporal trees, respectively ( Fig. 5 A and C, and Fig. S13A ). Notably, normal epithelial cells were detected at a stage preceding layer 1 in SM35. The expression of layer 1 marker genes (e.g., S100A7, SPRR2A , and SPRR3 ) progressively decreased along pseudo-time, whereas layer 3 marker genes ( MMP11 and TNC ) exhibited an upward trend ( Fig. S13B and C ). Cell-cycle analysis revealed a higher prevalence of cells in G2/M and S phases within layer 3 spots, indicating heightened mitotic activity ( Fig. 5 B and D). Furthermore, we inferred the CNV score by inferCNV, and layer 3 demonstrated significantly higher CNV scores compared to layers 1 and 2 ( Fig. 5 E and F), consistent with the greater frequency CNV events in late-stage tumour cells relative to early-stage counterparts. 39 , 40 Fig. 5 Tumour layer dynamics correlated with genomic evolution and clinical outcomes. (A–D) Pseudo-temporal ordering (A, C) and cell-cycle phase distribution (B, D) across tumour layers and normal epithelia in samples SM35 and SM27. (E–F) CNV scores across different tumour layers in sample SM35 (E) and SM27 (F), respectively. P -values were calculated using the Kruskal–Wallis test. (G) IHC staining of TNC (top) and SPRR3 (bottom) in collected CC samples from the TJH cohort, respectively. (H) Associations of TNC expression with overall survival (left) and progression-free survival (right) in the TJH cohort. P -values were calculated using a univariate Cox regression model. (I) Prognostic impact of the layer 3 gene signature on overall survival (left) and disease-specific survival (right) in TCGA cohort. P -values were calculated using a univariate Cox regression model. ∗∗∗∗ P < 0.0001. CNV: copy number variation; IHC: immunohistochemistry; CC: cervical cancer; TJH: Tongji Hospital; TCGA: The Cancer Genome Atlas. Tumour layer dynamics correlated with genomic evolution and clinical outcomes. (A–D) Pseudo-temporal ordering (A, C) and cell-cycle phase distribution (B, D) across tumour layers and normal epithelia in samples SM35 and SM27. (E–F) CNV scores across different tumour layers in sample SM35 (E) and SM27 (F), respectively. P -values were calculated using the Kruskal–Wallis test. (G) IHC staining of TNC (top) and SPRR3 (bottom) in collected CC samples from the TJH cohort, respectively. (H) Associations of TNC expression with overall survival (left) and progression-free survival (right) in the TJH cohort. P -values were calculated using a univariate Cox regression model. (I) Prognostic impact of the layer 3 gene signature on overall survival (left) and disease-specific survival (right) in TCGA cohort. P -values were calculated using a univariate Cox regression model. ∗∗∗∗ P < 0.0001. CNV: copy number variation; IHC: immunohistochemistry; CC: cervical cancer; TJH: Tongji Hospital; TCGA: The Cancer Genome Atlas. Next, we investigated the transcriptional regulatory mechanisms across different layers using single-cell regulatory network inference and clustering (SCENIC). 30 In both samples SM35 and SM27, five TFs ( ELF3 , PRDM1 , GRHL1 , TEAD2 , and SOX2 ) displayed elevated transcriptional activity and expression levels within the same layers ( Fig. S14A and B ). ELF3 , pivotal in regulating differentiation and homoeostasis, and known to inhibit EMT, was up-regulated in layer 1 ( Fig. S14C ). 41 Additionally, the tumour suppressor TFs, GRHL1 and PRDM1 , also exhibited increased expression in layer 1 ( Fig. S14E and F ). 42 , 43 TEAD2 , crucial for modulating the Hippo-pathway, was up-regulated in layer 3 ( Fig. S14G ). 44 SOX2 , linked to the maintenance of cancer cell stemness, was also found to be up-regulated in layer 3 ( Fig. S14D ). 45 We constructed a gene regulatory network incorporating these five TFs and their target genes ( Fig. S14H ). Subsequent functional enrichment analysis of the target genes indicated that ELF3 and GRHL1 were involved in pathways related to epidermis development, skin development, and epithelial cell differentiation. In contrast, PRDM1 was linked to innate immune responses, defence response to viruses, and defence response to symbionts. TEAD2 primarily participated in pathways related to alcohol metabolic process, actin cytoskeleton organisation, and response to xenobiotic stimulus, whereas SOX2 predominantly influenced the development of neuro projections and organisation of cell projections ( Fig. S14I ). These findings suggest that cells in layer 1 express genes associated with keratinocyte differentiation and tumour suppression, potentially counteracting tumour progression, whereas cells in layer 3 upregulate oncogenic drivers that promote tumourigenesis. Based on the layer-specific findings, we further assessed their prognostic relevance in CC samples from the TJH and TCGA cohorts. SPRR3 and TNC were identified as marker genes for layers 1 and 3, respectively. IHC staining revealed varying expression levels of these proteins in CC samples, categorised as negative/weak (+), positive (++) or strong positive (+++) ( Fig. 5 G). Prognostic analyses demonstrated that patients with high TNC expression exhibited significantly worse overall survival (HR = 6.50, 95% CI = 1.49–28.44, P = 0.013, the univariate Cox regression model) compared to those with low expression, as depicted in Fig. 5 H. In contrast, SPRR3 expression levels did not influence patient prognosis in the same cohort ( Fig. S14J and K ). Furthermore, in the TCGA cohort, patients with an elevated layer 3 signature displayed markedly poorer overall survival (HR = 1.85, 95% CI = 1.15–2.98, P = 0.011, the univariate Cox regression model) and a trend toward worse disease-specific survival (HR = 1.70, 95% CI = 0.99–2.93, P = 0.055, the univariate Cox regression model) ( Fig. 5 I). In general, these findings suggest that layer 3 represents a more aggressive and advanced tumour state, correlating with poorer patient prognosis in CC, a pattern consistent with the features of peripheral tumour cells in the Bayesian state-dependent evolutionary phylodynamic model. 40 Building upon the observed transcriptional variations, we further investigated the metabolic dynamics associated with the tumour foci development. We quantified the activities of four functional modules (proliferation, metabolism, immune, and signalling) based on established certain cancer hallmarks across three distinct layers ( Fig. 6 A). Intriguingly, layer 1 exhibited marked activation of both P53 and KRAS signalling pathways, molecular hallmarks of early tumourigenesis. 39 The evolutionary trajectories of six metabolism-related pathways were mapped for samples SM35 ( Fig. S15A ) and SM27 ( Fig. S15B ), revealing a progressive decline in glycolysis and HAEM metabolism from layers 1 and 2 to layer 3. In contrast, bile acid metabolism, fatty acid metabolism, oxidative phosphorylation, and xenobiotic metabolism exhibited a gradual increase in activity ( Fig. 6 A). Leveraging ST annotation results ( Fig. S10B ), we selected six spatially demarcated tumour foci to identify dynamic metabolite changes ( Fig. 6 B). As demonstrated in Fig. 6 C, 209 negative (R = 0.72, P < 2e-16, Pearson correlation coefficient) and 195 positive metabolites (R = 0.51, P < 2e-16, Pearson correlation coefficient) that displayed intensity shifts from the foci centre to the periphery in both samples SM35 and SM27 samples (see Methods , Table S11 ). Fatty acyls constituted the largest proportion (43%), followed by glycerophospholipids ( Fig. 6 D). Fig. 6 Fatty acid degradation was enriched in invasive tumour layers (layer 3). (A) GSVA scores for certain cancer hallmarks across different tumour layers in samples SM35 and SM27. (B) Layered tumour foci in SM35 (top) and SM27 (bottom). (C) Negative (left) and positive (right) metabolites spatially correlated with tumour layers. P -values were calculated using the Pearson correlation coefficient. (D) Composition of metabolite classes. (E) Fatty acid biosynthesis, elongation, and degradation scores across different tumour layers in samples SM35 (top) and SM27 (bottom). P -values were calculated using the Student's T-test. (F) Schematic of fatty acid degradation. (G) Spatial distributions of FA (16:1), FA (17:1), FA (18:1), FA (20:1), FA (20:2), and FA (24:1) in selected tumour foci with clear spatial demarcation within samples SM35 (top) and SM27 (bottom), respectively. Different coloured rings represented distinct tumour foci. P -values were calculated using the Student's T-test. (H) Layer-specific expression of fatty acid degradation enzymes ( CPT1A , CPT2 , HADH , HADHA , and ECHS1 ) across different tumour layers in samples SM35 (top) and SM27 (bottom), respectively. P -values were calculated using the Kruskal–Wallis test. ∗ P < 0.05, ∗∗ P < 0.01, ∗∗∗ P < 0.001, ∗∗∗∗ P < 0.0001, ns, not significant. GSVA: gene set variation analysis. Fatty acid degradation was enriched in invasive tumour layers (layer 3). (A) GSVA scores for certain cancer hallmarks across different tumour layers in samples SM35 and SM27. (B) Layered tumour foci in SM35 (top) and SM27 (bottom). (C) Negative (left) and positive (right) metabolites spatially correlated with tumour layers. P -values were calculated using the Pearson correlation coefficient. (D) Composition of metabolite classes. (E) Fatty acid biosynthesis, elongation, and degradation scores across different tumour layers in samples SM35 (top) and SM27 (bottom). P -values were calculated using the Student's T-test. (F) Schematic of fatty acid degradation. (G) Spatial distributions of FA (16:1), FA (17:1), FA (18:1), FA (20:1), FA (20:2), and FA (24:1) in selected tumour foci with clear spatial demarcation within samples SM35 (top) and SM27 (bottom), respectively. Different coloured rings represented distinct tumour foci. P -values were calculated using the Student's T-test. (H) Layer-specific expression of fatty acid degradation enzymes ( CPT1A , CPT2 , HADH , HADHA , and ECHS1 ) across different tumour layers in samples SM35 (top) and SM27 (bottom), respectively. P -values were calculated using the Kruskal–Wallis test. ∗ P < 0.05, ∗∗ P < 0.01, ∗∗∗ P < 0.001, ∗∗∗∗ P < 0.0001, ns, not significant. GSVA: gene set variation analysis. Given these findings of ST and SM data, we focused on fatty acid metabolism-related processes in the tumour development. As illustrated in Fig. 6 E and Fig. S16A and B , we evaluated three key pathways: fatty acid biosynthesis, elongation, and degradation. Notably, only fatty acid degradation showed progressive activation from layers 1 to 3 in both samples SM35 and SM27, whereas the other two processes remained relatively stable. Samples with Stereo-seq confirmed elevated fatty acid degradation activity in high layer 3 foci ( Fig. S16C and D ). Correlative mapping demonstrated a positive association between fatty acid degradation activity and the layer 3 signature score in TJH37 (R = 0.1, P = 2.59e-02; Pearson correlation coefficient) and TJH90 (R = 0.19, P = 2.27e-07; Pearson correlation coefficient). Conversely, foci enriched with layer 1 signatures showed an inverse metabolic pattern (TJH 37: R = −0.42, P = 8.94e-22; TJH90: R = −0.37, P = 1.30e-24; Fig. S16E and F ). These findings suggest that fatty acid degradation is indeed associated with tumour spatial development. We subsequently focused on enzyme-encoding genes (e.g., CPT1A , HADH , and ECHS1 ) and specific metabolites associated with this pathway ( Fig. 6 F). In brief, fatty acids of varying chain lengths are initially activated through binding to coenzyme A (CoA) in the cellular cytoplasm. Carnitine palmitoyltransferase 1 (CPT1) then converts long-chain fatty acyl-CoA into acylcarnitine, while carnitine transferase (CAT) facilitates its transport across the mitochondrial inner membrane. Subsequently, long-chain acylcarnitine is reconverted into long-chain acylcarnitine CoA by CPT2 before undergoing β-oxidation. Ultimately, long-chain acylcarnitine CoA is further degraded into acetyl-CoA, FADH2 and NADH. 46 Notably, the concentrations of multiple fatty acids (FAs), specifically FA (16:1), FA (17:1), FA (18:1), FA (20:1), FA (20:2), and FA (24:1), were significantly higher in the inner tumour foci (layer 1) than in the peripheral region (layer 3) ( Fig. 6 G). To further quantify this spatial pattern, we grouped tumour spots from the spatial metabolomics data into three categories (Near, Medium, and Far) based on their Euclidean distance to the tumour boundary, which correspond respectively to layer 3, layer 2, and layer 1 as defined in the transcriptomic analysis. Consistent with the layer-based comparison, metabolite intensities increased progressively from the Near to the Far group, confirming that fatty acid levels are indeed elevated in the innermost tumour compartment. Furthermore, the expression levels of key enzyme-encoding genes ( CPT1A , CPT2 , HADH , HADHA , and ECHS1) were significantly higher in layer 3 than in other layers in samples SM35 (top) and SM27 (bottom) ( Fig. 6 H). Considering the consistent molecular activity and related metabolites of fatty acid degradation pathway, the increased degradation of fatty acids in layer 3 suggests a critical metabolic adaptation in the evolutionary trajectory of tumours to support progression. To functionally investigate whether the fatty acid degradation pathway, which was specifically activated in the outer tumour layer (layer 3), contributes to the malignant behaviour characteristic of that layer, we performed experimental verification in vitro and in vivo. We first designed siRNA targeting key enzymes genes ( HADH , HADHA , and ECHS1 ) involved in fatty acid degradation in two CC cell lines ( Fig. S17A ). Genetic perturbation significantly reduced fatty acid degradation activity, evidenced by elevated NADP + /NADPH ratios, elevated content of free fatty acid, decreased content of ATP ( Fig. S17B–D ). As shown in Fig. 7 A and Fig. S17E , the proliferative capacities of SiHa and HeLa cells were markedly reduced following the knockdown of HADH , HADHA , and ECHS1 , respectively. Additionally, colony formation was also attenuated ( Fig. 7 C and Fig. S17G ). Migratory assays revealed that the migratory abilities of SiHa and HeLa cells were significantly impaired upon downregulation of HADH , HADHA , and ECHS1 ( Fig. 7 B and Fig. S17F ). Further experimental methods confirmed key characteristics of layer 3, including proliferation and EMT, was altered in conjunction with the suppression of fatty acid degradation activity ( Fig. 7 D and Fig. S17H ). Fig. 7 Genetic or pharmacological inhibition of fatty acid oxidation suppresses CC aggressiveness. (A) Proliferation assays were performed following knockdown of HADH (left), HADHA (middle), and ECHS1 (right) in SiHa cells. Each group had three independent biological replicates. P -values were calculated using the One-way ANOVA. (B) Knockdown of HADH (top), HADHA (middle), and ECHS1 (bottom) could remarkably inhibit the capacities of migration in SiHa cells. Transwell assay had three independent biological replicates. P -values were calculated using the One-way ANOVA. (C) Knockdown of HADH (top), HADHA (middle), and ECHS1 (bottom) could remarkably inhibit the capacities of colony formation in SiHa cells. Colony formation had three independent biological replicates. P -values were calculated using the One-way ANOVA. (D) WB analysis of EMT markers (N-cadherin, E-cadherin, Snail) after silencing HADH (left), HADHA (middle), and ECHS1 (right) in SiHa cells, respectively. (E) Dose-dependent suppression of SiHa cells by mitotane (top) and perhexiline (bottom) treatment. Each group had three independent biological replicates. P -values were calculated using the One-way ANOVA. (F–G) Different concentrations of mitotane (F) and perhexiline (G) could remarkably inhibit the capacities of colony formation in SiHa cells, respectively. Colony formation had three independent biological replicates. P -values were calculated using the One-way ANOVA. (H) Different concentrations of mitotane (left) and perhexiline (right) could remarkably inhibit the capacities of migration in SiHa cells. Transwell assay had three independent biological replicates. P -values were calculated using the One-way ANOVA. (I) WB analysis of EMT markers (N-cadherin, E-cadherin, Snail) after mitotane (left) and perhexiline (right) treatment in SiHa cells. ∗ P < 0.05, ∗∗ P < 0.01, ∗∗∗ P < 0.001, ∗∗∗∗ P < 0.0001, ns, not significant. CC: cervical cancer; WB: western blotting; EMT: epithelial-to-mesenchymal transition. Genetic or pharmacological inhibition of fatty acid oxidation suppresses CC aggressiveness. (A) Proliferation assays were performed following knockdown of HADH (left), HADHA (middle), and ECHS1 (right) in SiHa cells. Each group had three independent biological replicates. P -values were calculated using the One-way ANOVA. (B) Knockdown of HADH (top), HADHA (middle), and ECHS1 (bottom) could remarkably inhibit the capacities of migration in SiHa cells. Transwell assay had three independent biological replicates. P -values were calculated using the One-way ANOVA. (C) Knockdown of HADH (top), HADHA (middle), and ECHS1 (bottom) could remarkably inhibit the capacities of colony formation in SiHa cells. Colony formation had three independent biological replicates. P -values were calculated using the One-way ANOVA. (D) WB analysis of EMT markers (N-cadherin, E-cadherin, Snail) after silencing HADH (left), HADHA (middle), and ECHS1 (right) in SiHa cells, respectively. (E) Dose-dependent suppression of SiHa cells by mitotane (top) and perhexiline (bottom) treatment. Each group had three independent biological replicates. P -values were calculated using the One-way ANOVA. (F–G) Different concentrations of mitotane (F) and perhexiline (G) could remarkably inhibit the capacities of colony formation in SiHa cells, respectively. Colony formation had three independent biological replicates. P -values were calculated using the One-way ANOVA. (H) Different concentrations of mitotane (left) and perhexiline (right) could remarkably inhibit the capacities of migration in SiHa cells. Transwell assay had three independent biological replicates. P -values were calculated using the One-way ANOVA. (I) WB analysis of EMT markers (N-cadherin, E-cadherin, Snail) after mitotane (left) and perhexiline (right) treatment in SiHa cells. ∗ P < 0.05, ∗∗ P < 0.01, ∗∗∗ P < 0.001, ∗∗∗∗ P < 0.0001, ns, not significant. CC: cervical cancer; WB: western blotting; EMT: epithelial-to-mesenchymal transition. Then, we also investigated whether the small-molecule inhibitors (mitotane and perhexiline) of fatty acid degradation exerted similar biological effects on CC malignant transformation. In lines with established mechanisms, both Mitotane and perhexiline significantly suppressed fatty acid degradation activity in SiHa and HeLa cells ( Fig. S18A–C ). Functionally, treatment with varying concentrations of mitotane and perhexiline consistently reduced the proliferative ( Fig. 7 E and Fig. S18D ) and colony-forming ( Fig. 7 F and G and Fig. S18F ) abilities of SiHa and HeLa cells. Additionally, cell migration was markedly inhibited ( Fig. 7 H and Fig. S18E ). Furthermore, the process of EMT was inhibited upon mitotane and perhexiline treatment, as evidenced by increased E-cadherin/Snail expression and decreased N-cadherin levels ( Fig. 7 I and Fig. S18G ). To extend these findings to more physiologically relevant models, we evaluated the therapeutic potential of mitotane and perhexiline in CC-derived organoids and animal models. As shown in Fig. 8 A and B, both mitotane and perhexiline displayed robust growth-suppressive effects on six PDOs by reducing the number and diameters of organoids, with IC 50 values ranging from 18.26 μM to 39.03 μM for mitotane, and 3.82 μM–6.65 μM for perhexiline ( Fig. 8 C). mIF staining revealed that the EMT process was attenuated following perhexiline treatment, accompanied by alterations in key EMT markers, including E-cadherin, N-cadherin, and Snail ( Fig. 8 D). In subcutaneous CC xenograft models, we found that both mitotane ( Fig. 8 E and F and Fig. S18H and I ) and perhexiline ( Fig. 8 G and H and Fig. S18J and K ) effectively suppressed tumour cell growth, as evidenced by significantly reduced tumour weight and volume in the treatment groups compared to the vehicle controls. These antitumour effects were further validated in a CC-derived PDX model ( Fig. 8 I and L). Taking together, these results indicate that inhibition of fatty acid degradation mitigates CC aggressiveness, concomitant with decreased proliferation and suppression of EMT traits. Fig. 8 Pharmacological inhibition of fatty acid oxidation suppresses CC progression in organoid and mouse models. (A) Morphological changes of CC-derived organoids following mitotane (top) and perhexiline (bottom) treatment. (B) Changes of number and diameter in CC-derived organoids following mitotane and perhexiline treatment. P -values were calculated using the One-way ANOVA. (C) The IC 50 values of mitotane (top) and perhexiline (bottom) in six CC-derived organoids. Each drug had at least three biological replicates in each constructed organoid. (D) The changes of E-cadherin (left), N-cadherin (middle), Snail1 (right) and Ki-67 expression following perhexiline treatment in CC-derived organoids. (E) Effect of mitotane treatment (n = 5) on tumour growth in the subcutaneous CC model in nude mice. P -values were calculated using the Wilcoxon rank-sum test. (F) Effect of mitotane treatment on tumour volume growth in each group of mice. P -values were calculated using the Student's T-test. (G) Effect of perhexiline treatment (n = 5) on tumour growth in the subcutaneous CC model in nude mice. P -values were calculated using the Student's T-test. (H) Effect of perhexiline treatment on tumour volume growth in each group of mice. P -values were calculated using the Student's T-test. (I) Effect of mitotane treatment (n = 5) on tumour growth in a CC-derived PDX model. P -values were calculated using the Student's T-test. (J) Effect of mitotane treatment on tumour volume growth in each group of the PDX model. P -values were calculated using the Student's T-test. (K) Effect of perhexiline treatment (n = 5) on tumour growth in a CC-derived PDX model. P -values were calculated using the Student's T-test. (L) Effect of perhexiline treatment on tumour volume growth in each group of the PDX model. P -values were calculated using the Wilcoxon rank-sum test. ∗ P < 0.05, ∗∗ P < 0.01, ∗∗∗ P < 0.001, ∗∗∗∗ P < 0.0001, ns, not significant. CC: cervical cancer; IC 50 : half maximal inhibitory concentration; PDX: patient-derived xenograft. Pharmacological inhibition of fatty acid oxidation suppresses CC progression in organoid and mouse models. (A) Morphological changes of CC-derived organoids following mitotane (top) and perhexiline (bottom) treatment. (B) Changes of number and diameter in CC-derived organoids following mitotane and perhexiline treatment. P -values were calculated using the One-way ANOVA. (C) The IC 50 values of mitotane (top) and perhexiline (bottom) in six CC-derived organoids. Each drug had at least three biological replicates in each constructed organoid. (D) The changes of E-cadherin (left), N-cadherin (middle), Snail1 (right) and Ki-67 expression following perhexiline treatment in CC-derived organoids. (E) Effect of mitotane treatment (n = 5) on tumour growth in the subcutaneous CC model in nude mice. P -values were calculated using the Wilcoxon rank-sum test. (F) Effect of mitotane treatment on tumour volume growth in each group of mice. P -values were calculated using the Student's T-test. (G) Effect of perhexiline treatment (n = 5) on tumour growth in the subcutaneous CC model in nude mice. P -values were calculated using the Student's T-test. (H) Effect of perhexiline treatment on tumour volume growth in each group of mice. P -values were calculated using the Student's T-test. (I) Effect of mitotane treatment (n = 5) on tumour growth in a CC-derived PDX model. P -values were calculated using the Student's T-test. (J) Effect of mitotane treatment on tumour volume growth in each group of the PDX model. P -values were calculated using the Student's T-test. (K) Effect of perhexiline treatment (n = 5) on tumour growth in a CC-derived PDX model. P -values were calculated using the Student's T-test. (L) Effect of perhexiline treatment on tumour volume growth in each group of the PDX model. P -values were calculated using the Wilcoxon rank-sum test. ∗ P < 0.05, ∗∗ P < 0.01, ∗∗∗ P < 0.001, ∗∗∗∗ P < 0.0001, ns, not significant. CC: cervical cancer; IC 50 : half maximal inhibitory concentration; PDX: patient-derived xenograft.

Intra-tumoural heterogeneity, a defining feature of malignant tumours, arises from the dynamic accumulation of molecular and genetic alterations during tumourigenesis, ultimately driving divergent phenotypic behaviours such as differential proliferation rates, invasive potential, and therapeutic resistance. 47 Despite its well-documented role in tumour progression, the spatial architecture of such heterogeneity in CSCC remains poorly characterised. In this study, we employed an integrative multi-omics approach, combining ST and spatial-MSI sequencing on a cohort of eight samples (two normal cervical tissues and six CSCC samples), demonstrating a high degree of concordance between these two omics modalities at an unprecedented resolution. Notably, our analyses uncovered profound transcriptomic and metabolomic stratification across distinct tumour foci. To our knowledge, this work represents the application of dual ST and spatial-MSI profiling in CSCC samples, providing a foundational resource for deciphering the spatial biology of cervical carcinogenesis. The evolutionary trajectories by which monoclonal tumour populations diversify into spatially and functionally distinct subclones remain incompletely understood. By focussing on two spatially continuous samples (SM35 and SM27), we delineated layer-specific malignant programmes and observed a clear spatial hierarchy. Functional enrichment analysis revealed tumour core regions (layer 1) exhibited prominent epithelial differentiation signatures, including keratinisation and epidermis development, whereas invasive fronts (layer 3) were characterised by metabolic reprogramming, cell cycle progression, and DNA repair, and represented the terminal state with increased CNV burden along pseudo-time. Although the high expression of differentiation-associated genes and tumour suppressor-related transcription factors in the tumour core may appear counterintuitive, this pattern is consistent with recent spatial transcriptomics studies in other malignancies. In oral squamous cell carcinoma, tumour core malignant cells displayed strong keratinocyte differentiation programmes and were considered to represent an early stage of tumourigenesis. 9 Similarly, in oesophageal squamous cell carcinoma, keratinisation- and cornified envelope-related functions peaked at the low-grade intraepithelial neoplasia stage and progressively attenuated with tumour progression, with these signatures being associated with favourable prognosis. 48 In line with these observations, the tumour suppressor-associated transcription factors identified in layer 1, including GRHL1 and ELF3 , are well-established regulators of keratinocyte differentiation and epithelial homoeostasis. 49 , 50 Their elevated expression in the central tumour region therefore likely reflects the retention of differentiation-associated regulatory programmes characteristic of early or less aggressive malignant states, rather than effective tumour suppression. Meanwhile, layer 2 emerged as a dynamic intermediate compartment along the spatial evolutionary trajectory, occupying a transitional pseudo-time position without forming an independent branch, and displaying intermediate CNV levels, cancer hallmark activities and metabolic features. Notably, fatty acid degradation activity increased progressively from layer 1 through layer 2 to layer 3, supporting a stepwise metabolic reprogramming process. Together, these findings indicate that layer 2 represents a transitional tumour state bridging a differentiation-dominant core and a metabolically reprogrammed, proliferative and invasive periphery, highlighting the continuous and progressive nature of spatial tumour evolution in CSCC. Malignant progression necessitates metabolic adaptations to fulfil biosynthetic and energetic demands. Although the Warburg effect dominates a classical paradigm, our spatial metabolomics analysis identified dysregulated lipid metabolism as a hallmark of invasive tumour fronts. Leveraging two spatial omics datasets, our study demonstrated that peripheral tumour cells exhibited higher fatty acid degradation activity, with up-regulation of β-oxidation enzymes (e.g., CPT1A , HADHA , ECHS1 ) coupled with depletion of long-chain fatty acids, a pattern further validated in independent CSCC samples using Stereo-seq. Mechanistically, mitochondrial fatty acid catabolism not only sustains ATP production but also mitigates oxidative damage during metastatic dissemination. 51 , 52 EMT serves as a pivotal driver of malignant progression by disrupting cell adhesion (e.g., downregulating E-cadherin), remodelling the cytoskeleton, and activating matrix metalloproteinases, thereby endowing cancer cells with invasive properties. 53 As another obviously biological characteristic of peripheral tumour, dysregulated EMT activation was identified in our prior CC-related studies. 54 , 55 Here, we observed that tumour cells in layer 3 exhibited higher EMT scores in both our cohorts and external samples using Stereo-seq. This metabolic adaptation appears closely linked to EMT, as invasive fronts exhibited elevated EMT scores concurrent with fatty acid degradation signatures, a nexus previously observed in colorectal carcinoma and pancreatic ductal adenocarcinoma. 56 , 57 This association compels us to investigate the underlying mechanisms that extend beyond a simple correlation. We hypothesise that fatty acid degradation may propel the EMT not by directly modulating classical adhesion molecules like E-cadherin or N-cadherin, but rather through intermediary metabolic and signalling relays. The bioenergetic and redox-supporting outputs of this pathway, namely ATP and NADPH, 57 could furnish the necessary foundation for cellular migration and invasion. Moreover, the degradation process itself or its byproducts may act as signalling molecules and engage key kinase cascades, such as the p38/MAPK pathway and other MAPK cascades, 58 , 59 which are established upstream regulators of the EMT transcriptional network. In our study, the regulatory interplay between fatty acid degradation and the acquisition of the EMT phenotype requires further elucidation in future investigations. The mechanistic link between fatty acid degradation and EMT prompted therapeutic exploration. To investigate the biological effects of fatty acid degradation on peripheral tumour foci (layer 3) in CC cells and CC-derived organoids, we employed two interference approaches: genetic silencing (siRNAs targeting HADH , HADHA , ECHS1 ) and pharmacological inhibitors (mitotane and perhexiline). Key related enzyme-encoding genes involved in fatty acid degradation ( CPT1A , ECHS1 , and HADHA ) have been shown to be up-regulated in multiple malignancies and associated with poor patient prognosis. 60 , 61 , 62 , 63 Additionally, mitotane and perhexiline were reported to interfere with the biological effects of CPT1 and CPT2 enzymes, exhibiting antitumour effects in malignant tumours. 64 , 65 In our study, both intervention methods disrupted redox balance (NADP + /NADPH), energy homoeostasis (ATP), and lipid pools (free fatty acids) in CSCC models. These interventions significantly suppressed the proliferative capacity and aggressive traits of CSCC cells, organoids and xenografts. Furthermore, the iconic feature of peripheral tumour foci, EMT, was also obviously reversed following these two interfere methods. Notably, our findings expand the potential repurposing of mitotane and perhexiline as therapeutic strategies for CSCC. Mitotane is a well-established adrenolytic agent whose clinical use is exclusively confined to the treatment of adrenocortical carcinoma, 66 functioning primarily through the disruption of mitochondrial function and steroidogenesis. Crucially, its mechanism involves the inhibition of CPT1, the rate-limiting enzyme for fatty acid degradation. Similarly, perhexiline is a clinically approved anti-anginal agent that inhibits myocardial fatty acid oxidation via direct antagonism of both CPT1 and CPT2 to shift cardiac metabolism toward glucose utilisation. 67 While preclinical evidence indicates that perhexiline can inhibit proliferation in models of glioblastoma and chronic lymphocytic leukaemia, 68 , 69 its primary clinical indication remains non-oncological. Our spatial multi-omics discovery that the aggressive peripheral tumour foci (layer 3) in CSCC are critically dependent on activated fatty acid degradation provides the mechanistic rationale for repurposing these agents. Their demonstrated efficacy in our preclinical CSCC models, attenuating proliferation, invasion, and EMT phenotypes-therefore represents a significant conceptual extension beyond their established clinical or mechanistic frameworks and highlights a spatially informed therapeutic vulnerability based on tumour architecture. Importantly, the rationale for repurposing these drugs is strongly supported by a dual foundation. Mechanistically, their known on-target inhibition directly counters the aggressive metabolic pathway we identified. Practically, their safety profiles are clinically manageable, with established protocols for monitoring mitotane (drug/hormone levels) and perhexiline (CYP2D6 genotype/plasma concentration). This combination of targeted mechanism and established clinical safety protocols establishes a viable foundation for future translational research in defined CSCC patient populations, such as individuals with refractory or metastatic disease. In addition to the significant findings highlighted above, this study possesses several limitations that warrant consideration. First, cohort size constraints necessitate validation in expanded clinical datasets to assess the universality of our spatial stratification model. We will initiate a multi-centre study to expand the sample size, encompassing diverse CC pathological subtypes and stages to validate and refine our spatial model. Second, the functional diversity of epithelial subpopulations could be further resolved via fluorescence-activated cell sorting using layer-specific markers. Finally, although our preclinical models establish proof-of-concept efficacy, advancing fatty acid oxidation inhibitors toward clinical application necessitates a two-pronged translational strategy: first, rigorous validation in comprehensive in vivo patient-derived xenograft models to confirm therapeutic efficacy and safety; second, the design and execution of a biomarker-guided Phase Ib/IIa clinical trial in patients with refractory or metastatic CSCC to evaluate the clinical feasibility of targeting this spatially defined metabolic vulnerability.

PW, CW, FR, and LW contributed to the study conception and design; SL, ZL, YL, BH, and MQ collected cervical cancer samples, conducted bioinformatic analyses and performed laboratory experiments. SL and ZL wrote the manuscript. BL, CC, TP, MX, YX, XL, XW, LL, WD, ZX, RL, and JZ supervised this work. SL and PW have directly accessed and verified the underlying data reported in this manuscript. All authors reviewed and approved the final manuscript.

Cervical cancer (CC) persists as a leading cause of cancer-related mortality among women worldwide, disproportionately affecting low- and middle-income countries. In 2022, China recorded approximately 150,700 new cases and 55,700 deaths from CC, 1 despite advancements in human papillomavirus (HPV) vaccination and screening programmes. Economic disparities limit access to preventive care, resulting in frequent late-stage diagnoses where 5-year survival rates plummet to 16%–58%. 2 Tumour heterogeneity, a hallmark of therapeutic resistance and metastatic potential, 3 , 4 drives this dismal prognosis. The precise mechanism underlying the development of heterogeneity in tumour cells during cancer initiation and progression remains somewhat elusive. Spatial biology combines the molecular characteristics of cells with spatial information, allowing us to analyse specific cancer foci and providing an opportunity for the study of cancer cell development in CC. 5 In recent years, spatial transcriptomics (ST) has been shown in many studies to help analyse tumour development. A cross-tissue ST study revealed that tumour cell clonal branching correlates strongly with spatial architecture, providing the first evidence for reconstructing tumour evolution through ST approaches. 6 Noteworthy is the use of spatial enhanced resolution omics-sequencing (Stereo-seq) and 10x Genomics Visium data, which have revealed significant cellular infiltration and invasion at the boundary of hepatocellular carcinoma. 7 , 8 Research also suggests that cancer cells from the leading edge of tumour exhibit a higher estimated epithelial-to-mesenchymal transition (EMT), along with increased invasion and metastasis capacities, compared to cells from the tumour core. 9 In fact, our prior work showed that CC also exhibits relatively large heterogeneity across tissues, with clear spatial clustering of tumour foci, and we found POSTN + myofibroblastic cancer-associated fibroblasts (myoCAFs) around the leading edge of tumours, potentially contributing to cervical squamous cell carcinoma (CSCC) progression. 10 However, the developmental pattern inside the tumour foci, from quiescent cores to invasive edges, remains unexplored in CSCC. Metabolic reprogramming, a cardinal cancer hallmark, intersects with transcriptional networks to sustain tumour development. 11 , 12 Recent advances in ST and spatial mass spectrometry imaging (spatial-MSI) now enable multi-omics mapping of intact tissues, offering unprecedented resolution to dissect tumour development mechanisms. 13 The study of glioblastomas has unveiled critical tumour–host interactions and adaptive transcriptional programmes through the integration of ST and spatial metabolomics (SM). 14 Similarly, Sun et al. constructed a transcriptional and metabolic atlas of gastric cancer using ST-SM integration, revealing localised metabolic rewiring during tumour progression. 15 In oral squamous cell carcinoma, combined ST and SM analyses identified dysregulated polyamine metabolism as a key driver of tumourigenesis. 16 In CC, lipid metabolism enzymes like FABP5 and FASN promote lympho-vascular invasion, 17 , 18 while dysregulated fatty acid oxidation modulates therapeutic resistance. 19 , 20 The spatial metabolic landscape of CC remains poorly understood. Integrating ST and SM could elucidate the spatially resolved mechanisms of CC progression and uncover precision therapeutic targets. In this study, we comprehensively mapped the transcriptional and metabolic profiles of CC by integrating ST and spatial-MSI across eight samples (six CSCC carcinomas and two normal cervical tissues), alongside 15 additional CC samples previously analysed by Stereo-seq for validation. Through systematic investigation of spatial heterogeneity in CC, we demonstrated that samples exhibiting high spatial continuity of tumour foci effectively recapitulate the transitional process from early neoplastic development to malignant cell maturation. Notably, this spatial state transition displayed a pronounced hierarchical organisation. Correlative analysis with SM data further revealed that dynamic alterations in metabolites and associated enzymatic activities implicate fatty acid degradation as a critical regulator of cancer cell state progression. Functional validation with CC cell lines, organoids and patient-derived xenograft (PDX) demonstrated that targeting this metabolic pathway suppressed malignant phenotypes, nominating its potential as a therapeutic vulnerability. Our study establishes spatial architecture as a key determinant of CC progression and provides a framework for mechanistically informed therapeutic intervention.

The authors declare that they have no competing interests.