Multimodal Analysis of Composition and Spatial Architecture in Human Squamous Cell Carcinoma

Abstract

To define the cellular composition and architecture of cutaneous squamous cell carcinoma (cSCC), we combined single-cell RNA sequencing with spatial transcriptomics and multiplexed ion beam imaging from a series of human cSCCs and matched normal skin. cSCC exhibited four tumor subpopulations, three recapitulating normal epidermal states, and a tumor-specific keratinocyte (TSK) population unique to cancer, which localized to a fibrovascular niche. Integration of single-cell and spatial data mapped ligand-receptor networks to specific cell types, revealing TSK cells as a hub for intercellular communication. Multiple features of potential immunosuppression were observed, including T regulatory cell (Treg) co-localization with CD8 T cells in compartmentalized tumor stroma. Finally, single-cell characterization of human tumor xenografts and in vivo CRISPR screens identified essential roles for specific tumor subpopulation-enriched gene networks in tumorigenesis. These data define cSCC tumor and stromal cell subpopulations, the spatial niches where they interact, and the communicating gene networks that they engage in cancer.

Keywords: CRISPR screen; MIBI; intra-tumoral heterogeneity; multi-omics; scRNA-seq; skin cancer; spatial transcriptomics; squamous cell carcinoma; tumor immunology; tumor microenvironment.

Graphical abstract

Figure S1

scRNA-seq and Whole-Exome Sequencing of Normal Skin and cSCC, Related to Figure 1 (A) Patient cohort characteristic table. (B) UMAP of all cells (k = 48,164) recovered labeled by patient on the top and tissue type (normal skin or tumor) on the bottom. (C) UMAP feature plots showing expression of marker genes by annotated cell types. (D) Types of mutations in select genes in patient cSCC samples. (E) Proportion of transitions and transversions mutations across patients. (F) Tumor mutational burden across patients. MB = megabase.

Figure 1

A Single-Cell Transcriptomic Atlas of Normal Skin and cSCC (A) Workflow of patient sample processing for scRNA-seq, MIBI, and ST with integration of a xenograft CRISPR screen. (B) Uniform manifold approximation and projection (UMAP) of scRNA-seq cells recovered from both normal skin and cSCC labeled by cell type. (C) UMAP of myeloid subsets, including dendritic cells (DCs), macrophages, and monocytic myeloid-derived suppressor cells (MDSCs) (marker genes in parentheses). (D) UMAP of subsets of natural killer (NK) and T cell subsets, including CD4⁺ and CD8⁺ conventional T cells and associated subsets, and regulatory T cells (Treg). Pre-Exh, pre-exhausted; T_EM, effector memory T cell; T_EMRA, recently activated effector memory T cell. (E) Bar plots of proportion of cell type by patient, tumor or normal origin, and total cell number. (F) Similar to (E), for NK cell and T cell subsets. See also Figure S1 and Table S1.

Figure 2

A Dysregulated Differentiation Hierarchy in Tumor Keratinocytes (A) Left, UMAP of normal keratinocytes (KCs) of the interfollicular epidermis and tumor KCs labeled by patient; middle, expression of basal, cycling, and differentiating genes found in both normal and tumor KCs; right, expression of tumor-specific genes. (B) UMAP classifying normal and tumor KCs into analogous basal, cycling, and differentiating clusters. Tumors contain a tumor-specific keratinocyte (TSK) subpopulation. (C) Expression of top 10 shared basal, cycling, differentiating, and TSK marker genes. TSK genes are minimally expressed in normal KCs. (D) Correlation matrix of overlapping differentially expressed genes (STAR Methods). (E) Dot plot of gene ontology (GO) terms for top 200 up- and downregulated genes in subpopulation differential expression for tumor versus normal KCs. TSK was compared to normal basal. (F) UMAP feature plots of expression for genes in hallmark epithelial-mesenchymal transition (EMT) signature in tumor KCs. (G) Violin plots of the hallmark EMT gene signature score in normal and tumor KC subpopulations. ^∗p < 2.2 × 10⁻¹⁶ by pairwise Wilcoxon rank-sum tests with Benjamini-Hochberg correction. (H) Heatmap of expression of EMT-associated transcription factors in tumor KCs. (I) Heatmap of single-cell regulon scores inferred by SCENIC (g, genes; extended, SCENIC-annotated additional genes). (J) Bar plots of percentage of cycling cells per KC subpopulation after regressing cell-cycle signature from cycling cells. (K) Bar plots of percentage of subpopulation representation of KC cycling cells after cell-cycle regression. Data shown in (J) and (K) are averages ±SEM. Normal, n = 10 patients; tumor, n = 7 patients. p values were determined by Mann-Whitney U test. For visualization, a random sampling of 100 cells per subpopulation are shown in (C), (H), and (I). See also Figure S2 and Tables S2 and S3.

Figure S2

Keratinocyte Annotation and Subpopulation Analysis, Related to Figure 2 (A) UMAP of all epithelial cell clusters (k = 29,143) with labeled cycling, eccrine cell, and pilosebaceous clusters. (B) UMAP of all epithelial cells labeled by tissue type and patient. (C) UMAP feature plots of expression of genes marking cycling, eccrine, or pilosebaceous clusters. (D) UMAPs of filtered normal and tumor keratinocytes labeled by patient before and after batch correction. CCA = canonical correlation analysis. (E) Bar plots showing proportion of each subpopulation across normal and tumor keratinocytes. (F) Violin plots of differentiation signature score (n = 387 genes, Lopez-Pajares et al., 2015) in normal and tumor keratinocyte subpopulations. p values were determined with pairwise Wilcoxon Rank Sum tests with Benjamini-Hochberg correction for multiple comparisons. (G) Clustered heatmap of all recovered transcription factor (TF) regulons from SCENIC analysis (n = 370 regulons) across tumor keratinocyte subpopulations. For visualization purposes, 100 random cells from each subpopulation are shown in the heatmap. (H) Cell cycle regression analysis with re-clustering into basal, differentiating, and TSK subpopulations with UMAP feature plots showing marker gene expression after cell cycle regression.

Figure S3

Spatial Transcriptomics Identifies TSK Localization and Patterns of Cluster Adjacency, Related to Figure 3 (A) Spatial transcriptomics (ST) spot size and resolution. (B) Violin plots of UMI counts per spot and genes per spot across tissue section replicates. (C) UMAP of all transcriptome spots labeled by patient (top) and replicate (bottom). (D) Tumor-associated spot clusters (clusters encompassing annotated tumor regions in sections), stromal or immune-associated, and non-tumor-adjacent stromal and adnexal spot clusters projected individually with labeled top differentially expressed genes. (E) Hematoxylin and eosin (H&E) staining of sections from Patients 5 and 9 with unbiased clustering of spots based on global gene expression within individual spots. Scale bar = 500 μm (F) Violin plots of TSK scores of individual spots derived from scRNA-seq data (sc-TSK score) for each cluster. Dotted boxes outline clusters with highest average sc-TSK score. (G) and (H) Overlap correlation matrix of genes differentially expressed in ST clusters across all patients (G). Highlighted similar spatial clusters were used to generate ST Cluster Signature (n = 6 genes), and violin plots of ST Cluster Signature score by cell types in scRNA-seq data (H). (I) Top, schematic of nearest neighbor analysis for spots. Bottom, heatmaps showing number of nearest neighbor identities for each cluster. ^∗indicates p < 0.001 by permutation test. (J) Visium platform ST spot size and resolution. (K) Violin plots of UMI counts per spot and genes per spot across tissue section replicates from Visium. (L) Coefficient of variation of sc-TSK score (COV_TSK) normalized to COV of KRT5 expression (COV_KRT5) across leading edge spots in patient ST data. (M) Left, spatial feature plot of sc-TSK score. Middle, spots annotated by leading edge (navy) and TSK-proximal stroma (red). Differential expression was performed with TSK-proximal stroma spots versus all other spots in tissue. Right, Violin plots of TSK-proximal stroma scores (n = 209 genes, Wilcoxon rank-sum test, adjusted p value < 0.05) by cell types in scRNA-seq data (STAR Methods).

Figure 3

Leading Edge Heterogeneity Revealed by Spatial Transcriptomics (A) Hematoxylin and eosin (H&E) staining of tissue sections and unbiased clustering of ST spots. Scale bar, 500 μm. (B) Violin plots of TSK scores of individual spots derived from scRNA-seq data (sc-TSK score) for each cluster. Dotted boxes outline clusters with highest average sc-TSK score. (C) Spatial feature plots of TSK-high cluster, sc-TSK score, and TSK marker MMP10 expression in tissue sections. (D) Violin plots of TSK-proximal signature score (intersection of differentially expressed genes in patient 2 cluster 4 and patient 10 cluster 2, n = 34 genes) for cell types in scRNA-seq data. (E) Left, H&E staining with leading edge of tumor annotated (dotted lines) and isolated leading edge spots labeled by cluster. Right, bar plots of total number and percentage of spots at leading edge per cluster. (F) Violin plots of non-TSK leading edge signature score (intersection of differentially expressed genes of non-TSK spots at leading edge from patient 2 and patient 10, n = 6 genes) by tumor subpopulations. (G) Violin plots of COL17A1 expression by tumor subpopulation in scRNA-seq data. (H) Immunohistochemical staining of COL17A1 in patient tumors. Scale bar, 50 μm. (I) Projection of non-TSK leading edge-associated clusters with dot plot of select gene ontology (GO) terms of differentially expressed genes in each cluster (n = 200 genes for patient 2 clusters, n = 20 for patient 10 cluster). (J and K) H&E, sc-TSK score, and MMP10 expression feature plots from (J) patient 4 and (K) patient 6 with data generated using the Visium ST platform. Scale bar, 500 μm. KC, keratinocyte; TAM, tumor-associated macrophage. See also Figure S3 and Table S4.

Figure 4

The Immune Landscape of cSCC (A) Violin plots of select immunosuppression-associated gene expression across cell types of the TME from scRNA-seq data. (B) Top, spatial transcriptomics (ST) cluster map from patient 2. Bottom, heatmap of select immunosuppression-associated gene expression across spatial transcriptomic spots grouped by cluster. (C) Feature plots of select immunosuppression-associated genes by ST. (D) Z-scored mean log expression heatmap of genes associated with cytotoxicity, exhaustion, and co-stimulatory function across T cell subsets in scRNA-seq data. (E) Heatmap of expression of T cell marker genes, and genes associated with cytotoxicity, exhaustion, and co-stimulatory function in ST spots grouped by cluster. (F) Left, Z-scored mean log expression heatmap of chemokine genes across cell types from scRNA-seq data. Right, Z-scored mean log expression heatmap of chemokine receptor genes across lymphocyte cell types from scRNA-seq data. Colored lines connect matching ligands for each chemokine receptor. Pre-Exh, pre-exhausted; T_EM, T effector memory; T_EMRA, T recently activated effector memory; Mig., migrating; DC, dendritic cell; CAF, cancer-associated fibroblast; KC, keratinocyte. See also Figure S4.

Figure S4

T Cell Subset Characterization and Spatial Positioning, Related to Figure 4 (A) UMAP of T cell subsets and NK cells with feature plots showing expression of subset marker genes. (B) Heatmap of top differentially expressed genes by NK cells and T cell subsets in scRNA-seq data. (C) Bar plots of proportion of cycling cells by NK or T cell subset. (D) Spatial feature plots showing expression of select genes grouped by function from Patient 2.

Figure 5

Spatial Architecture of Lymphocyte Subsets in cSCC (A) Select MIBI fields of view (FOVs) for patient samples with expression of highlighted features. (B) Heatmap of feature expression across cell types identified by MIBI. (C) Top, bar plots of proportion of non-tumor cell types across all FOVs. Bottom, bar plots of total numbers of non-tumor cells identified in each FOV. (D) Correlation heatmap of non-tumor cell types across all FOVs. (E) Scatterplot and correlation of CD8 T cell and regulatory T cell (Treg) correlation across FOVs. (F) Density plots of CD8 T cell distance to Tregs across four FOVs, demonstrating variation in CD8 to Treg co-localization. FDR, false discovery rate by permutation (STAR Methods). (G) Expression of Treg marker FoxP3, CD8, and tumor cell markers E-cadherin and pan-keratin. (H) Ranked bar plots of CD8 proportion in each FOV labeled by co-localization pattern with Tregs. (I) Co-localization patterns of Tregs and either CD4, CD8, or macrophages. (J) Non-tumor cells flagged by location relative to tumor and stromal compartments. (K) Heatmap of relative abundance of cell types in each compartment. Values represent proportion of total non-tumor cells in compartment contributed by each cell type. (L) B cells infiltrate tumor parenchyma in select FOVs. All FOVs in figure are 800 × 800 μm. All scale bars, 100 μm. See also Figure S5 and Table S5.

Figure S5

Multiplexed Ion Beam Imaging Acquisition and Analysis, Related to Figure 5 (A) Workflow for image acquisition and standard processing to achieve marker quantification for single cells. (B) Heatmap representing scRNA-seq expression for genes exhibiting similar expression patterns in MIBI and scRNA-seq. (C) Plot of cell positions and tumor or normal status in large, tiled image from Patient 8 (2mm x 2mm). Regions outlined by squares are analyzed for composition of non-tumor cells (similar to Figure 5C). (D) Distribution of mean distance to Tregs for each CD4 T cell in blue. In red, CD4 cell identities were permuted among all non-tumor, non-CD4, non-Treg cells and the mean distance to Treg was re-calculated. (E) Similar to (D), analyzing macrophage to Treg distances. (F) Comparison of CD8 versus Treg differential expression and correlation with FOXP3 across ST spots for patients 2 and 4. (G) Spatial feature plots of Treg marker FOXP3 and CD8 effector/chemokine genes demonstrating co-localization. (H) Similar to Figure 5K. Cells were labeled as stromal (tumor-excluded), leading edge (tumor-stroma border), and infiltrated in tumor parenchyma.

Figure 6

Cellular Crosstalk Landscape Associated with Leading Edge Niches (A) Schematic of combined scRNA-seq and ST ligand-receptor analysis of tumor keratinocyte (KC) subpopulations at the leading edge and TME cells. (B and C) Bar plots of significant ligand-receptor (L-R) pairs (p < 0.001, permutation test; STAR Methods) when tumors express ligands (B) or receptors (C) matched to cell types in scRNA-seq data. (D) Left, heatmap of scRNA-seq average log fold change (logFC) across normal and tumor KC subpopulations of NicheNet top predicted ligands expressed by tumor KCs that modulate TME cell types. Bottom, heatmap of scRNA-seq average logFC of ligand-matched receptors expressed by TME cell types. Middle, heatmap of significant ligand-receptor pairs between ≥1 tumor KC subpopulation and TME cell type pair in scRNA-seq. (E) UMAP scRNA-seq (N = 7 patients) and spatial feature plots (patient 4) of select ligands expressed by tumor KC subpopulations and cognate receptor expression by TME cell types. (F) Left, heatmap of scRNA-seq average logFC of NicheNet top predicted ligands expressed by TME cells that modulate the TSK signature. Bottom, heatmap of scRNA-seq average logFC of ligand-matched receptors across normal and tumor KC subpopulations (only tumor shown). Middle, heatmap of significant ligand-receptor pairs between ≥1 TME cell type and tumor KC subpopulation pair. ^∗indicates spatial transcriptomic proximity p < 0.001 (permutation test; STAR Methods) in ≥1 patient in (D) and (F). Average logFCs in heatmaps from (D) and (F) were calculated using Wilcoxon rank-sum test. (G) UMAP scRNA-seq (N = 7 patients) and spatial feature plots (patient 4) of select ligands expressed by TME cell types and cognate receptor expression by tumor KC subpopulations. (H) Scatterplot of TCGA (31 cancer types) correlation values between single genes and CAF-signature expression (y axis) (FDR <0.05; STAR Methods) and effect of CNV loss on CAF-signature expression (x axis). TSK ligands are labeled and colored by cancer type. Known CAF activator TGFB1 is highlighted in red. (I) Kaplan-Meier plot of progression-free survival after PD-1 inhibitor treatment of patient SCC (HNSC and LUSC) tumors exhibiting high and low TSK-associated gene expression (p value, chi-square test). See also Figure S6.

Figure S6

Ligand-Receptor Crosstalk in the TSK Niche, Related to Figure 6 (A) UMAP scRNA-seq (N = 7 patients) and spatial feature plots (N = 2 patients) of ligand-receptor and cell type pair expression highlighting TSK signaling to cancer-associated fibroblasts (CAFs) and endothelial cells. (B) UMAP scRNA-seq (N = 7 patients) and spatial feature plots (N = 2 patients) of ligand-receptor and cell type pair expression highlighting CAF or endothelial cell signaling to TSK. (C) Heatmap of average TME cell type signature expression across TME cell types in scRNA-seq data (D) Correlation heatmap of TME cell type signatures (columns) each correlated with the TSK signature across TCGA cancers (rows). (E) Boxplots showing RNA expression across TCGA cancers of specific genes in copy number variant (CNV) lost or stable tumors. P value was determined by Kolmogorov-Smirnov test. (F) UMAP feature plots of select ligands expressed by SCC cells and cognate receptors expressed by murine TME cell types. (G) UMAP feature plots of select ligands expressed by murine TME cell types and cognate receptors expressed by SCC cells.

Figure 7

In Vivo CRISPR Analysis of Tumor Keratinocyte Vulnerabilities (A) Workflow of xenograft mouse models of human SCC cancer cells and CRISPR screen of tumor keratinocyte subpopulation-specific genes. (B) UMAP feature plots of signature scores derived from patient scRNA-seq data for basal, cycling, differentiating, and TSK scores in xenograft tumor cells. High-scoring subpopulations of cells are highlighted by dotted circle. (C) Differentially expressed gene (DEG) overlap correlation matrix across patient tumor subpopulations and unbiased clusters of xenografted SCC cell lines (STAR Methods). (D) UMAP of murine TME cells isolated from xenograft tumors. Top, labeled by cell type. Bottom, labeled by xenografted SCC cell line. (E) DEG overlap correlation matrix across patient and murine TME cell types demonstrating conserved cell-type expression programs across patient and xenograft TME cells. (F) Ratio of absolute coefficient of variation for subpopulation signatures between xenografted CAL27 cells and in vitro cultured CAL27 cells. (G) Left, heatmaps of mean log fold-change of sgRNAs targeting depleted and enriched genes in SCC cell lines (STAR Methods). All genes were significant by false discovery rate (FDR from STARS algorithm) <0.05 in ≥1 cell line. Right, heatmap of gene expression in patient scRNA-seq. (H) Co-essentiality network plot of depleted genes with STARS FDR <0.10 in ≥1 xenograft tumor line. Larger nodes indicate inclusion in screen and are colored by subpopulation designation while smaller node and font labels indicate connecting genes not included in screen. Gray larger nodes indicate broadly expressed genes. Edges represent a co-essential relationship at a previously determined FDR <0.10 (Wainberg et al., 2019). See also Figure S7.

Figure S7

Analysis of In Vivo and In Vitro CRISPR Screen, Related to Figure 7 (A) Heatmap displaying average log fold-change of each tumor subpopulation-identifying gene relative to other subpopulations in corresponding patient sample or SCC cell line. (B) Heatmap displaying average log fold-change of each TME cell type-identifying gene relative to other cell types in either patient tumors or murine xenografts. (C) UMAP feature plots displaying patient-derived tumor subpopulation signatures for SCC cells harvested from xenograft tumors or in vitro culture. (D) Heatmap representing expression of genes targeted in CRISPR screen. Genes are grouped as tumor subpopulation-specific or broadly expressed (STAR Methods). (E) Volcano plots representing STARS false discovery rates (FDRs) and log2 fold-change (log₂FC) values for each in vitro xenograft screen. Red dotted lines represent FDR 0.10. (F) Log2 fold-change (log₂FC) values of select genes from CRISPR screens in vivo versus in vitro. Error bars represent standard deviation across tumor or cell culture biological replicates. Significance determined by permutation (STAR Methods). (^∗∗: Permutation FDR < 0.05, ^∗: Permutation FDR < 0.1). (G) Comparison of mean CNV gain and loss proportions across 31 TCGA tumor types for all genes. TSK-expressed CRISPR hits and control recurrently altered genes are highlighted. (H) Kaplan-Meier plots for TCGA tumor types with patients stratified by expression of TSK signature. P value derived from chi-square test.