Abstract

   While cervical cancer (CC) prognosis depends on tumor staging, the
   spatiotemporal evolution of tumor microenvironment (TME) heterogeneity
   during metastatic progression remains poorly characterized at
   single-cell resolution. We employed an integrative multi-omics
   approach, combining single-cell RNA sequencing (scRNA-seq; n = 11),
   spatial transcriptomics (ST), and bulk RNA-seq data from the TCGA-CESC
   cohort (n = 304), to systematically map TME remodeling across CC
   progression stages. scRNA-seq was performed on primary lesions from
   patients with localized (n = 3), regional (n = 4), and metastatic (n =
   4) diseases, with in-depth analyses focusing on cellular
   characteristics, cell type composition alterations, functional changes,
   differentiation trajectories, and cell-cell interaction networks. These
   findings were further validated using spatial transcriptomics, bulk
   RNA-seq data, and multiple immunohistochemistry (mIHC) experiments.
   ScRNA-seq data revealed that the TME of the metastatic group displayed
   a distinct immunosuppressive phenotype. Three key subclusters closely
   linked to TME remodeling in this group were identified. Notably, a
   novel metastasis-associated epithelial subpopulation (Epi0_AGR2),
   characterized by both epithelial-mesenchymal transition (EMT) and
   chemokine secretory phenotypes, was discovered. Gene Set Variation
   Analysis (GSVA) revealed that transforming growth factor β (TGF-β)
   signaling activation served as its primary transcriptional driver.
   Additionally, a neutrophil subset with pro-tumor and immunosuppressive
   properties, as well as a cancer-associated fibroblasts (CAFs) subset
   that promoted angiogenesis, were enriched in the metastatic group.
   Cell-cell interaction analysis and spatial mapping further revealed the
   formation of coordinated Epi0-neutrophil-CAFs niches, which established
   TGF-β-CXCL1/2/8-OSM/OSMR feedforward loops. Importantly, a
   computational model derived from the TME metastatic niche signature
   demonstrated significant prognostic stratification in the TCGA cohort
   (HR = 2.5179, p = 0.0144). In all, this study provides the first
   comprehensive delineation of stage-specific TME dynamics in CC,
   revealing TGF-β-driven cellular cooperativity as a metastatic switch.
   The joint framework establishes a potential clinically translatable
   tool for precision prognosis and therapeutic targeting.

   Keywords: cervical cancer, tumor progression, multi-omics analysis,
   TME, niche

Introduction

   Cancer metastasis is associated with tumor progression and lethality in
   cases of cervical cancer (CC) [49]^1. The five-year survival rates of
   patients with stage III and stage IV CC were only 10-40% [50]^2. The
   recent introduction of immune checkpoint inhibitors (ICIs) as a
   treatment for advanced CC patients resulted in a significant
   improvement in their progression-free survival and overall survival
   [51]^3^-[52]^6. The response to ICIs has been shown to be closely
   associated with the initial tumor microenvironment (TME) [53]^7^,
   [54]^8. Therefore, a comprehensive understanding of the TME in advanced
   CC, and the exploration of the underlying mechanisms that promote
   cancer progression and immunosuppression, will offer new insights on
   the treatment and management of CC.

   Metastasis is not a simple linear cascade of events but involves
   multiple parallel and intersecting pathways [55]^9^-[56]^11. The
   phenotypes of metastasis-initiating cells (MICs) depend on the
   interaction among various cell types [57]^12^-[58]^14, as well as
   modulations of multiple signaling pathways during distinct metastatic
   phases, including local invasion, dissemination, and colonization
   [59]^15. Among multiple signaling pathways involved in the process,
   transforming growth factor β (TGF-β) signaling has a crucial role in
   metastasis and cancer progression. It can not only stimulate epithelial
   mesenchymal transition (EMT) and dissemination of tumor cells [60]^15^,
   [61]^16, but also create an immunosuppressive TME via inhibiting the
   expansion and function of immune cells [62]^17. Taniguchi et al.
   reported the presence of an invasion niche in a squamous cell carcinoma
   mouse model which creates a positive feedback loop promoting tumor cell
   invasion and malignant progression via TGF-β [63]^18. Thus, further
   exploration is required to unveil the multicellular interactions and
   spatiotemporal relevance of TGF-β signaling during the progression and
   metastasis of CC.

   Single-cell RNA sequencing (scRNA-seq) has emerged as a powerful tool
   for characterizing TME and dissecting intratumor heterogeneity [64]^19.
   Several scRNA-seq studies predominantly focused on investigating the
   diversity and heterogeneity of TME in primary CC or CC initiation
   [65]^20^-[66]^24. However, how tumor cells crosstalk with other
   important cell types to reprogram the TME of CC in both temporal and
   spatial dimensions during progression remain unclear. In this study, we
   integrated scRNA-seq profiling of localized, regional, and metastatic
   TME to elucidate the dynamic alterations in the proportion and the
   biological function of different cellular components. We identified a
   TGF-β signaling-responsive malignant subset, as well as pro-tumor
   neutrophils and cancer-associated fibroblasts (CAFs) that secreted
   TGF-β. These cells were found to be enriched in metastatic TME.
   Combining scRNA-seq data with spatial transcriptomics (ST) data further
   elucidated the colocalization and cell-cell interactions between these
   subsets. In addition, we provided a clustering strategy based on bulk
   RNA-seq that related these key clusters with metastasis and patients'
   prognosis.

Materials and methods

Data collection

   The raw scRNA-seq data were obtained from our previously deposited
   database (the GSA human genome sequence archive: HRA001742 and
   HRA007492). Primary CC patients receiving no anti-tumor treatment at
   the time of sampling were included, patients with secondary CC or
   associated with other primary tumors were excluded. Detailed sources
   and clinical information of the patients in the scRNA-seq cohort are
   listed in [67]Table S1. Our cohort included three localized primary
   tumors, four regional primary tumors and four metastatic primary tumors
   of CC. The processed spatial transcriptomics data were obtained from
   Science Data Bank ([68]https://doi.org/10.57760/sciencedb.11624),
   published by Fan et al. [69]^21. Normalized The Cancer Genome Atlas
   (TCGA) expression matrix and clinical information for cervical squamous
   cell carcinoma and endocervical adenocarcinoma (CESC) cohorts were
   obtained from the cBioPortal website ([70]http://www.cbioportal.org/).

scRNA-seq data preprocessing

   Sequencing data were aligned to the GRCh38 human reference genome and
   then processed into the initial Unique Molecular Identifier (UMI)
   matrix using the Cell Ranger Single-Cell toolkit (version 7.0.1). The
   Seurat workflow (version 4.2.1) was utilized for downstream analyses
   with default parameters unless otherwise stated. Quality control
   measures were implemented to exclude cells with fewer than 1000 UMIs or
   fewer than 500 features or a mitochondrial gene fraction greater than
   25%. The NormalizeData function was then used to normalize the library
   size and log-transform the expression matrix to prepare it for further
   downstream analyses.

Data integration and cluster identification

   Two rounds of unsupervised clustering were performed to identify major
   cell types and subclusters, and to eliminate residual doublets. The
   FindVariableFeatures function was used to identify the highly variable
   genes (HVGs). After scaling the data, Principal Component Analysis
   (PCA) was performed using the RunPCA function with the top 2000 HVGs
   for dimensionality reduction. For visualization, Uniform Manifold
   Approximation and Projection (UMAP) was implemented to reduce
   dimensionality in the RunUMAP function with the top 20 principal
   components. As we observed batch effects between patients in the UMAP
   plot, we used the canonical correlation analysis (CCA) function to
   integrate expression data from different samples for batch effect
   correction.

   Cell clusters were annotated based on cluster-specific marker genes
   found by the FindAllMarkers function. Any clusters expressing more than
   two canonical cell type markers were considered as contaminated
   doublets and removed. Repeated dimension reduction and unsupervised
   clustering was performed if any doublets were removed. The first round
   of cell clustering and annotation identified major cell types including
   T and natural killer (T/NK) cells, myeloid cells, epithelial cells,
   B/plasma cells and stromal cells. To further analyze these major cell
   types, T/NK cells, myeloid cells, epithelial cells, B/plasma cells and
   stromal cells were reintegrated, re-clustered and re-annotated using
   the procedure described above.

Distributed preference analysis

   To characterize the tissue distribution of different clusters, odds
   ratios (OR) were calculated to indicate preferences as previously