Introduction

Cervical cancer (CC) is the fourth leading cause of cancer-related death among all types of women malignancies (Sung et al., 2021). The two major histological types of CC are squamous cell carcinoma of cervix (SCC) and adenocarcinoma of cervix (ADC), accounting for 70% and 25% of all cases, respectively (Cohen, Jhingran, Oaknin, & Denny, 2019). ADC patients tend to have a higher rate of rapid progression, relapse and insensitivity to chemotherapy, radiotherapy and even immunotherapy, making it a more aggressive phenotype than SCC (Sasieni, Castanon, & Cuzick, 2009). Unfortunately, the reasons are not yet well understood due to its rarity. Furthermore, the standard treatment for both types of CC is not stratified based on pathologic types, which can result in treatment failures of ADC patients. Therefore, there is an urgent need for a comprehensive understanding of ADC.

In 2022, the National Comprehensive Cancer Network (NCCN) guidelines recommended pembrolizumab (Keytruda ®, a PD1 blockade) as the first-line combination therapy for late-stage, recurrent and metastatic CC patients with positive PD-L1 expressions (Colombo et al., 2021). It is a milestone that immunotherapy has been updated as a first-line treatment for CC. While it has been reported that 30%-50% of CC patients exhibit positive PD-L1 expression, however, only an estimated 20%-30% are expected to respond positively to PD-1/PD-L1 blockade therapy (Omenai, Ajani, & Okolo, 2022; Sun et al., 2020). Therefore, CC is characterized as an immunosuppressive cancer type, due to modulation of the tumor immune microenvironment (TIME) (Cao et al., 2023). TIME, consisting of tumor cells, immune cells, cytokines and their interactions, predominantly determines the efficacy of immunotherapy (Lv et al., 2022). The mechanisms of tumor immune evasion within the TIME are similar across different types of cancer, involving the deactivation of immune surveillance pathways or recruitment of immunorepressive cells (Vinay et al., 2015). However, the specific characteristics of the TIME in ADC, which is known as insensitive to immunotherapy, remain poorly understood.

Persistent infection with human papillomavirus (HPV) is the primary cause of CC (Mann, Singh, & Kumar, 2023). Approximately 85% of CC patients are identified with HPV infection, while 5%-15% of patients test negative for HPV (X. Yang et al., 2023). It is reported that HPV-negative patients tend to experience worse clinical outcomes, including rapid metastasis to distant organs, and an increased susceptibility to developing resistance to immunotherapy (Fernandes, Viveros-Carreño, Hoegl, Ávila, & Pareja, 2022). However, current therapeutic strategies for CC do not take into account the HPV infection status as a stratification factor. Therefore, it is imperative to investigate the precise mechanisms underlying the development of HPV-negatively infected tumors. In addition, an important and unavoidable clinical issue in CC treatment is post-surgical upstaging, which affects the disease outcome and treatment approach.

To extensively investigate the cellular and molecular heterogeneity of CC, particularly ADC, the utilization of single-cell RNA sequencing (scRNA-seq) technique is crucially needed. In the current study, we focused on the remodeling of TIME in ADC and performed scRNA-seq analysis on 11 ADC samples, with another 4 SCC samples included as controls. By high-throughput analysis, we are able to characterize the landscape of TIME in ADC, with different types of tumor-promoting cells being recruited in the vicinity of the tumor, thus creating an immunosuppressive environment in ADC. The cellular crosstalk analysis provides valuable insights into the interactions between different cell types within the TIME. Furthermore, we have conducted a screening of novel signature genes as prognostic biomarkers for ADC. These biomarkers have the potential to accurately predict lymph node metastasis and may guide personalized treatment decisions more effectively.

Methods

Clinical sample collection

The 15 cases of fresh cervical cancer tissues were collected from Hunan Cancer Hospital, with informed consent obtained from 15 independent patients. The inclusion criteria are as follows: 1) The diagnosis of CC, with pathological type of ADC or SCC, should be confirmed by pathologists; 2) Pre-surgical FIGO stages should be within the range from IB (which means the tumor mass should be detectable via CT or MRI) to IIA (which means no parametrial involvement, nor distant metastasis), in which patients have explicit indication for radical surgery following NCCN guidelines; 3) No neo-adjuvant treatments should be given before the surgery, including chemotherapy and radiotherapy; 4) The surgery options should be type-C radical hysterectomy and include pelvic lymphadenectomy with/ without para-aortic lymphadenectomy. Detailed information of the 15 patients was listed in Supplementary Table 1. This study was supervised and approved by the Ethics Committee Board of Hunan Cancer Hospital.

Sample preparation for single-cell isolation

Immediately after surgical removal, the fresh tumor tissues were washed by saline for 3 times and stored in cold storage solution (MACS® Tissue Storage Solution) for transportation. To prepare high-quality samples, the tissues were cut into pieces to 1mm³ and were incubated at 37°C for 30 minutes. 0.25% trypsin solution was used to digest the tissue pieces at 37°C for 10 minutes. Then single-cell samples were filtered through a 70-μm meshed filter and suspended in a red blood cell lysis buffer (MACS® Red Blood Cell Lysis Solution) for dissolution. Before instrumental sequencing, suspension of single cells were visualizing by TC20^TM Automated Cell Counter to assess the cell viability. Each sample was confirmed to have a cell viability rate of over 80%.

Single-cell RNA sequencing

The suspension of single cells was transferred onto the 10X Chromium Single-Cell instrument to generate single-cell beads in the emulsion (GEMs). Then the scRNA-seq libraries were constructed by using the Chromium Next GEM Single Cell 3’ Reagent Kits v3 (10X Genomics, USA) and sequenced by using the sequencer Novaseq6000 (Illumina, USA). All procedures were following the manufacturer’s protocol.

Dimension reduction and cell clustering

The gene expression matrixes were generated via the Cell Ranger toolkit (version 5.0.0) of 10X Genomics platform. The package of Seurat (version 4.0.5 R) was used to normalize, convert and format data. With the function of RunPCA, the principal component analysis (PCA) was used for data dimension reduction. The functions of FindNeighbors and FindClusters were used for cell clustering. For visualization, the RunUMAP function was performed. In each unsupervised cell cluster, the gene markers were identified by the function of FindAllMarkers with comparison to other clusters. For data processing, cells with low qualities, which showed <200 expressed genes or >25% mitochondrial UMIs (unique molecular identifiers), were excluded.

Cell type annotation

The cell typing was conducted and annotated according to the selected gene markers via CellMarker database and publications (C. Li et al., 2022; Qiu et al., 2023; Xue et al., 2022). The expression of differential expressed genes (DEGs) was used to determine each cell type following the rules: Log₂Foldchange should be > 0.2 and adjusted-P values should be < 0.05.

Pseudotime trajectory analysis

To identify the developmental changes of epithelial cells, the pseudotime trajectory analyses were performed by using Monocle2 R package (version 2.22.0), in which Seurat was used as input. Calculated by the Monocle algorithm, the top 50-1000 genes were conducted for differential GeneTest. Then, the DDR-Tree and default parameters were generated to visualize the translational relationship and developmental orders among different epithelial sub-clusters.

Cell malignancy scoring

Among these epithelial cells, the malignant score of each single cell was predicted by using scCancer R package (Guo et al., 2021). The malignant cells were distinguished from non-malignant cells by inferring large-scale copy-number variations (CNVs), which was calculated by the interCNV R package as described (Guo et al., 2021; Xue et al., 2022).

Signaling pathway enrichment analysis

The biological signaling pathway of each cell cluster was identified by performing Gene Ontology (GO) analyses and Hallmark Pathway enrichment analyses, according to the Molecular Signature Database (MSigDB version 7.5.1). Then DEGs of each cell cluster were converted to ClusterProfiler (version 4.2.2) package for analysis of predicted functions.

CytoTRACE

In order to predict the differential state of cells, scRNA-seq derived data were pass to the CytoTRACE software (version 0.3.3) with R package for further analyses.

Survival analysis

The CESC dataset from TCGA database was downloaded from Xena, with clinical characteristics, survival endpoints and treatments were included. The Kaplan-Meier survival curves were generated by using R software.

Cellular communication analysis

The possible interactions, with ligand-to-receptor communicating signaling pathway, between two types of cells were analyzed using the R packages of CellChat (version 1.5.0). Seurat was used for normalization and the rank with significance was calculated. The probability of ligand-to-receptor interaction in each signaling pathway was plotted by P-value and intensity.

Immunohistochemistry (IHC)

The paraffin-embedded ADC samples (whether post-surgical or biopsy samples) were sectioned into 4μm slides. The protocol of IHC strictly followed our previous published works (Peng et al., 2019). The primary antibody of SLC26A3 was purchased from Santa Cruz Biotechnology (Cat# sc-376187) and diluted at the ratio of 1:20. The protein expression of SLC26A3 was determined using the method of H-scoring: H-score = intensity score × percentage of positivity. The intensity was scored as: Negative=0, Weak=1, Moderate=2, Strong=3. Samples were grouped as the higher expression group (H-score ≥ 150) and the lower expression (H-score < 150) group of SLC26A3.

Multiplexed immunofluorescence (multi-IF)

The multi-IF was performed using the Opal 6-Plex Manual Detection Kit (NEL811001KT, Akoya), following the manufacturer’s instructions. Generally, the slides were stained with each primary antibody sequentially following the staining cycle as: (1) antigen retrieval by citrate (pH=6) at water bath for 30 minutes; (2) independent primary antibody incubation at 4°C atmosphere overnight or 37°C atmosphere for 1 hour; (3) common HRP-crosslinked secondary antibody incubation at 37°C atmosphere for 30 minutes; (4) independent opal reactive fluorophore solution at intended wavelength (opal 570 was excluded in order to save for FOXP3 signal). Each independent primary antibody started a new cycle. After cycles finished, the slides were prepared for incubation of FOXP3 monoclonal antibody with eFluor™ 570 (Thermo Fisher, Cat#41-4777-82) conjugated, as well as DAPI (Sigma) for nucleus staining.

Statistics

Specific analyses of scRNA-seq data were decoded and visualized by using the R software. Statistics were processed using IBM-SPSS Statistics 20.0 software and R software as well, including Student’s t-test, two-sided Wilcoxon test, etc. Correlations analysis was performed by the Chi-square test. Survival analysis was performed by the Kaplan-Meier method. P values < 0.05 were considered as statistically significant.

Results

Single-cell transcriptomic landscape of cervical cancer

We collected 15 samples of cervical cancer tissues from 15 individual patients immediately after undergoing radical hysterectomy surgery. Among these cases, 11 were ADC types, with 5 being HPV-negative and 6 being HPV-positive. The remaining 4 cases were SCC types, with an even distribution of HPV-positive and HPV-negative status (Supplementary Table 1). Prior to surgery, all 15 patients were diagnosed as early FIGO stages (stage I ∼IIA) based on evaluation of CT, MRI, and physical examination. However, after surgery, 5 cases turned out to be upstaged with pathological evidence of lymph node metastasis (FIGO stage IIICP), requiring theses patients to receive radiation therapy according to NCCN guidelines. Therefore, in addition to investigate the difference between histological types and HPV infection status, we also focused on the issue of post-surgical upstaging by comparing different FIGO stages.

In order to analyze the genomic profiling, we resolved these samples at the single-cell level using the 10X genomics chromium platform. The cell viability and sample quality were ensured beforehand, and a total of 102,012 qualified cells were utilized for further analysis (Fig. 1A). The technique of Uniform Manifold Approximation and Projection (UMAP) was utilized to reduce the dimension of gene profiling data for visualization. The obtained cells were primarily categorized into 12 cell clusters, which were further re-clustered into 9 distinct cell types, based on specific gene markers (C. Li et al., 2022) (Fig. 1B-C). These cell types included T cells (38,627 cells in total, proportioning 37.87%, marked with PTPRC, CD3D, CD3E and CD3G), epithelial cells (32,199 cells in total, proportioning 31.56%, marked with EPCAM, SLP1 and CD24), neutrophils (12,586 cells in total, proportioning 12.34%, marked with CSF3R), macrophages (6,798 cells in total, proportioning 12.34%, marked with CD68 and CD163), fibroblasts (4,258 cells in total, proportioning 4.17%, marked with COL1A2 and DCN), plasma cells (4,044 cells in total, proportioning 3.96%, marked with JCHAIN), endothelial cells (1,389 cells in total, proportioning 1.36%, marked with ENG and VWF), B cells (1,326 cells in total, proportioning 1.30%, marked with CD79A and MS4A1) and mast cells (785 cells in total, proportioning 0.77%, marked with MS4A2) (Fig. 1D-E and Supplementary Table 2).

Figure 1.

According to our scRNA-seq data, ADC exhibited a significant higher proportion of neutrophils than SCC (20.05% vs. 4.23%). The proportion of T cells was also higher in ADC (41.13%) compared to SCC (35.48%). Conversely, SCC showed a nearly threefold enrichment of epithelial cell (47.43%) compared to ADC (16.48%) (Fig. 1F). When considering the HPV infection status, HPV-negative patients demonstrated decreased proportions of epithelial cells (7.78% vs. 28.38%) and neutrophils (16.66% vs. 24.68%) compared to HPV-positive patients. The proportion of T cells enriched in HPV-negative cases (45.77%) was higher than that in HPV-positive cases (32.42%). HPV-negative patients exhibited prominently higher proportions of plasma cells (10.13% vs. 1.99%) and B cells (2.84% vs. 0.77%) compared to HPV-positive patients, even though their proportions were lower than T cells, epithelial cells and neutrophils (Supplementary Fig. 1A). Furthermore, we conducted a comparison between early-stage and late-stage cases, revealing that late-stage patients have a higher proportion of epithelial cells (35.50%) than early-stage cases (28.82%) (Supplementary Fig. 1B). In summary, these data describes primarily a complex microenvironment in cervical cancers, especially in the ADC type, with varied composition of immune and tumor cells under different conditions.

The sub-clusters of epithelial cells exhibit elevated stem-like features to promote the aggressiveness of ADC

The aforehand data has shown that the majority of the CC tumor is composed of epithelial cells. In order to explore the characteristics of these epithelial cells within the TIME of ADC, we further classified them into 12 sub-clusters based on the subsets of differently expressed genes (DEGs) in each group (Fig. 2A). Each sub-cluster was named based on the most prominently up-regulated gene in each panel of signature genes (Fig. 2B-C). Several epithelial cell clusters were strongly enriched in ADC than SCC, such as Epi_04_TFF2, Epi_06_TMC5, Epi_07_CAPS and Epi_08_SCGB3A1. Notably, Epi_09_SST, Epi_10_CYSTM1, Epi_11_REG1A and Epi_12_RRAD were exclusively identified in ADC, not in SCC (Fig. 2D). The stratified enrichment of different cluster between ADC and SCC has captured our attention for further investigation. Additionally, Epi_09_SST, Epi_10_CYSTM1 and Epi_11_REG1A were solely enriched in HPV-positive cases, although no distinctive clusters had been found in HPV-negative patients (Supplementary Fig. 2A). When comparing samples from patients in the early stages with those who have lymph node metastasis, 3 clusters (Epi_02_IGLC2, Epi_05_CCL5 and Epi_12_RRAD) were slightly increased among the late-stage patients. Surprisingly, we found that the sub-cluster Epi_10_CYSTM1 was exclusively present in stage IIIC patients, making it a potential target for identification of the biomarkers for late stage patients (Supplementary Fig. 2B).

Figure 2.

We further evaluated the features of each epithelial cluster, especially their roles in regulating tumorigenesis. To assess the differentiation potential of each cell cluster, we employed the Cellular Trajectory Reconstruction Analysis using gene Counts and Expression (CytoTRACE) method (Fig. 2E). The results showed that clusters of Epi_07_CAPS, Epi_09_SST, Epi_10_CYSTM1 and Epi_12_RRAD ranked the top four in CytoTRACE score, indicating a significant potential for pluripotency and differentiation, also referred to stem-like features. Interestingly, Epi_09_SST, Epi_10_CYSTM1 and Epi_12_RRAD were exclusively enriched in ADC. The proportion of cluster Epi_07_CAPS, which acquired the highest level of stemness, was overwhelmingly higher in ADC, although Epi_07_CAPS was also identified in SCC (Fig. 2D). Subsequently, pseudotime analysis was conducted to elucidate the differentiation trajectories. The results showed that Epi_07_CAPS, Epi_09_SST and Epi_10_CYSTM1 were positioned towards the end of pseudotime developmental trajectory, which was in consistence with their high stem-like characteristics and indicated a higher degree of malignancy in ADC than SCC (Supplementary Fig. 2C). In addition, we performed the malignancy scoring analysis, which demonstrated that cluster Epi_10_CYSTM1 exhibited the highest degree of malignancy compared to other clusters (Fig. 2F). Interestingly, cluster Epi_10_CYSTM1 displayed a remarkably diverse developmental profile and was exclusively identified in ADC. Conversely, Epi_07_CAPS and Epi_12_RRAD, which exhibited high stem-like characteristics, were ranking lower in terms of malignancy scoring when compared to other clusters. This observation may partially indicate that high stemness cluster Epi_10_CYSTM1 is essential for ADC to present more aggressive features.

Next, we proceeded to further investigate cluster Epi_10_CYSTM1, particularly in the context of ADC. To this end, Gene Ontology Biological Process (GOBP) analysis on specific signaling pathways based on the gene expression pattern observed in each cell cluster was performed (Fig. 2G and Supplementary Fig. 2D). The results revealed that ADC-specific epithelial cell clusters, particularly the most malignant cluster Epi_10_CYSTM1, tended to be more active in regulating immunity-related pathways, which indicated epithelial sub-cluster might be associated with dysfunction of immunity. Notably, Epi_10_CYSTM1 cluster also played a crucial role in regulating cell-to-cell interaction (Fig. 2G). This gives us a hint that ADC cells exert their tumor-promoting functions via interacting with other types of cells, such as neutrophils or immune cells. Therefore, further investigation is warranted to explore the intricate crosstalk networks between epithelial cells and other cell types.

ADC-specific epithelial cluster-derived gene SLC26A3 is a potential prognostic marker for lymph node metastasis

Given that Epi_10_CYTSM1 represents the most aggressive cell cluster in ADC, we further examined the relevance of this cluster and the clinical features of CC. As an alternative, we also examined Epi_12_RRAD (Supplementary Fig. 3A), which was another malignant cluster predominantly found in ADC. We firstly checked the DEGs of these two clusters. In Epi_12_RRAD cluster, we identified Insulin-like growth factor 2 (IGF2) and Alcohol dehydrogenase 1C (ADH1C), which were exclusively expressed in Epi_12_RRAD, as two most specific genes of this cluster (Supplementary Fig. 3B). However, the immunohistochemistry (IHC) staining showed that IGF2 was almost negative across cases with different clinical stages (Supplementary Fig. 3C). On the other hand, ADH1C was strongly expressed in both early-stage and late-stage cases (Supplementary Fig. 3C). Likewise, the survival analysis of Epi_12_RRAD-enriched gene signatures did not demonstrate a significant difference between the two groups (Supplementary Fig. 3D). These results suggest that the Epi_12_RRAD cluster may not play a significant role in the progression of ADC.

In contrast, the panel genes enriched in Epi_10_CYSTM1 might potentially predict poor prognosis of CC patients based on TCGA database (Fig. 3C). In the cluster of Epi_10_CYSTM1, the DEG with the highest expression is CYSTM1 (Fig. 3A). However, CYSTM1 was also expressed in other cell clusters and showed low specificity (Fig. 3B). On the other hand, SLC26A3, ORM1 and ORM2 were specifically enriched in Epi_10_CYSTM1 and could be considered as potential candidate biomarkers for representation of this cluster. Nevertheless, ORM1/OMR2 were not satisfied to distinguish the severity of CC cases, due to an irrelevance between expression intensity and clinical stages (Supplementary Fig. 3E). Interestingly, the gene of SLC26A3 derived from Epi_10_CYSTM1 showed the most potential of prognostic value for CC patients (Supplementary Fig. 3E).

Figure 3.

An important and unavoidable clinical issue in CC treatment is post-surgical upstaging, which affects the disease outcome and treatment approach. In order to examine the fact of post-surgical upstaging (or misdiagnose of early-stage cases), we merged data from two independent patient cohorts obtained from two clinical centers: Xiangya hospital (Cohort 1) and Hunan cancer hospital (Cohort 2). The criterion for objective selection was outlined below (Supplementary Fig. 3F). The rates of upstaging for each center are 10.85% (23 out of 212 in Cohort 1) and 15.31% (214 out of 1398 in Cohort 2), respectively. These rates fall within the range reported by previous studies (Dabi et al., 2018; Thelissen et al., 2022). Interestingly, both cohorts of data revealed that the misdiagnose rate was significantly higher in HPV-negative patients than HPV-positive patients, especially among ADC patients. From the data of Cohort 2, it was observed that ADC patients were more likely to stage up post-surgically. In Cohort 1, we can also observed the same tendency although without significant difference between different histological types, probably due to a smaller sample size (Table. 1-2).

Table 1.

Table 2.

The identification of Epi_10_CYSTM1 as the only cell cluster found in patients with stage IIICp raises the possibility that this cluster may be a potential marker to diagnose patients with lymph node metastasis. Therefore, we detected SLC26A3 expression on biopsy and post-surgical ADC samples to examine whether it was associated with lymph node metastasis. The results from IHC staining showed that ADC patients with stage IIIC had a higher expression of SLC26A3 (Fig. 3D, E and Table 3, 4). In summary, our results propose that the Epi_10_CYSTM1 cluster could function as a diagnostic marker to predict both poor prognosis and lymph node metastasis in ADC patients. Moreover, SLC26A3 could be employed as a marker for the Epi_10_CYSTM1 cluster, aiding in the diagnosis of lymph node metastasis to prevent post-surgical upstaging in ADC patients in the future.

Table 3.

Table 4.

Enrichment of regulatory T cells indicates an immunosuppressive status of ADC

T cells play a crucial role as the primary effectors of tumor immunity. Our data reveals that T cells constitute the predominant cell population within the TIME of cervical cancer (Fig. 1A). A total of 52,297 T cells were categorized into 8 main clusters (Fig. 4A). The gene expression patterns of T cells in CC were in accord with those observed in other types of cancers (Sorin et al., 2023; Xue et al., 2022). Based on canonical genes (C. Li et al., 2022), these cells were re-clustered as the following 4 types of T cells: exhausted T cells (marked with CD8, TIGIT, PDCD1, LAG3, HAVCR2), cytotoxic T cells (marked with CD8 and CD3), regulatory T cells (Treg, marked with CD4, FOXP3 and IL2RA), activated T cells (marked with CD8, C69 and CD3G) (Fig. 4A-C). The gene subset of cytotoxic T cells exhibited a higher expression of KRT17, KRT19, while Tregs showed a higher expression of genes such as CCL4, GZMB and CCL3 (Fig. 4D). When comparing the checkpoint pathway states among different T cell clusters, we observed relatively higher levels of LAG3 and PDCD1 (PD1) in exhausted T cells, while CTLA4 and TIG1 were more highly expressed in Tregs. These genes were identified to present immunosuppressive functions (Fig. 4E). On the other hand, certain genes, such as NKG7, PRF1, INFG and TNFRSF, which indicate the activation of immune checkpoint pathways, exhibited low expression levels in cytotoxic T cells and activated T cells (Fig. 4E). These findings lead us to hypothesize that the establishment of an immunosuppressive TIME in ADC may involve the recruitment of regulatory T cells (Tregs) to the tumor area and the subsequent inactivation of a substantial proportion of cytotoxic T cells.

Figure 4.

We further compared different status of CC with different attribution of T cell subtypes. It was demonstrated that ADC had a significantly higher proportion of Tregs compared to SCC (Fig. 4F). Meanwhile, ADC patients who tested negative for HPV showed a higher enrichment of Tregs and a lower proportion of activated T cells compared to those with HPV-positive status (Supplementary Fig. 4A). These findings shed light on the underlying mechanism for immunotherapy insensitivity of ADC, especially in HPV-negative status. Furthermore, analyses of GOBP have shown that Tregs were more active in regulating immune activities in ADC than SCC, indicating a greater tendency for Tregs to be activated in the TIME of ADC (Fig. 4G). In ADC, the proportion of activated T cells was higher compared to SCC. However, the GOBP analysis suggested that their leading roles were to promote the apoptosis of immune cells than to exert cytotoxicity (Supplementary Fig. 4D). In addition, in HPV-negative cases, activated T cells demonstrated a stronger capacity to interact with other cell types, including neutrophils and epithelial cells (Supplementary Fig. 4E). Consequently, it is further proved that the enrichment of Tregs in ADC does execute immunosuppressive functions.

Tumor associated neutrophils (TANs) surrounding ADC tumor area contribute to the formation of a malignant microenvironment

The large-scale enrichment of specific gene subsets in neutrophils in CC prompted us to further examine their roles. A total of 12,586 neutrophils were detected and the composition proportion in ADC (20.05%) was almost five times higher than that in SCC (4.23%) (Fig. 1F). Based on this observation, we hypothesize that ADC-specific neutrophil clusters may be related to poor prognosis and increased malignancy. All the cells were initially grouped into 6 clusters based on the gene expression modules (Fig. 5A). In the context of cancer, tumor associated neutrophils (TANs) can exhibit dual functions, either promoting tumor (pro-tumor TANs) or inhibiting tumor progression (anti-tumor TANs). Based on the consensus markers identified by Jaillon et al. (Jaillon et al., 2020) and Xue et al.(Xue et al., 2022), we performed further clustering of the neutrophils in our dataset, resulting in the identification of 3 sub-types of TANs: pro-tumor TANs (Neu_03, marked with CD66, CD11 and MT1X), anti-tumor TANs (Neu_02, marked with CD66, CD11, CCL5, KRT5 and ELF3), TANs with isg (interferon stimulated genes) and other undefined types (Neu_01, distinguished by CD66, CD11, IFITM2, S100A8 and LST1) (Fig. 5B-C). Generally, the Neu_01 sub-cluster, which consists of TANs with isg with the subset of specific genes, accounted for the largest proportion among all neutrophils. Notably, the presence of pro-tumor TANs was significantly higher in ADC compared to SCC (Fig. 5D), as well as in HPV-negative cases compared to HPV-positive cases (Supplementary Fig. 5A). On the other hand, the proportion of anti-tumor TANs in ADC was predominantly lower than that in SCC (Fig. 5D).

Figure 5.

We further conducted the GOBP analysis on gene subsets derived from pro-tumor TANs to explore their potential functions. The findings revealed that pro-tumor TANs in ADC were more active in regulating cell-to-cell interactions and chemotaxis. Specifically, they were involved in pathways such as chemokine-mediated or cytokine-mediated signaling, neutrophil activation, and receptor internalization (Fig. 5E). These results suggest that pro-tumor TANs may engage in extensive communication with various cell types. Interestingly, pro-tumor TANs exhibited enhanced activation in mediating inhibitory pathways related to immunity (Fig. 5E). These pathways include the negative regulation of leukocyte proliferation and lymphocyte activation, as well as the suppression of IL-10 production, which is known to promote immune responses (Saraiva, Vieira, & O’Garra, 2020; Sarvaria, Madrigal, & Saudemont, 2017).

The role of pro-tumor TANs in association with disease prognosis was further verified using TCGA database. The patients with higher expression of gene panels specific to pro-tumor TANs exhibited a poorer prognosis compared to those with lower expression (Fig. 5F). These findings suggest that the increased malignancy observed in ADC, particularly in HPV-negative cases, may also be attributed to the enrichment of pro-tumor TANs and evasion of anti-tumor TANs.

Cellular heterogeneity of plasma/B cells in ADC

In our study, we identified 6 distinctive clusters of plasma/B cells based on the set of gene enrichment: Plasma/B_01(marked with IGHA2, IGHG2 and IGLC2), Plasma/B_02 (marked with HLA-DRA, HLA-DPB and CD37), Plasma/B_03 (marked with KRT17 and S100A2), Plasma/B_04 (marked with CCL5, CXCL13, IL-32 and CCL4), Plasma/B_05 (marked with CXCL8), Plasma/B_06 (marked with HMGB2) (Fig. 6A-C). Among these clusters, the most abundant one was Plasma/B_01_IGHA2 (Fig. 6D and Supplementary Fig. 6A-B). Interestingly, the proportion of Plasma/B_01_IGHA2 was found to be significantly higher in ADC compared to SCC, and it was also slightly higher in HPV-negative cases than in HPV-positive cases (Fig. 6D and Supplementary Fig. 6A). Generally, B-lymphocytes are known to have tumor-inhibitory properties. However, the ectopic gene module of plasma/B cells in ADC suggested that some sub-clusters might present a tumor-promoting role. This hypothesis was further supported by the results of GOBP analysis, which indicated that Plasma/B_01_IGHA2 gene signatures were more active in negatively regulating B cell activation, complement activation and receptor signaling pathway (Fig. 6E). When we combined the top 50 signature genes in Plasma/B_01_IGHA2-enriched subset and performed the survival analysis using TCGA database, we observed that higher expression of the signature genes predicted a poorer prognosis of CC patients (Fig. 6F). Taken together, the cellular heterogeneity of plasma/B cells in the TIME of ADC is well identified at the single-cell level.

Figure 6.

Crosstalk among tumor cells, Tregs and neutrophils establishes the immunosuppressive TIME in ADC

In the above description, we have outlined that certain ADC-enriched cell sub-clusters, such as Epi_10_CYSTM1, Tregs, and pro-tumor TANs, played oncogenic roles in regulating various cellular activities, including maintenance of stemness and evasion of immune surveillance. To gain deeper insights into the interplay among these clusters and their influence on ADC progression, we harnessed the Cellchat tool (Jin et al., 2021). Our primary focus revolved around unraveling the mechanisms governing Treg recruitment to the tumor microenvironment and the maintenance of stemness within ADC-enriched cell sub-clusters.

Given that the 3 ADC-specific sub-clusters (Epi_07_CAPS, Epi_10_CYSTM1 and Epi_12_RRAD) of epithelial cells exhibited quite high levels of stem-like properties, and particularly that Epi_10_CYSTM1 presented the highest degree of malignancy among them, we further analyzed their interactions with other types of cells. Firstly, in comparison to SCC, ADC-enriched epithelial sub-clusters exhibited a higher propensity for interaction with Tregs through ligand-receptor signaling pathways, including ALCAM (ALCAM>>CD6) and MHC-II (HLA-DR>>CD4) (Fig. 7A-E), which have been reported to be essential for Treg recruitment, expansion and stabilization to establish the immunosuppressive microenvironment (Chalmers et al., 2022; Ferragut, Vachetta, Troncoso, Rabinovich, & Elola, 2021; Freitas et al., 2019). To further uncover the underlying reason for the enhanced stemness observed in ADC epithelial cells, we focused on the signaling received by Epi_10_CYSTM1 cell cluster. Interestingly, the Cellchat analysis revealed that, compared to other pathways, such as CD46>>JAG1 and GZMA>>PARD3, the interaction of TGFβ>>TGFBR between Tregs and epithelial cells (particularly in cluster Epi_10_CYSTM1) was remarkedly activated in ADC but not in SCC (Fig. 7F-G). TGFβ is a canonical secreted protein that induces carcinogenesis, such as cancer cell proliferation, invasion and self-renewal (Derynck, Turley, & Akhurst, 2021; Massagué, 2008; Tauriello, Sancho, & Batlle, 2022). Thus, we speculate that the recruited Tregs might secrete TGFβ to stimulate the tumor epithelial cells of ADC, leading to increased stemness. Additionally, pro-tumor TANs were also recruited to the tumor area and interacted with Epi_10_CYSTM1, as well as Epi_07_CAPS and Epi_12_RRAD cells, via the ANNEXIN-to-FPR1/FPR2 signaling pathway (Supplementary Fig. 7A-B). Previous studies have reported that the ANNEXIN pathway can promote the chemotaxis of neutrophils (Araújo et al., 2021). Likewise, TANs expressing isg, which have been shown to promote tumorigenesis as mentioned earlier, were also recruited to the tumor area through the ANNEXIN pathway. These findings suggest that the aggressive clusters of TANs may contribute to the formation of an ADC-specific TIME.

Figure 7.

To approve these hypotheses, we conducted immunofluorescence (IF) analysis and observed a greater recruitment of FOXP3⁺ cells in close proximity to regions exhibiting high SLC26A3 expression (Epi_10_CYSTM1 cell cluster) (Fig. 7H). Furthermore, when compared to regions with low SLC26A3 expression, tumors displaying high SLC26A3 levels exhibited significantly increased expression of Aldehyde Dehydrogenase 1 Family Member A1 (ALDH1A1), a known marker for cancer stem cells in CC (Douville, Beaulieu, & Balicki, 2009) (Supplementary Fig. 7C). Notably, we identified SLC26A3^high tumor cells expressing ALCAM in close proximity to CD6+ Tregs, suggesting potential involvement of ALCAM>>CD6 signaling in Treg recruitment (Fig. 7I). Additionally, we observed reduced E-cadherin expression in tumor cells within SLC26A3high regions, possibly induced by TGFβ secreted by adjacent Tregs (Fig. 7J and Supplementary Fig. 7D). Altogether, these findings confirm the intricate crosstalk between the Epi_10_CYSTM1 cell cluster and Tregs, which plays a pivotal role in establishing an immunosuppressive TIME and sustaining heightened stemness in tumor cells.

Discussion

Cervical cancer is often regarded as less aggressive compared to ovarian cancer and endometrial cancer, primarily due to the effectiveness of HPV vaccines in prevention. When diagnosed at early stages, the 5-year survival rate for cervical cancer patients exceeds 90%, surpassing that of other gynecologic malignancies (Gennigens et al., 2022; Schiffman, Castle, Jeronimo, Rodriguez, & Wacholder, 2007). However, if cervical cancer metastasizes to the intraperitoneal or retroperitoneal lymph nodes, the survival rate will decrease to 30%∼60% (Cohen et al., 2019). Worse still, the 5-year survival rate for recurrent CC patients is less than 20% (H. Li, Wu, & Cheng, 2016). Specifically, when considering the ADC type alone, the prognosis is even worse. Unfortunately, our understanding for ADC is limited. Most previous studies have concentrated on SCC type of CC, and obtained heterogenetic information derived from bulk RNA-seq method. Therefore, it is worthwhile to shift our attention towards an in-depth investigation on the disparities between ADC and SCC at the single-cell resolution.

More recently, researchers have started to utilize scRNA-seq approach to investigate the tumor microenvironment of CC. For instance, Cao et al compared 5 cases of SCC tumor tissues with paired adjacent normal tissues and described the immune landscape in which specific clusters of T and B cells exhibited immune-exhausting or activating processes (Cao et al., 2023). Another study conducted by Qiu et al aimed to elucidate the distinct molecular patterns of immune reactions between ADC and SCC in the context of different HPV infection status at the single-cell level (Qiu et al., 2023). However, due to limited samples tested in these studies, our knowledge on the TIME in CC is yet inadequate. In the current study, we have generated a genomic atlas of ADC, which may depict the landscape of TIME built by both ADC tumor cells and neighboring cells. Apart from the predominant presence of epithelial cells, the cell clusters in ADC are majorly composed of T cells, tumor-associated neutrophils (TANs), and plasma/B cells. The gene subsets identified within ADC-abundant cell sub-clusters, including Epi_10_CYSTM1, Tregs, pro-tumor TANs and P/B_01_IGHA2, exert potential roles in promoting both the carcinogenesis and immunosuppression of ADC. In the TIME of ADC, the unique module of cellular interaction has been decoded and validated that the immunosuppressive Tregs and the tumor-promoting TANs are recruited to tumor cells via chemo-attraction. Finally, our hypothesis from scRNA-seq data is further confirmed using clinical samples, providing evidence that SLC26A3 is a reliable biomarker for ADC (Supplementary Fig. 8).

It has been reported that immunotherapies are less effective on ADC patients, especially those who are not infected with HPV (Attademo et al., 2020; Ferrall, Lin, Roden, Hung, & Wu, 2021). As the two major components of ADC cell types, T cells and epithelial cells have been shown to be the dominant regulators in the TIME due to their reciprocal interaction dynamics. We identified several ADC-specific epithelial cell clusters that exhibit strong stemness and aggressiveness. Among these clusters, Epi_10_CYSTM1 demonstrated the highest level of malignancy, but the underlying mechanism behind this behavior remains unknown. The Cellchat analysis revealed that Tregs might quite actively crosstalk to Epi_10_CYSTM1. To confirm this, immunostaining for co-localization demonstrated an increased aggregation of Tregs around tumor cells exhibiting overexpression of SLC26A3. It has also been widely reported that tumor cells acquire higher stemness present more aggressive features of resistance to immunotherapy (Dianat-Moghadam et al., 2022). This study confirmed that TGFβ produced by Tregs induce the stemness of Epi_10_CYSTM1 cells by activating downstream oncogenic pathways, inducing EMT on tumor epithelial cells (Nixon, Gao, Wang, & Li, 2023). Interestingly, Epi_10_CYSTM1 cells express ALCAM, which is a soluble ligand, and bind to its receptor CD6 on Tregs, activating chemotaxis and facilitating Treg recruitment. The infiltration of a large number of Tregs in ADC significantly accelerates immune evasion (Supplementary Fig. 8). Therefore, the presence of a higher number of Tregs infiltrating the tumor area of ADC may explain the increased resistance to immunotherapy. It is worth noting that, in comparison to the proportions of Tregs, the majority of T cell fractions in the CC tumor area are comprised of cytotoxic T cells and exhausted T cells. However, their distribution is relatively equal, indicating the dysfunctional status of the majority of infiltrated T cells in the CC tumor area.

Another interesting phenomenon observed in our sing-cell sequencing data is the involvement of TANs in the regulatory interactions within the TIME of ADC. It has been clearly investigated that pro-tumor TANs, characterized by CD66, CD11, CD170 and PD-L1, exert the oncogenic roles by promoting tumor cell proliferation, angiogenesis, inducing genetic instability, and most importantly causing immunosuppression (Jaillon et al., 2020). In this study, pro-tumor TANs are predominantly enriched in ADC tissues compared to SCC, especially in HPV-negative ADC cases. This suggests that ADC-specific epithelial cell sub-clusters actively engage in the recruitment of pro-tumor TANs through the ANXA1-FPR1/FPR2 ligand-receptor interaction.

The anti-tumor function of B cells has been identified previously (Bod et al., 2023; Helmink et al., 2020). However, our study has presented a specific cluster of plasma/B cells with tumor-promoting functions. The overexpression of gene signatures from Plasma/B_01_IGH2 cells is associated with poor regression in CC patients. Interestingly, Plasma/B_01-IGHA2 was found to be closely related to the PC Cluster (Plasma cell Cluster) in Cao’s study due to the shared gene set (Cao et al., 2023), such as IGHA1, IGHA2, IGHG2, IGHG4, etc. However, their roles were contradictory. According to that study, PC cluster was derived from the local expansion of TIME and showed potential effectiveness in anti-tumor immunity. In our case, Plasma/B_01-IGHA2 had a reverse effect and played a tumor-promoting role in ADC. This difference in function is likely attributed to the distinct cell types and their interactions in ADC, instead of SCC. Although the cell-to-cell interaction did not yield positive results, probably owing to the small cell proportion, the underneath dynamics of B cells still need further investigation.

Most importantly, in this study SLC26A3 is identified as a promising diagnostic biomarker to predict lymph node metastasis in ADC patients at stage IIIC. In 2018, the International Federation of Gynecology and Obstetrics (FIGO) updated a new criterion of stage IIIC by including lymph node status, which indicates that the therapeutic strategy is completely different from early stages. There are two stratification of stage IIIC: IIICp (p stands for pathological findings of LN, which means misdiagnosis) and IIICr (r stands for radiological findings of LN, which means pre-surgical evidence) (Bhatla, Aoki, Sharma, & Sankaranarayanan, 2018). Unfortunately, radiological imaging tools do have some limitations, as CT or MRI can only detect suspiciously metastatic LN when the minimum diameter is ≥ 10mm (Saleh et al., 2020). However, over 10% of CC patients who are initially diagnosed with early-stage and deemed suitable for radical surgery, end up being upstaged due to lymph node metastasis discovered after the surgery (Dabi et al., 2018; Thelissen et al., 2022). For these patients, they require additional radiotherapy such as external beam irradiation besides surgery (Abu-Rustum et al., 2020). However, if the presence of lymph node metastasis could have been accurately diagnosed in biopsy sample, the optimal therapeutic approach for these patients would be radiotherapy alone. Moreover, combination of both surgery and radiotherapy may lead to more complications for patients than radiotherapy alone, without offering additional benefits (Cushman et al., 2018; Kashima et al., 2023; J. Yang et al., 2015). To solve this problem, a reliable biomarker is required to diagnose LN metastasis in conjunction with radiographic tool. Combining two large cohorts of clinical data, we identified that CC patients with ADC types were more vulnerable to be misdiagnosed for clinical staging. Derived from scRNA-seq data, SLC26A3 could be employed as a marker for the Epi_10_CYSTM1 cluster, which was found to be exclusively enriched in stage IIIC cases of ADC. The IHC staining further confirmed that ADC patients at stage IIIC exhibited higher expression of SLC26A3 compared to early stages (FIGO stage I∼ II). Furthermore, we tested on biopsy samples and it implied that SLC26A3 could be a significant predictor of LN metastasis. Our study for the first time provides evidence that ADC patients are at a higher risk of developing lymph node metastasis, which can be challenging to detect solely through radiological imaging techniques. However, by incorporating lymph node imaging with biomarkers such as SLC26A3, it may significantly reduce the misdiagnosis rate among FIGO IIICp patients.

In summary, we characterized the genomic landscape of the TIME in ADC at a single-cell resolution. Our study revealed the presence of specific epithelial cell clusters in ADC that exhibited high aggressiveness. The enriched gene signatures within these clusters may serve as potential therapeutic targets to improve the efficacy of immunotherapy. The interactions between tumor epithelial cells and other types of cells, such as T cells and neutrophils, comprehensively depict a highly immunosuppressive microenvironment of ADC. Moreover, our investigation extended to the issue of upstaging at the single-cell level, with the identification of potential biomarkers from specific tumor cell clusters showing promise as diagnostic indicators for the precise diagnosis of cervical cancer patients.

Additional information

Author contributions

Yang Peng, Conceptualization, Methodology, Formal analysis, Investigation, Data curation, Visualization, Writing – original draft preparation, Project administration, Funding acquisition; Yilin Li, Methodology, Validation, Investigation, Resources, Data curation, Visualization; Jixing Ao, Methodology, Validation, Software, Resources, Investigation, Data curation; Jia Shen, Software, Validation, Investigation, Data curation, Visualization; Xiang He, Software, Investigation, Data curation, Visualization; Dihong Tang, Methodology, Formal analysis, Writing – review & editing; Chaonan Chu, Formal analysis, Writing – review & editing; Congrong Liu, Conceptualization, Methodology, Formal analysis, Writing – review & editing, Supervision, Project administration, Funding acquisition; Liang Weng, Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Writing – original draft preparation, Writing – review & editing, Supervision, Project administration.

Ethics approval and consent to participate

Specimen collection and pathological data used in this study were approved by the Ethics Committee of Hunan Cancer Hospital (No. KYJJ-2020-036). Written informed consent was obtained from patients for use of samples. This study was also in accordance with the Declaration of Helsinki and none of any procedures conducted interfered with the treatment plan of patients.

Data availability

In this study, the single-cell RNA sequencing data generated is available in the Genome Sequence Archive for Human (GSA-Human, from China National Center for Bioinformation) by accession number: PRJCA023028 (https://ngdc.cncb.ac.cn/gsa-human/). These original sequencing data can be obtained only for non-profitable use, under request from the corresponding authors.

Competing interests

All authors declare no competing interests of publication.