Single-cell analysis of human basal cell carcinoma reveals novel regulators of tumor growth and the tumor microenvironment

Abstract

How basal cell carcinoma (BCC) interacts with its tumor microenvironment to promote growth is unclear. We use singe-cell RNA sequencing to define the human BCC ecosystem and discriminate between normal and malignant epithelial cells. We identify spatial biomarkers of tumors and their surrounding stroma that reinforce the heterogeneity of each tissue type. Combining pseudotime, RNA velocity-PAGA, cellular entropy, and regulon analysis in stromal cells reveals a cancer-specific rewiring of fibroblasts, where STAT1, TGF-β, and inflammatory signals induce a noncanonical WNT5A program that maintains the stromal inflammatory state. Cell-cell communication modeling suggests that tumors respond to the sudden burst of fibroblast-specific inflammatory signaling pathways by producing heat shock proteins, whose expression we validated in situ. Last, dose-dependent treatment with an HSP70 inhibitor suppresses in vitro vismodegib-resistant BCC cell growth, Hedgehog signaling, and in vivo tumor growth in a BCC mouse model, validating HSP70's essential role in tumor growth and reinforcing the critical nature of tumor microenvironment cross-talk in BCC progression.

Fig. 1.. Cellular characterization of human BCC subtypes using scRNA-seq.

(A) Schematic representation of in toto epithelial and stromal tissue isolation and processing from human PTS and BCC tissues for 3′-droplet–enabled single-cell RNA sequencing (scRNA-seq). (B) Two-dimensional clustering of single cells isolated from individual human BCC subtypes. IDs represent subtype and donor. BCC subtypes are color-coded on the basis of subtype and donor and include the following: superficial, nodular, and infiltrative (BCC-I); superficial and nodular (BCC-II); unknown/“hybrid” (BCC-III); and infiltrative with perineural invasion (BCC-IV). Ten distinct meta-clusters are identified at distinct proportions across BCC subtypes and annotated with their putative identities. The putative identity of each cell meta-cluster is defined on the bottom and color-coded accordingly. (C) Copy number variant analysis of putative malignant epithelial cells with InferCNV. Blue indicates low modified expression, corresponding to genomic loss; red indicates high modified gene expression, corresponding to genomic gain. Internal reference cells refer to nonepithelial, nonimmune control cells. Observations refer to putative malignant epithelial cells. Genomic regions (chromosomes) are labeled and color-coded. (D) Clustering of corrected and integrated PTS and BCC datasets is grouped by condition and donor using scMC. Conditions and donor are labeled and color-coded. (E) Two-dimensional clustering reveals cellular heterogeneity of integrated human PTS and BCC datasets. Ten distinct metaclusters are identified at various proportions across BCC subtypes and annotated with their putative cell type identities. The putative identity of each cell meta-cluster is defined on the right and color-coded accordingly per cell type. (F) Proportion of cell types grouped by condition. (G and H) Feature plots showing bona fide genes (G) and BCC-specific epithelial markers (H). Gray, low normalized gene expression based on normalized counts; black, high normalized gene expression based on normalized counts.

Fig. 2.. Comparison of epithelial cells reveals regulators of malignancy in human BCC.

(A and B) Clustering of 30,058 corrected KRT14⁺ epithelial cells from human PTS and BCC subtypes grouped by condition and donor. Fifteen putative KRT14⁺ epithelial cell identities, including 1 proliferating epithelial and 3 interfollicular epithelial cells, and 11 basal/basaloid epithelial cells were identified and defined. PTS/BCC agglomerative clustering shows relationships between KRT14⁺ epithelial cells. Cells are color-coded accordingly. (C) Dot plot of top two marker genes identified by differential gene expression among epithelial cells. Gray, low average gene expression; purple, high average gene expression. Size of circle represents the percentage of cells expressing gene markers of interest. (D to F) Protein immunostaining of select BCC–epithelial cell markers shows cluster specificity and distinct spatial localization in human primary clinical tumors. Inset shows magnified area of BCC nest. White arrows point at epithelial cells expressing KRT14. Yellow arrows point at epithelial cells coexpressing protein of interest and KRT14. Tissues were counterstained with DAPI. Scale bars, 100 μm. (G to I) Heatmap of condition and donor-specific active gene regulatory networks demonstrates differentially active FOX, HOX, and SOX regulons in BCC epithelial cells compared to PTS (Z score > 0). Yellow, low regulon activity; blue, high regulon activity; white, absent regulon activity.

Fig. 3.. Analysis of stromal cells highlights FIB and FIB-like cell heterogeneity in human BCC.

(A) Clustering of 7080 corrected FIB/FIB-like (FIB/FIB-like) cells from human PTS and BCC subtypes grouped by condition and subtype. Four putative PDGFRA⁺ FIB and two putative RGS5⁺ FIB-like cell identities were identified and defined. (B and C) Quantification and agglomerative clustering of FIB/FIB-like cells. Bar graph represents cell average per donor per cluster ± SEM. Unpaired Student’s two-tailed t test. n.s., not significant. (D) Dot plots of canonical/marker genes in FIB/FIB-like cells. Blue, low-average gene expression; red, high-average gene expression. (E to H) Feature plots and in situ RNA/protein staining show FIB marker specificity/distinct spatial localization in human primary clinical tumors. Inset shows magnified area in BCC nests. White arrows point at FIBs expressing gene/protein of interest. Tissues were counterstained with KRT14 (RNA/protein) and DAPI. Scale bars, 100 μm. (I) Pseudo-bulk dot plots of ECM remodeling genes. Red, low-average gene expression; blue, high-average gene expression. Size of circle, percentage of expressing cells. (J) Heatmap of differentially expressed genes in FIB IV cells. Yellow, down-regulated genes; red, up-regulated genes. (K and L) Gene expression (K) and cellular density (L) plots of TMEM119, WNT5A, or TMEM119;WNT5A cells. Purple, low cellular density; yellow, high cellular density. (M) RNA in situ hybridization of WNT5A in human primary clinical tumors. Inset shows magnified area of BCC cells. White arrows point at WNT5A⁺ FIBs. Tissues were counterstained with KRT14 and DAPI. Scale bars, 100 μm. (N and O) Heatmap showing active regulons in FIB/FIB-like cells (Z score > 0). Yellow, low regulon activity; blue, high regulon activity; white, absent regulon activity. Regulon activity was used for dimensionality reduction in a two-dimensional embedding. White arrows mark BCC-specific FIBs (IV). Purple, low regulon activity; yellow, high regulon activity; density plots, AUC distribution per regulon.

Fig. 4.. RNA dynamics analyses reveal differential stromal developmental trajectories in human BCC.

(A and B) Three-dimensional Waddington energy (i.e., entropy) landscape of human PTS and BCC. Blue, low entropy; blue, high entropy. (C) Quantification of cellular energy. Color of circles corresponds to distinct FIB/FIB-like clusters. Dashed lines connect FIB/FIB-like clusters and are color-coded on the basis of type of condition (i.e., PTS versus BCC). (D to G) Modeling of initial states in FIB/FIB-like cells suggests distinct developmental trajectories in PTS and BCC stroma. Arrows representing direction of cells’ flow of PAGA-velocity graph were projected as vector field on a two-dimensional embedding. In PTS, bidirectional path of FIBs is represented by trajectory 1 ($\overrightarrow{V}$_{PTS, Trajectory 1}) and 2 ($\overrightarrow{V}$_{PTS, Trajectory 2}). In BCC, unidirectional path of FIBs is represented by trajectory 1 ($\overrightarrow{V}$_{BCC, Trajectory 1}). (H and I) Rolling-wave plots identify pseudo-dependent TFs overexpressed in PTS (H) and BCC (I) along developmental trajectory 1 and grouped depending on their dynamics (k = 3 in PTS; k = 4 in BCC). Pseudotime levels are based on normalized counts. Blue, down-regulated TFs; red, up-regulated TFs. (J) Comparison of significant pseudo-dependent TFs overexpressed in PTS and BCC developmental trajectories in specific groups. TF dynamics are color-coded on the basis of condition. (K) Significant pathway ontologies associated with PTS and BCC FIB developmental trajectory 1 (P_adj < 0.05). Specific pathway ontologies in BCC are color-coded on the basis of significance. Adjusted P value scale shown on the right is based on a rainbow scale. Purple, low significance; red, high significance. (L) TGF-β, inflammation, and noncanonical WNT pathway scores based on normalized counts overlaid on two-dimensional embedding with RNA velocity streams reveal specific pathway programs associated with PTS and BCC stromal developmental trajectory 1. Yellow, low score; black, high score.

Fig. 5.. Epithelial-FIB communication modules in human BCC.

(A) Heatmap of active signaling pathways in epithelial-FIB cross-talk from human PTS and BCC samples. Blue, active signaling pathway; red, inactive signaling pathway. (B) Ranking of active signaling pathways in PTS and BCC based on their overall information flow within the inferred cellular networks. Signaling pathways are colored according to condition where they are enriched, whereas those in black are enriched equally across conditions. (C) Joint clustering of active signaling pathways from PTS and BCC into a shared two-dimensional manifold according to their functional signaling relationship similarity (k = 4). Circles represent PTS signaling pathways; squares represent BCC signaling pathways. Each shape represents the communication network of one signaling pathway. A magnified view of each cluster with labeled active signaling pathways is shown on the right. (D) Circle plots show ncWNT signaling in sending and receiving cells. Nodes are colored similarly as senders. Size of cell clusters is representative of the number of active cells in signaling network. ncWNT is active in BCC but not in PTS. Cell types participating in signaling pathway network are labeled. Bar graphs show relative contribution of specific ligand-receptor pairs for ncWNT signaling in BCC. WNT5A ligand is the only active ligand in the ncWNT signaling network. (E) Dot plots show cross-talk probability between FIBs (senders) and epithelial cells (receivers) via ncWNT signaling. Blue, low cross-talk probability; red, high cross-talk probability. Size of circle represents the percentage of cells with high cross-talk probability. Ligands are colored aqua blue; receptors are colored magenta. (F to I) Circle plots and network centrality analysis for CXCL (F), IL6 (G), IFN-I (H), and TNF (I) signaling. Only cell clusters participating in signaling network are labeled. Inactive pathway indicates that the pathway is not active. Heatmaps represent network centrality. White, low importance; green, high importance.

Fig. 6.. Heat shock proteins are prominent regulators of BCC.

(A) Heatmap of pseudo-bulk HSP70-coding gene expression in human PTS versus BCC epithelial cells. (B) In situ expression of HSP70 protein shows distinct spatial localization in human primary clinical tumors. Inset shows magnified area in BCC nest. White arrows point at HSP70⁺ epithelial cells in BCC nest. Scale bars, 100 μm. (C) Western blot and quantification of RFI (Relative Fluorescence Intensity) against Hsp70 in ASZ001 murine cells treated with HSP inhibitor Ver155008. β-Actin served as loading control. Mann-Whitney test (**P = 0.007). (D to F) HSP inhibitor Ver155008 negatively affects growth of ASZ001 murine cells (D) (two-way ANOVA test; ****P < 0.0001) and down-regulates Gli1 mRNA (E) and protein expression (F) in vitro in a concentration-dependent manner (unpaired Student’s two-tailed t test with Welch’s correction; *P < 0.05 and ***P < 0.001). Tubb served as loading control. Experiments were repeated at least three times, and data are represented as the means ± SEM. (G and H) HSP inhibitor Ver155008 significantly induces apoptosis via Casp3 (H) and negatively affects proliferation (G) via Mki67 in ASZ001 murine cells in vitro in a concentration-dependent manner. Bar graphs represent the mean of nine replicate wells ± SEM. Unpaired Student’s two-tailed t test (*P < 0.05 and **P < 0.01). (I) Schematic representation of microtumor development and HSP inhibitor treatment in Gli1-Cre^ERT2;Ptch1^fl/fl mice. Ver155008- and vehicle-treated dorsal skin tissues were collected and assessed for microtumors. (J) H&E of Ver155008- and vehicle-treated Gli1-Cre^ERT2;Ptch1^fl/fl mouse dorsal skin tissues. Scale bars, 100 μm. Quantification of microtumor surface area in vehicle- and Ver155008-treated Gli1-Cre^ERT2;Ptch1^fl/fl mouse dorsal skin tissues. Bars represent average individual microtumor area ± SEM. Surface area decreased in a concentration-dependent manner compared to vehicle-treated control mice. Unpaired Student’s two-tailed t test (*P < 0.05).