A CSCC cellular map from expression and genetic features by scRNA-seq
To comprehensively interrogate the CSCC ecosystem, we first conducted scRNA-seq of tumors from 14 treatment-naïve patients and 3 normal cervix samples obtained from healthy donors (Figure 1A, Table S1). After undergoing quality control and doublet removal, we obtained a total of 163,880 cells and constructed a cell type map which comprised 50 minor cell subtypes derived from 14 major cell types annotated using canonical marker genes, including CD4 T cells (6 subtypes), CD8 T cells (7 subtypes), NK cells (5 subtypes), B cells (2 subtypes), plasma B cells (2 subtypes), macrophages (7 subtypes), DC cells (5 subtypes), mast cells (1 subtype), neutrophil cells (1 subtype), fibroblasts (7 subtypes), myocytes (3 subtypes), endothelial cells (3 subtypes), and glandular cells (1 subtype), and squamous (1 subtype) (Figures 1B, and S1A-B).
In accordance with previous findings13–15,several cell types were selectively expression in the CSCC samples, including exhausted T cells13 (CD4_Th-CXCL13, CD8_Tex-GZMK, CD8_Tex-PDCD1), myeloid cells14 (Macro-SPP1, Macro-ISG15, cDC2-CD1A, DC-LAMP3, Mast-KIT), and cancer associated fibroblasts15 (CAF) (CAF-DES and CAF-MMP11) (Figure S1C).
Consistent with the CSCC pathological diagnosis, epithelial populations were mainly composed of squamous rather than glandular cells in all tumors (Figure 1B). To further distinguish malignant cells from normal squamous cells, we inferred large-scale chromosomal copy-number variations (CNVs) based on averaged expression profiles across chromosomal intervals in each squamous cell with normal epithelial cells from healthy donors (HDs) being used as control. These inferred CNVs, which were consistent with their corresponding WES CNVs (Figures S1D), effectively separated malignant from normal squamous cells. CNVs profiles and mutational landscape were consistent with prior studies16 (Figures S1D-E), supporting the quality of the scRNA-Seq data. Of note, a subset of tumors, such as CC05, displayed significant intra-tumoral heterogeneity in their CNV patterns at the single cell level (Figure S1F).
Patterns of expression heterogeneity within the malignant compartment of CSCC
Non-negative matrix factorization (NMF) provides a powerful tool for uncovering intra- and inter-tumoral heterogeneity from scRNA-seq data that contributes to tumor metastasis17,18 and drug resistance15. To explore transcriptional states of the CSCC malignant squamous cells, we applied NMF to each malignant cell (see Methods). First, we detected 79 robust expression programs across all analyzed patients (Table S2), then hierarchical clustering was performed to distill these 79 programs into 8 meta-programs (MPs) that reflect common patterns of intra-tumoral expression heterogeneity among different tumors (Figure 1C). Each MP was functionally characterized through enrichment analysis of signature genes (Table S3). Notably, our CSCC NMF programs were highly consistent (albeit not identical) with those previously observed in pan-cancer15 and head and neck squamous cell carcinoma (HNSCC)18 (Figure S1G).
MP1 and MP2 that resembled stress response programs identified in pan-cancer and HNSCC datasets (Figure S1G) and were enriched for oxidative phosphorylation (MT-ND4, MT-ATP6, MT-CO3) and immediate early genes implicated in cellular stress responses (FOSB, ATF3, JUN), respectively, were designated as oxidative phosphorylation and stress response MPs respectively (Figure 1C-D). MP3 that was concordant the pan-cancer and HNSCC partial EMT program (pEMT)18 (Figure S1G), consisted of genes associated with ECM (CAV1, COL17A1, LAMB3 and LAMA3) and had features of the pEMT state (ITGA6, TGFBI) (Figures 1C-D) was named MP pEMT. MP4 and MP5 that strongly resemble HNSCC proliferation programs, were characterized by enrichment of cell cycle related genes corresponded to S phase (PCLAF, RRM2, TYMS) and mitotic phase (UBE2C, TOP2A, AURKB) respectively were named MP DNA replication and MP mitosis.
It is noteworthy that three additional meta-programs (MP6, MP7, and MP8), with signature genes primarily composed of epithelial genes such as keratins (KRT5, KRT6B, KRT14, and KRT16), mucins (MUC1, MUC4, MUC16, and MUC22), and small proline-rich proteins (SPRR2A, SPRR2D, and SPRR3), were identified (Figures 1C-D). The expression of genes from these three MPs varied coherently across malignant cells (Figure S1H).
MP6 (Epi-Krt) that demonstrated high concordance with the HNSCC squamous epithelial differentiation (Epi. Diff. 2) signature18, was characterized by enrichment of keratinization related genes (KRT5, KRT6B, KRT14, KRT16) consistent with terminal differentiation of cervix keratinizing epithelia (Figures 1D, and S1G). MP8 (Epi-Sen (senescence)) was similar to the HNSCC squamous epithelial intermediate differentiation (Epi. Diff. 1) state18 and displayed features of cell senescence (Figure 1D, and S1G). Finally, MP7 (Epi-Imm) exhibited high correlation with pan-cancer immune response program15 containing multiple immune associated molecules (HLA-DRB, B2M, and IDO1). MP7 also displayed co-expression of squamous epithelial markers, such as KRT5 and KRT6, together with glandular markers LCN2, MUC1, and WFDC2, a subtype which has not been reported previously. The presence of three distinct epithelial programs may reflect epithelial differentiation states that could have implications for patient prognosis and treatment strategies.
Epi-Krt and Epi-imm meta-programs are inversely and directly correlated with immune cell infiltration, respectively
Correlation analysis demonstrated positive associations between MP4 and MP5 (Pearsons’ r = 0.32, p < 0.001), MP1 and MP2 (Pearsons’ r = 0.34, p < 0.001), MP7 and MP8 (Pearsons’ r = 0.36, p < 0.001) (Figure 1E). A negative correlation between MP3 (pEMT) and MP6 (Epi-Krt) (Pearsons’ r = -0.43, p < 0.001) recapitulated patterns observed in HNSCC in which the pEMT state was linked to metastasis18. Notably, MP6 was also negatively associated with MP7 (Pearsons’ r = -0.40, p < 0.001) (Figure S1I).
To prioritize the most relevant MPs associated with the TIME, and to determine the directionality of their correlation, we correlated the summarized expression (eigengene) values of each MP with percentages of TIME cells inferred from scRNA-seq data within individual patients. Interestingly, MP6 was positively associated with fibroblast cell (CAF-DES and CAF-MMP11) content and negatively with immune cell (CD4_Th-CXCL13, CD8_Tex-GZMK and CD4-MKI67) content (Figure 1E, Table S4), consistent with an association with immune exclusion. Conversely, MP7 was positively associated with immune cells and negatively with stromal cells (Figure S1J, Table S4), consistent with a role in immune infiltration. These findings highlight the opposite roles of MP6 and MP7 in terms of immune cell infiltration and exclusion.
To better characterize heterogeneity and plasticity of malignant cells in CSCC, we employed an unsupervised Leiden-based clustering after correcting batch effect across patients with the bbknn method19. Uniform manifold approximation and projection (UMAP) identified five unique squamous cell states that were all detected in each tumor (Figure 1F). To established the relationships between identified clusters and meta-programs, MP activities were quantified (see Methods) and then projected onto UMAP space (Figure S1K). While the MP1 and MP2 programs were present in all tumor cells, the MP3, MP4, MP5, MP6, MP7 and MP8 programs were present in different tumor cell subpopulations in UMAP space (Figure S1K). We subsequently annotated tumor subclusters based on corresponding enriched MPs (Figure 1F).
To examine the distinctions and similarities between normal cervical epithelial cells and tumor cells, we re-clustered normal cervical epithelial cells from HDs resulting in three sub-populations each expressing known representative genes, including (1) basal cells (TP63), (2) glandular cells (MUC5B), and (3) keratinized cells (KRT6A) (Figure 1G and S1L)20. Although exhibiting shared marker genes (KRT5, KRT6A, S100A8, S100A9), the MP6 (Epi-KRT) tumor subpopulation demonstrated higher levels of these genes when compared with their normal counterparts (HD-Squamous cells), consistent with expression of these genes changing with transformation from squamous cell to tumor (Figure 1H and S1M). Interestingly, MP7 (Epi-immune) tumor cells that co-expressed squamous cell markers (KRT5, KRT6A) as well as partial glandular-specific genes (such as LCN2, MUC1 and WFDC2, but not MUC5B, and AGR220) did not have a normal cervical counterpart (Figure 1H). Similar to the MP7 tumor cell type, MP8 (Epi-Sen) did not have an equivalent in the normal cervix consistent with this phenotype being uniquely acquired during transformation to the malignant state.
Epi-Imm is associated with abundant tumor infiltrating T cells skewed toward a cytotoxic/exhaustion state
Beside immune cell abundance, the spatial distribution of immune infiltrates, especially CD8+ T, has been shown to be important for TIME stratification and response to immune therapy21. To assess the relationship between MPs and TIP5, we performed CD8 immunohistochemistry staining (IHC) analysis to assign a TIP to each tumor evaluated by scRNA-seq based on the location of CD8+ T cells in the TIME: (1) an infiltrated phenotype in which CD8+ T cells infiltrate the tumor epithelium, (2) an immune excluded phenotype with CD8+ immune cells primarily embedded in the surrounding tumor stroma, and (3) an immune desert phenotype represented by tumors devoid of CD8+ immune cells (Figure 2A, Table S5)5.
Based on the results, we categorized tumors into infiltrated or excluded/desert TIPs (Figure 2B). CSCCs exhibiting an infiltrated phenotype had a higher number of CD8+ T cells compared to excluded/desert tumors, as inferred from scRNA-seq analysis (Figures 2C-D). The variation of immune cells in different TIPs stimulated deeper analyses to assess immune cell phenotypic states. Density estimation and immune cell enrichment based on the ratio of observed to expected cell numbers (Ro/e) highlighted clear TIP-specific differences, especially in T cell states (Figure 2E-G, and S2A). In particular, CD8_Tex-GZMK, CD8_Tex-PDCD1, and CD8-MKI67 cells, which are expected to exhibit cytotoxic capacity but co-expressed inhibitory receptors and exhaustion markers, were enriched in infiltrated CSCCs and depleted in excluded/desert tumors (Figure S2B), potentially indicative of antitumor immune response along with immune evasion feedback occurring in infiltrated CSCCs but not in excluded/desert CSCCs.
Importantly, for tumor cell states based on scRNA-seq clustering, we found the excluded/desert CSCCs demonstrated higher proportion of MP6 cells and the infiltrated CSCC had a higher proportion of MP7 cells (Figures 2H-I), demonstrating a relationship of tumor meta-programs with TIP.
Spatial transcriptomics confirmed inverse associations of Epi-Krt and Epi-Imm with immune cell infiltration
scRNA-seq results in a loss of spatial context preventing parsing of spatial organization and interactions of neighborhood cells22. To further evaluate the relationship between Epi-Krt/Epi-Imm and immune cell infiltration, we assessed the spatial organization of tumor Epi-Krt/Epi-Imm populations and their neighboring TIME cellular composition using spatial transcriptomics (ST) data from 15 CSCCs with stereo-seq, a method with high resolution and sensitivity in our previous published study 23.
Integrating ST data with scRNA-seq data, we performed probabilistic label transfer (see Methods) and annotated each spot to specific cell types with the highest cell type score (Figures 3A-B, S3A). Impressively, cellular architecture showed high consistency with CSCC histologic architecture assessed by HE analyses by a pathologist. To compare the spatial heterogeneity of tumor Epi-Krt/Epi-Imm states, we scored tumor spots with MP6 and MP7 signature genes and applied hierarchical clustering to classify tumor spots as MP6 or MP7 in each section individually to decrease technical bias across patient tumors (Figure 3C, S3B-C). The distribution of MP6 and MP7 spots varied across different slides. For instance, while some slides demonstrated regional MP6 and MP7 areas, slides from patients TJH14 and TJH90 were dominated by a MP6 region, while slides from patients TJH17 and TJH48 were dominated by a MP7 region, indicating that the two states were spatially restricted and differentially present in individual tumors (Figure S3D-E).
To study the composition of TIME cells associated with MP6 and MP7 regions, we defined a cell neighborhood score by calculating the average cell type score surrounding corresponding spots (Figure 3D-E, see Methods). In all 15 ST slides, MP7 spots demonstrated elevated immune cell neighborhood scores, including CD4 T, CD8 T, B, Macro and DC cells, when compared to MP6 spots (Figure 3F and S3F). Moreover, correlation analysis of the MP6, MP7 and cell neighborhood scores illustrated a clear positive correlation of MP7 with immune cells neighborhoods, while the MP6 (only positively correlated with tumor neighborhood scores) demonstrated a significant negative correlation with immune cells neighborhoods (Figure 3G), further supporting the mutually exclusive association between MP6 and immune infiltration and the coordinate association between MP7 and immune cells.
Next, to further explore the relationship between tumor state and immune infiltration histologically, we examined the expression of genes chosen to distinguish various tumor cell states and T cell infiltration. Consistently, tumors with high SERPINB3 and FABP5 (MP6 dominant tumors) showed sparse staining with MUC4 and, importantly, lacked CD8 T cell infiltration (Figure 3H). In contrast, tumors with high MUC4 (MP7 dominant tumors) showed fewer SERPINB3 and FABP5-positive tumor cells but abundant CD8 T cell infiltration (Figure 3H).
CNV diversity is not sufficient to define phenotypic diversity in CSCC
The CNV diversity has been demonstrated to be one driver of tumor heterogeneity24. As described above, considerable inter- and intra-tumoral heterogeneity exists in the single-cell CNV patterns. We then compared the assignment of cells to CNV-based genomic subclones with their patterns of transcriptional heterogeneity. Surprisingly, cell clones defined by CNV patterns did not directly align with expression subtypes (See Figure S1F). In parallel, cells harboring highly similar CNV patterns harbor multiple distinct transcriptional states (See Figures S1F), indicating the phenotypic heterogeneity in the CSCC as defined by MP subtype is not simply a result of genetic heterogeneity.
Similarly, in a single patient slide, there was a spatial distribution of a single malignant clone as defined by CNVs inferred from ST data (Figure S3G). However, when MP6 and MP7 states were projected onto the malignant cells in the slide, there was no clear evidence that these states were separate clones based on CNVs inferred from ST data (Figure S3H). This provides further support for the concept that the MP states were not solely driven by genomic aberrations at least as indicated by variations in CNVs.
Multimodal analyses of tumor states demonstrate TGFβ signaling in MP6 (Epi-Krt) tumors
ICB therapy has provided patient benefit in various cancer types, however, benefit for patients where tumors have low or absent tumor-infiltrating lymphocytes (TILs) has been limited25. Understanding the role of tumor cell heterogeneity in the exclusion/desert TIME phenotype is critical for optimizing benefit from ICB. To address this question, we probed MP6 states using multi-omics at single cell and spatial resolution.
Tumor cells in the MP6 subcluster exhibited a gene expression module (e.g FABP5, SPRR1B, SERPINB3, KRT6B and S100A8) that was enriched in keratinization process (Figure 4A, S4A). In contrast, MP7 tumor cells exhibited high expression of genes associated with antigen presentation (HLA-A, HLA-E, CD74, and HLA-DRB1), immune response genes (B2M, LCN2, and CXCL9), immune inhibitory molecules (IDO1, CD274, CD276, and BTN3A1) and multiple mucins (MUC1, MUC16, and MUC20) that are enriched in antigen processing and presentation and inflammatory response pathways (Figure 4A, S4A-B). Consistently, differential expression analysis between MP6 and MP7 spots of ST demonstrated elevated levels of FABP5, SPRR1B, SERPINB3, KRT6B and S100A8 as well as S100A9 in MP6 spots and HLA-B, HLA-C, HLA-DPB1, HLA-DQA1, HLA-DQB1, IDO1, CD74, and CXCL9 as well as CXCL10 in MP7 spots (Figure 4B and Table S6).
To determine whether elevated RNA levels resulted in concordant increases in protein, we applied high-sensitivity mass spectrometry to spatially defined tumor regions of interest (ROIs) with various degrees of TILs. 60 spatially distinct ROIs (around 100 × 100 μm2 with 60-80 cells per ROI), including 20 immune-rich and 40 immune-poor ROIs, were isolated by laser-capture microdissection (LCM) for proteomics profiling. We identified 3,541 unique proteins from at least one ROI, and up to 2,194 proteins were quantified in > 50% of all ROIs (Figure 4C). MP6-specific genes, such as FABP5, SERPINB3, KRT6B, and S100A8, were higher in immune-poor ROIs consistent with the transcriptomic data (Figure 4D, S4C and Table S7). Moreover, although only 3 of 25 MP7-specific genes were detected in spatial proteomics data, B2M and multiple MHC molecules were elevated in immune-rich ROIs (Figure 4D and S4C).
Importantly, pathway activity evaluation by PROGENy26 revealed multiple tumor-related pathways, such as androgen, p53 and PI3K, and especially the TGFβ pathway to be elevated in MP6 cells (Figure 4E). TGFβ signaling can mediate immune suppression and evasion of immunosurveillance27, and FABP5 has been reported to augment TGFβ signaling in tumor and fibroblast cells28, indicating a potential role of FABP5 and TGFβ signaling pathway on the interaction between MP6 tumor and TIME which may contribute to the MP6 related immunosuppressive TIME. PROGENy analysis of spatial transcriptome and proteome data demonstrates consistently enhanced TGFβ signaling in MP6 spots and immune-poor ROIs (Figure 4F-G). Notably, the PROGENy TGFβ was preferentially located on the periphery of MP6 tumor regions (Figure 4G), suggesting that the TGFβ pathway could alter the peritumoral microenvironment, such as CAFs and immune phenotypes.
Surrounding MMP11+ CAF contribute to the immune exclusionary microenvironment at the periphery of MP6 tumor spots
Recent studies revealed distinct CAFs phenotypes linked to diverse TIPs5,29. Desbois et al. demonstrated TGFβ-mediated CAF associated with the T cell-excluded TIP in ovarian cancer5. However, the spatial distribution and associations of CAF subtypes with tumor cells are poorly understood. Here, we identified 7 fibroblast subclusters, 3 of which (CAF-CD34, CAF-DES and CAF-MMP11) were preferentially expressed in CSCC (Figure 5A-B). Compared to the previously established functional CAFs signatures30, CAF-CD34 displayed high IL1 CAF signature expression, and CAF-MMP11 showed characteristics of myCAF and TGFβ CAF (Figure 5A). Meanwhile, PROGENy analysis confirmed activated TGFβ signaling in CAF-MMP11 (Figure 5C). Moreover, CAF-MMP11 exhibited gene expression patterns (e.g. MMP11, POSTN, CTHRC1, COL8A1, COL10A1) which have been previously associated with TGFβ-induced reactive stroma31, and pathways involved in EMT and extracellular matrix organization (Figure 5D-E). Furthermore, as previous studies showed that TGFβ-dependent LRRC15+ CAFs can directly suppress CD8+ T cell function, leading to a suppressive TIME, and limit responsiveness to checkpoint blockade29,30, our CAF-MMP11 cells expressed LRRC15 (Figure 5F) furthering supporting the immune suppressive context of CAF-MMP11.
The positive correlation between MP6 score and CAF-MMP11 cell proportion in TIME (Table S4) and the observation of activated TGFβ signaling at the periphery of MP6 tumor cell spots (Figure 4G) promoted us to investigate their in situ spatial interactions in ST. We thus scored the expression of the MMP11 gene module (CTHRC1, POSTN, INHBA, SFRP2, SULF1, COL8A1, COL11A1 and LRRC15) on fibroblastic spots at spatial resolution. Fibroblasts near MP6 tumor cells exhibited elevated CAF-MMP11 signatures (Figure 5G). By evaluating the expression of MP6 or MP7 signatures in fibroblastic neighborhoods, we observed a significantly positive correlation between CAF-MMP11 score and MP6, but not MP7 in neighborhoods (Figure 5H). These observations demonstrated CAF-MMP11 to be geographically segregated around MP6 tumors.
Spatial cellular interaction analysis demonstrated MP6/7-specific cell-to-cell communication
Cell-to-cell interactions within the TIME are critical for diverse aspects of tumor biology32,33. Intercellular communication, activated through ligands diffusing from sender cells to surrounding receiver cells, is spatially constrained34. Although scRNA-seq provides a powerful discovery tool for identification of ligand-receptor interactions, it overestimates or underestimates communication activity since the spatial information is lost9. Hence, we characterized differential intercellular interactions between MP6 and MP7 tumors with TIME directly within the tissue. For each spot in ST, we assessed the expression of 2847 ligand-receptor pairs from the CellChat database35 and derived combined expression values when both partners were expressed in a spot and its adjacent area. Next, we focused on the contrast between MP6 and MP7 spots to identify important ligand-receptor interactions associated with each tumor state. Finally, the contribution of cell types to each ligand-receptor interaction were modeled by correlation analysis between ligand/receptor expression and cell type probabilistic score of each spot (Figure 6A).
Interactions upregulated in MP6 were dominated by tumor cells, together with some contributions by macrophages, fibroblasts, and myocytes (Figure 6B-C). CEACAM1-CEACAM5, DSC2-DSG2, and COL9A6-SDC4/SDC1/ITGAV/CD44 interactions between tumor-tumor cells were increased in MP6 tumor spots (Figure 6B and Table S8). DSC2-DSG2 encode the major components of the desmosome that help resist shearing forces36. CEACAM1-CEACAM5 encode adhesion proteins important in organization of adherens junctions and tight junctions37, indicating tumor-tumor cohesive cell interactions in MP6 spots could constitute a physical barrier limiting the entry of immune cells
Interestingly, intercellular TGFβ signaling was associated with MP6 spots suggested by GDF15:TGFBR2 with MP6 expressing GDF15 interacts with TGFBR2 on adjacent CAFs, myocytes, and endothelial cells (Figure 6D). Simultaneously, MP6 associated CAFs, endothelial cells, myocytes, DC, and macrophages expressed TGFβ inducible ligands such as THBS1, FN1, VWF, and POSTN, which orchestrate stromal extracellular matrix fiber alignment38–40, with the matching receptors of ITGB6, ITGA3, ITGAV, SDC1, SDC4 receptors in MP6 tumor cells supporting tumor progression35. Consistent with the TGFβ signaling being mainly located at the periphery of MP6 spots, the ligand-receptor pairs, including GDF15-TGFBR2, and FN1-ITGB6/ITGA3/ITGAV/SDC1/SDC4 (Figure 6D), were also segregated around MP6 spots.
MP7 spots were characterized by increases in cellular communication and also diverse communication between MP7 tumor cells and NK, B cell, CD4 and CD8 T cells in addition to CAFs, endothelial cells, and myocytes (Figure S5A-B). Several of the identified ligand–receptor pairs, including CXCL9-CXCR3, CXCL9-CKR1, and CXCL11-CXCR3 have well-established roles in recruitment of tumor-infiltrating immune cells41. Moreover, MP7 tumor cells communicated with DCs, macrophages, and lymphocytes by presenting MHC-I (HLA-A/B/C) and MHC-II (HLA-DQA1/DQB1/DRB1/DPB1) (Figure S5B). Interestingly, the complement C3 pathways enriched in MP7 regions may contribute to an inhibitory TIME via dysfunctional T-cell activity (Figure S5A).
The inverse relationship of MP6 tumor cell phenotype and immune infiltration is present in multiple CSCC datasets
To determine the generalizability of the malignant MP subtype and to develop evidence for tumor MP state-specific associations with immune phenotypes in clinical settings, we first explored our published multi-omics CSCC study42. Considering challenges with single sample gene set enrichment analysis (ssGSEA) of bulk RNAseq that cannot fully distinguish contributions of malignant cells and other cells in the TIME43, we applied a new strategy, BayesPrism derived tumor cell specific gene profiles using scRNA-seq as reference44, to deconvolute the MP6 and MP7 identities in each tumor. Using this approach, the intra- and inter-tumoral heterogeneity observed in our single-cell analysis were replicated in the bulk RNA-seq data (Figure 7A, top panel). Importantly, paired proteomics and phosphoproteomics analysis demonstrated that multiple MP6 or MP7-specific transcripts were differentially expressed at protein and phosphorylation level (Figure S6A), indicating the generalizability of the malignant MP subtypes in CSCC.
A relative enrichment of MP6 (K2I score, calculated by Epi-Krt (MP6) minus Epi-Imm (MP7) score, See Methods)) was calculated and tumors were arranged according to K2I scores. Decreased immune activation (IFNG and CD8) 45and immune evasion (CD274)46 were observed in high K2I scoring tumors (Figure 7A and Table S9). We subsequently used xCell to deconvolute the TIME in tumors. Multiple immune cell types were inversely associated with the K2I score, while cell types associated with immune excluded TIME (MSC, myocytes, and keratinocytes) were elevated in tumors with high K2I scores (Figure 7A).
To determine whether the relationship between MP6 and immune exclusion was recapitulated histologically, 20 primary CSCCs from our previous published cohort with bulk RNA-seq data were randomly selected for IHC and multiplex immunofluorescence (mIF) evaluation of immune cell infiltration (Figure 7B, S6B and Table S10). CD8+ T cell infiltration was negatively correlated with K2I score (Figure7B-C). Specifically, tumors with high K2I scores, represented by sample 001A, showed sparse CD8 T cell infiltration. By contrast, samples with low K2I scores, represented by 0015 and 0018 were infiltrated by abundant CD8 T cells (Figure 7B and S6B). Importantly, the majority of findings from our analysis of scRNA-seq data were duplicated in a larger cohort of TCGA squamous cervical cancer dataset (Figure S6C). Collectively, scRNA-seq analysis, in concert with multi-omics bulk data and histological findings, supports the existence of the MP6 and MP7 states the association between these states and the degree of immune cell infiltration.
To extend these findings, we biopsied and performed bulk RNA-seq with paired pre- and post- one cycle NACT (neoadjuvant chemotherapy: cisplatin plus nabpaclitaxel) of 21 patients from the NACI clinical trial (NCT04516616) where CSCC patients were treated with one cycle of NACT followed by two cycles of neoadjuvant chemo-immunotherapy combining NACT and camrelizumab (anti-PD-1). Tumors were evaluated as pCR (pathological complete response or non-pCR groups following radical surgery. Strikingly, K2I scores were markedly decreased after one cycle of platinum-based NACT, with this decrease being most marked in the pCR group (Figure 7D and Table S11).
We subsequently asked whether MP6 and MP7 tumor states would inform response to ICB therapy in another squamous tumor lineage: squamous non-small cell lung cancer (NSCLC) 47–49. In the squamous NSCLC subtype of 2 large ICB therapy cohorts (OAK and POPLAR) with patients treated with anti-PD-L1 antibody or chemotherapy (Docetaxel), low K2I scores were associated with better survival in patients treated with immunotherapy but not in patients treated with chemotherapy in both cohorts (Figure 7E-F). Importantly, no association between K2I score and patient survival was found in either the immunotherapy or chemotherapy arms in non-squamous lung cancer (Figure S6D-E), suggesting the predictive value of K2I scores for response to immunotherapy in squamous cancers in both cervix and lung.
FABP5 is critical for MP6 state maintenance and contributes to activated TGFβ signaling and interactions with MMP11+ CAF
We next examined features that define the MP6 state. In addition to keratinization genes, this program included SERPINB3 and S100A8/A9, both of which that are associated with advanced stage and unfavorable prognosis in CSCC50,51. Suppression of these molecules has resulted in attenuated cellular proliferation and invasion50,52. Interestingly, one of the top scoring differentially expressed transcripts and proteins in MP6 was FABP5 (Figure 1C and 4A), which encodes the fatty acid binding protein that has been linked to psoriasis and various types of cancer53. Silencing FABP5 caused a substantial decrease in MP6 signature transcripts and protein levels in both cervical cancer and HNSCC cell lines (Figure 8A-B, and S7A-B). Notably, FABP5-IN-1, a selective FABP5 inhibitor54, suppressed cervical cancer patient-derived organoid (PDO) formation along with a decrease of S100A8/A9 and SERPINB3 by immunofluorescent staining (Figure 8C-D, and S7C-D). These data indicate that FABP5 may not only mark MP6 cells, but also may contribute to the acquisition and/or maintenance of the MP6 state.
In addition, genetic depletion and pharmacological inhibition of FABP5 significant reduced TGFβ, especially the TGFβ2 isoform, as well as downstream signaling targets including BMP1, ID1, ID3, and SNAIL1 in multiple cell lines (Figure 8E-F, and S7E). Immunoblotting further confirmed that FABP5 inhibition led to decreased TGFβ2 protein levels and reduced downstream SMAD2 and SMAD3 phosphorylation (Figure 8G-H, and S7F). Thus, TGFβ production and downstream signaling in MP6 tumors is at least partly mediated through FABP5. We further proceeded to culture murine fibroblast cell (L929 cells) with conditioned medium from TC-1 tumor cells with or without FABP5 knockout (KO). Parental TC-1 cells stimulated a state transition of L929 fibroblast cells to MMP11+ CAF-like cells with accompanied elevated expression of TGFβ inducible genes, including LRRC15, POSTN, BMP7, and ID1/3. Strikingly, the induction of TGFβ-inducible genes by the conditioned medium was found to be markedly decreased in FABP5 KO cells (Figure 8I). These results indicate MP6 tumors produce FABP5 dependent TGFβ production that promotes transition of CAF states to a state that drives immune exclusion.
Finally, we examined whether the modulation of MP6 states by FABP5 inhibitors can impact tumor growth and T cell infiltration in vivo. As expected, FABP5-IN-1 treatment significantly decreased tumor burden, as well as increased T cell infiltration, as demonstrated by the increased CD4 and CD8 T cells detected by flow cytometry and immunohistochemical analysis of tumors at study termination (Figure 8J-L, and S7G-H). Strikingly, Knockout of FABP5 in TC-1 cells almost completely prevented tumor formation in C57BL/6 mice (except for mouse B1) but not in immune deficient BALB/c-nu mice (Figure 8M-N and Figure S7I). Despite tumor growth in the B1 tumor, there was substantial CD4 and CD8 T cells infiltration in the FABP5 KO tumors in mouse B1when compared to tumors from parental TC-1 cells (Figure 8O). These data provide further support for the role of FABP5 in tumor immune exclusion.
The close correlation between MP7 and immune infiltration in situ suggests that MP7 tumor cells and immune elements in microenvironment may interact and potentially bidirectionally through cytokine production. Interferon gamma (IFN-γ) produced by NK and T cells in tumor microenvironment play a critical role in tumor immune surveillance55. Additionally, interferons have been used as a treatment for some types of cancer with the goal of stimulating anti-tumor immune responses56. To test our predictions, we treated squamous cervical cancer and HNSCC cells with IFN-γ as a microenvironmental perturbation. As predicted, IFN-γ suppressed the MP6 phenotype underscoring the potential interplay between tumor states and TIME. Interestingly, IFN-γ also substantially induced MP7 specific genes, as exemplified by MHCI/II, B2M, and IDO along with multiple glandular epithelium markers (MUC4, and MUC16) (Figure S7J). Moreover, IFN-γ induced MP6 and MP7 state switches were reproduced in cervical PDO which markedly decreased their formation (Figure S7K-L), raising the possibility that IFN-γ secretion by tumor infiltrating NK or T cells might induced tumor cells state transitions potentially in response to TIME cues.