c-FOS drives reversible basal to squamous cell carcinoma transition.

While squamous transdifferentiation within subpopulations of adenocarcinomas represents an important drug resistance problem, its underlying mechanism remains poorly understood. Here, using surface markers of resistant basal cell carcinomas (BCCs) and patient single-cell and bulk transcriptomic data, we uncover the dynamic roadmap of basal to squamous cell carcinoma transition (BST). Experimentally induced BST identifies activator protein 1 (AP-1) family members in regulating tumor plasticity, and we show that c-FOS plays a central role in BST by regulating the accessibility of distinct AP-1 regulatory elements. Remarkably, despite prominent changes in cell morphology and BST marker expression, we show using inducible model systems that c-FOS-mediated BST demonstrates reversibility. Blocking EGFR pathway activation after c-FOS induction partially reverts BST in vitro and prevents BST features in both mouse models and human tumors. Thus, by identifying the molecular basis of BST, our work reveals a therapeutic opportunity targeting plasticity as a mechanism of tumor resistance.

INTRODUCTION

Tumor evolution in growth and differentiation dependencies remains the main challenge for successful targeted tumor therapies. Efforts at early detection, precise identification of tumor type, and clear understanding of initial cancer sensitivities fail as tumors toggle between cellular identities, resulting in relapse and recurrence (Boumahdi and de Sauvage, 2020; Shen et al., 2020). Basal-to-squamous carcinoma transdifferentiation is a common epithelial tumor lineage reprogramming, occurring spontaneously or upon targeted therapy in adenocarcinomas of the lung (Hou et al., 2016; Nakagawa et al., 2003; Shundo et al., 2011), stomach (Mori et al., 1987), esophagus (Kay, 1968; Takeuchi et al., 2019), pancreas (Madura et al., 1999; Rahemtullah et al., 2003), prostate (Lee, 2019), and skin (Tan et al., 2017). Poor patient outcomes associated with these alterations in tumor cell identity and sensitivities present clear unmet clinical needs, requiring further mechanistic understanding (Boumahdi and de Sauvage, 2020; Shen et al., 2020).

In the skin, basal cell carcinoma (BCC) growth depends on the hyperactivation of the Hedgehog (HH) pathway, motivating the successful clinical introduction of Smoothened inhibitors (SMOi) (Migden et al., 2015; Rubin and de Sauvage, 2006; Sekulic et al., 2012). However, advanced BCCs harbor a high mutational burden and heterogeneity and typically resist or relapse after SMOi treatments (Bonilla et al., 2016; Chang and Oro, 2012; Sekulic et al., 2012), providing an ideal model system for understanding tumor evolution. Several BCC resistance mechanisms have been identified, including canonical HH pathway activation (Atwood et al., 2015; Sharpe et al., 2015; Yauch et al., 2009); non-canonical mechanisms to bolster HH pathway output through enhanced PI3K, aPKC signaling, or G-actin-mediated activation of the myocardin-related transcription factor (MRTF)-serum response factor (SRF) complex (Atwood et al., 2013; Buonamici et al., 2010; Whitson et al., 2018); and a growth-arrested Wingless and Int-1 (Wnt)-dependent persister state (Biehs et al., 2018; Sánchez-Danés et al., 2018). BCC transdifferentiation into squamous cell carcinoma (SCC), which we refer to as basal to squamous cell carcinoma transition (BST), represents a distinct and important form of resistance for BCCs, switching away from the universal dependency on the HH pathway (Ransohoff et al., 2015; Saintes et al., 2015). BST frequently occurs in SMOi-treated advanced BCCs (Ransohoff et al., 2015; Saintes et al., 2015) but also spontaneously in aggressive sporadic basosquamous carcinomas (BSCs) (Garcia et al., 2009; Tan et al., 2017; Wermker et al., 2015). Indeed, we recently observed this lineage and molecular plasticity with HH to Ras/MAPK pathway switching in both resistant BCCs (rBCCs) and BSCs (Chiang et al., 2019; Kuonen et al., 2019; Zhao et al., 2015), although the molecular basis of the switch remains mysterious. A more detailed understanding of BST will help uncover new strategies to prevent tumor lineage transition and tumor resistance and to improve the patient’s prognosis.

Here, we use multi-dimensional genomics of patient and mouse tumors to identify the transcription factor (TF) network, chromatin accessibility (“epigenetic”) remodeling, and transcriptional output defining BST. We show that the surface markers TACSTD2, LYPD3, and LY6D, previously associated with MRTF-mediated resistance (Yao et al., 2020), extensively cover the morphological and transcriptional cell states of BST. Mechanistically, we identify c-FOS as a central player of BST in mouse and human tumors and show that c-FOS-mediated BST may be partially reverted by EGFR or MAPK inhibitors, offering a therapeutic opportunity and additional insights to prevent tumor plasticity in patients.

RESULTS

BCC-RMs identify BCC tumor clones with SCC features

Because of the occurrence of squamous transdifferentiation in many adenocarcinomas (Hou et al., 2016; Kay, 1968; Lee, 2019; Madura et al., 1999; Mori et al., 1987; Nakagawa et al., 2003; Rahemtullah et al., 2003; Shundo et al., 2011; Takeuchi et al., 2019; Tan et al., 2017) and its particular relevance to BCC resistance and aggressiveness (Kuonen et al., 2019; Zhao et al., 2015), we sought to identify the transcriptional reprogramming underlying BST. We first built bulk transcriptomic signatures for BCC and SCC based on the differential gene expression analysis between human BCCs (Bonilla et al., 2016) and SCCs (Chitsazzadeh et al., 2016) (Figures S1A and S1B; Table S1). We then superimposed these signatures onto single-cell datasets of pooled naive BCCs (nBCCs) and rBCCs processed for tumor epithelial cells as previously shown (Yao et al., 2020; Yost et al., 2019). BCC tumor cells harbored significant heterogeneity in terms of BCC and SCC signature enrichment, irrespectively of being naive or resistant (Figures 1A and S1C). However, compared to nBCCs, rBCCs displayed a wider spectrum of BCC and SCC signature enrichment (Figures 1B and S1C), reflecting the various strategies BCC tumor cells may adopt to escape HH inhibition. Typically, resistant GLI1pos BCC6 maintained a BCC-like transcriptomic profile, while resistant GLI1neg BCC8 showed an SCC-like transcriptomic profile (Figure S1D).

Figure 1.

BCC-RMs identify BCC tumor clones with SCC features

(A) scRNA-seq projection plot of tumor cells obtained from both nBCC and rBCC analyzed for SCC (red) versus BCC (blue) enrichment scores. Human gene signatures and differential gene expression between human BCCs and SCCs are shown in Table S1 and Figures S1A and S1B. Identification of naive versus SMOi-resistant tumor cells on the scRNA-seq projection plot are shown in Figure S1C.

(B) Violin plots showing BCC and SCC enrichment scores in rBCC versus nBCC. p values calculated using two-sided Wilcoxon rank sum test.

(C) scRNA-seq projection plot shown in (A) analyzed for the BCC-RM TACSTD2, LYPD3, and LY6D.

(D) Principal-component analysis (PCA) of bulk RNA-seq data obtained from human normal skin (gray, n = 4), BCCs (blue, n = 6), BSCs (green, n = 3), well-differentiated SCCs (orange, n = 3,) and poorly differentiated SCCs (red, n = 5). Mean values (centroids) are represented by bigger dots.

(E) DAPI, TACSTD2, LYPD3, LY6D, GLI1, CD44, and MUC1 protein expression indicated by immunofluorescence staining and quantified in Figure S1E on early (Gorlin) and advanced BCCs, BSCs, and well-differentiated and poorly differentiated SCCs. H&E staining shows representative pictures of tumor phenotypes. Scale bars:100 μm.

(F) Representative view of PCA euclidean distances and relative markers expression from early BCCs to poorly differentiated SCCs.

We previously identified TACSTD2, LYPD3, and LY6D (BCC-resistant markers [BCC-RMs]) as markers of MRTF-mediated resistance to SMOi in advanced BCCs (Yao et al., 2020). The epigenetic proximity of BCC-RMpos resistant cells with basal and suprabasal transit-amplifying cells rather than hair follicle germ or bulge cells (Yao et al., 2020), as well as the conserved expression of the BCC-RM in the suprabasal layers of the inter-follicular epidermis (Kriegbaum et al., 2016; Shvartsur and Bonavida, 2015; Stepan et al., 2011), suggests that BCC-RMs may reflect keratinocyte commitment. Accordingly, we observed a substantial overlap between BCC-RM expression and SCC signature enrichment in BCCs (Figures 1C and S1D). To further explore the correlation between BCC-RM expression and SCC differentiation, we looked at bulk transcriptomic data obtained from histologically confirmed human early (Gorlin) BCCs, advanced BCCs, BSCs, and well-differentiated or poorly differentiated SCCs (Bonilla et al., 2016; Chitsazzadeh et al., 2016). At the transcriptional level, BCC tumors appeared closer to well-differentiated compared to poorly differentiated SCC tumors by principal-component analysis (PCA) (Figure 1D). With mixed basal and squamous differentiation features, BSC tumors appeared between BCCs and well-differentiated SCC tumors (Figure 1D). Importantly, we observed a progressive reduction in GLI1 (Chiang et al., 2019; Kuonen et al., 2019) and a progressive increase in the SCC markers CD44 and MUC1 (Beer et al., 2000) and various epidermal differentiation markers along the squamous differentiation path (Figures 1E, S1E, and S1F). As previously suggested by single-cell RNA sequencing (scRNA-seq), BCC-RM expression was specifically found in advanced BCCs, BSCs, and well-differentiated SCCs, with a gradual increase along squamous differentiation (Figures 1E, S1E, and S1F). Altogether, we conclude that differential BCC-RM expression defines the transcriptional and morphological transition from advanced BCC to well-differentiated SCC (Figure 1F).

Modulation of epithelial TGF-β and/or Ras/MAPK signaling in BCC-RMpos clones drive BST

To elucidate the molecular drivers of BST, we first aligned BCC-RM-associated clusters and then examined the associated signaling pathways and TFs to further assay in our functional assays. To get an insight into the dynamics of naturally occurring BST in BCC-RMpos tumor clones, we focused on BCC3, a naive tumor exhibiting a wide spectrum of BCC and SCC transcriptomic differentiation, the latter overlapping with BCC-RM expression (Figure 2A). A similar pattern of BCC and SCC transcriptomic differentiation, where the latter overlaps with BCC-RM expression, was found in the nBCC8 (shown in Figure S2A). Interestingly, RNA velocity analysis indicated that the BCC-RMpos clusters (clusters 3, 7, and 4) had the largest nuclear/cytoplasmic RNA dynamics from BCC toward SCC transcriptomic enrichment (Figure 2B). Among BCC-RMpos clusters, cluster 3 (1058 top-enriched genes; Table S2) displayed similarity with human BCCs, while clusters 7 and 4 (149 and 542 top-enriched genes, respectively; Table S2) showed stronger similarity with human BSCs and well-differentiated SCCs (Figure 2C).

Figure 2.

Modulation of epithelial TGF-β and/or Ras/MAPK signaling in BCC-RMpos clones drive BST

(A) scRNA-seq projection plots of tumor cells obtained from a nBCC analyzed for SCC (red) versus BCC (blue) enrichment score (left panel) and corresponding BCC-RM TACSTD2, LYPD3, and LY6D expression (right panels). scRNA-seq projection plots of tumor cells obtained from an additional nBCC sample are shown in Figure S2A.

(B) RNA velocity analysis on tumor cells obtained from the nBCC shown in (A). Note the bidirectional nature of the BCC cluster C3 and SCC clusters C4 and C7. Main directionalities are emphasized using magnified red arrows.

(C) Relative enrichment of the top-enriched gene lists obtained from clusters 3, 4, and 7 in human BCC (blue, n = 6), BSC (green, n = 3) and well-differentiated SCC (orange, n = 3) bulk RNA-seq. Left panel shows, at higher magnification, the scRNA-seq projection plot of clusters 3, 7, and 4 obtained from the nBCC shown in (A) and (B). Boxes represent the mean and the distribution of values from minimum to maximum. Markers associated with the indicated clusters are listed in Table S2.

(D) GSEA for cancer-related canonical pathways in clusters 4 and 7 compared to cluster 3 (full bars) and human well-differentiated SCCs compared to BCCs (empty bars). Blue and red bars indicate consistently reduced or enriched pathways, respectively.

(E) TACSTD2, LYPD3, LY6D, GLI1, CD44, and MUC1 protein expression indicated by immunofluorescence staining on ASZ suggest their advanced BCC states. Scale bars: 25 μm.

(F) Relative enrichment for cluster 3 (green) and clusters 4 and 7 (purple) signatures in RNA-seq data obtained from ASZ upon Ras/MAPK activation and/or TGF-β inhibition (through TGF-βRIshRNA) (n = 2 biological replicates). RNA-seq profiles are shown in Figure S2F.

(G) H&E, GLI1, ARL13B (for primary cilia), CD44, and MUC1 staining of ASZ xenografted tumors upon Ras/MAPK activation and/or TGF-β inhibition (through TGF-βRIshRNA). White arrowheads indicate primary cilia. Scale bars: 100 μm.

(H) Quantification of GLI1, CD44, and MUC1 protein expression in (G) by pixel intensity measurements (n ≥ 40 cells, measured in n ≥ 3 fields, across n ≥ 3 tumors for each condition). Quantification of primary cilia (ARL13B) in (G) by the percentage of ciliated tumor cells (measured in n ≥ 6 fields, across n ≥ 3 tumors for each condition). Horizontal bars indicate the mean ± SD.

(I) Relative enrichment for BCC (blue) and SCC (red) signatures in RNA-seq data obtained from ASZ upon Ras/MAPK activation and/or TGF-β inhibition (through TGF-βRIshRNA) (n = 2 biological replicates). Mouse gene signatures and differential gene expression between mouse BCC and SCC are shown in Table S3 and Figures S2I and S2J. GSEA plots are shown in Figures S2K and S2L.

Horizontal bars in (F) and (I) represent the mean. p values in (C), (F), (H) and (I) calculated using one-way ANOVA with Tukey’s post-test. *p < 0.05, **p < 0.01, ***p < 0.001.

To elucidate key BST signaling pathways, we used gene set enrichment analysis (GSEA), comparing clusters 4 and 7 to cluster 3 in nBCC3 (Figure 2C, left panel), clusters 2 and 5 to cluster 0 in nBCC8 (Figure S2B, left panel; Table S2), as well as human well-differentiated SCCs to human BCCs for cancer-related canonical pathways. In accordance with the current understanding of keratinocyte-derived tumors (Cammareri et al., 2016; Campbell et al., 2018; Hosseini et al., 2017; South et al., 2014), we mostly found a reduced representation of the HH, Wnt, tumor growth factor (TGF)-β, and Notch pathways and an increased representation of Ras/MAPK and tumor necrosis factor (TNF)-α pathways (Figures 2D and S2B). Inspired by these data, we exposed an “advanced” GLI1pos, BCC-RMint (TACSTD2pos, LYPD3pos, LY6Dneg), CD44neg, MUC1neg BCC cell line (Figure 2E)—which, according to human data, would be prone to undergo BST—to inhibitors of the HH (vismodegib, 1 μM), Wnt (XAV-939, 10 μM), TGF-β (SB431542, 10 μM), and Notch (RO4929097, 10 μM) pathways or activators of the Ras-MAPK pathway 12-O-tetradecanoylphorbol-13-acetate (TPA, 50ng/ml) for 6 days (Figure S2C). Remarkably, TGF-β inhibitor SB431542 and TPA had a similar effect in reducing nuclear GLI1 while inducing CD44 and MUC1 (Figure S2C), supporting their involvement in BST. Moreover, TGF-β inhibition and Ras activation are reported to drive SCC outgrowth in mouse models (Cammareri et al., 2016; Guasch et al., 2007). Based on these encouraging results, we engineered a doxycycline-tamoxifen-inducible RasV12 (Kuonen et al., 2019) and a TGF-βR1-targeting short hairpin RNA (shRNA) (Figures S2D and S2E). Remarkably, RNA-seq profiles obtained upon Ras/MAPK activation and TGF-β inhibition showed reduced enrichment of the cluster 3 signature and increased enrichment of the clusters 4 and 7 signature (Figures S2F and 2F), supporting the human relevance of our experimental model. While control BCC cells typically form well-defined infiltrative-like BCC when xenografted (Kuonen et al., 2018), Ras/MAPK activation and/or TGF-β inhibition both induced BST in vivo, as shown by paler and larger nuclei, eosinophilic cytoplasmic changes, keratohyalin granules, and keratin pearls (Figure 2G). As predicted from the human data, in both our in vivo xenograft and our in vitro models, Ras/MAPK activation and/or TGF-β inhibition reduced nuclear GLI1 and ARL13B (for primary cilia) while inducing CD44 and MUC1 (Figures 2G, 2H, S2G, and S2H). To further confirm BST at the transcriptomic level, we built bulk transcriptomic signatures for mouse BCC and SCC, based on the differential gene expression analysis between mouse SCC (7,12-dimethylbenzanthracene (DMBA)-treated non-obese diabetic severe combined immunodeficiency (NOD-SCID) model; Nassar et al., 2015) and BCC (Ptch+/− ;K14-Cre-ER;p53 model; Whitson et al., 2018) tumor models (Figures S2I–S2K; Table S3). Based on these signatures, Ras/MAPK activation and/or TGF-β inhibition induced a switch from BCC to SCC transcriptional profile (Figures 2I and S2L), although more efficiently together than individually. Altogether, BCC-RMs point to key transcriptional reprogramming of BST, driven by Ras/MAPK activation and/or TGF-β inhibition.

Transcriptional and chromatin accessibility changes identify activator protein (AP-1) as a critical regulator of BST

To identify the TF network controlling the dynamic chromatin landscape (Li et al., 2019) associated with BST, we compared RNA-seq and assay for transposase-accessible chromatin using sequencing (ATAC-seq) changes upon Ras/MAPK and/or TGF-β inhibition (Figures S2F and S3A). Significant chromatin remodeling occurred in regions enriched with hair follicle and epidermal TF recognition motifs (Figures 3A and S3B) (Adam et al., 2020; Appleford and Woollard, 2009; Folgueras et al., 2013; Fortunel et al., 2019; Gerdes et al., 2006; Hwang et al., 2008; Kadaja et al., 2014; Leishman et al., 2013; Li et al., 2019; Sastre-Perona et al., 2019; Soares and Zhou, 2018; Yang et al., 2015; Young et al., 2017). Looking specifically at the genomic proximity of differentially expressed genes upon BST induction (downregulated [DN] and upregulated [UP] genes; see Table S4), hair-follicle-related TF motifs such as Sox9 and Lhx showed reduced accessibility, mainly in the proximity of downregulated genes. In contrast, epidermal-related TF motifs like E-twenty-six-specific sequence (ETS) had increased accessibility, mainly in the proximity of upregulated genes (Figure 3B). Strikingly, AP-1 appears as the most significantly enriched motif in both closing and opening distinct chromatin regions during BST (Figures 3A and S3B). AP-1-enriched peaks that closed during BST were associated with the hair-follicle-associated downregulated genes, while a distinct set of AP-1-enriched peaks that opened were associated with the epidermal upregulated genes (Figures 3B and 3C). This observation suggests that BST involves the differential regulation of AP-1-related chromatin-accessible sites associated with hair follicle and epidermal gene subsets.

Figure 3.

Transcriptional and chromatin accessibility changes identify AP-1 as a critical regulator of BST

(A) Enrichment of TF motifs in the regions with decreased (left) or increased (right) chromatin accessibility in ASZ upon Ras/MAPK activation and TGF-β inhibition. The y axis is −log10(p value) of a motif enrichment score. The x axis is the ranking number of sorted motifs. Motifs belonging to one TF family are labeled with the same color. ATAC-seq profiles are shown in Figure S3A. Similar analysis upon Ras/MAPK activation and TGF-β inhibition individually is shown in Figure S3B.

(B) Chromatin accessibility changes with enrichment for the specified motifs and c-FOS ChIP-seq peaks found in the proximity of genes downregulated (722 genes, left panels; see Table S4) or upregulated upon BST (783 genes, right panels; see Table S4). Numbers indicate the proportion (%) of genes having at least one ATAC peak with decreased (blue) or increased (red) accessibility in the proximity.

(C) Proportion (%) of genes having at least one ATAC peak enriched for the AP-1 motif, with decreased or increased accessibility in the proximity. p values calculated using proportional Z test. ***p < 0.001.

(D) c-FOS protein expression indicated by immunofluorescence staining on ASZ cells upon Ras/MAPK activation and/or TGF-β inhibition. Scale bars: 25 μm. c-FOS protein expression indicated by western blot is shown in Figure S3D.

(E) Quantification of c-FOS protein expression shown in (D) by pixel intensity measurements (n ≥ 40 cells, measured in n ≥ 3 fields for each condition). Horizontal bars and error bars represent the mean ± SD. p values calculated using one-way ANOVA with Tukey’s post-test. ***p < 0.001.

(F) Venn diagram showing the number of peaks (and GREAT-associated genes) that overlap between c-FOS ChIP-seq peaks and chromatin-accessible peaks enriched for AP-1 motif. Identified peaks coordinates and GREAT-associated genes are shown in Table S5.

(G) Proportion (%) of down- versus upregulated genes having a c-FOS ChIP binding site in the proximity. p value calculated using proportional Z test. ***p < 0.001.

(H) Relative enrichment of the AP-1-associated genes with or without overlapping c-FOS ChIP-seq peaks (1589 and 7230 genes, respectively) in mouse BCC (blue, n = 4) and SCC (orange, n = 7) bulk RNA-seq. Boxes represent the mean and distribution of values from minimum to maximum. p values calculated using unpaired two-sided Student’s t test. ***p < 0.001.

AP-1 TFs encompass the JUN family (JUN, JUNB, JUND), the FOS family (c-FOS, FOSB, FOSL1, FOSL2), and the related ATF and MAF subfamilies, with their transcriptional regulatory activity being highly dependent on the homo- or hetero-dimeric combinations formed (Lopez-Bergami et al., 2010). This led to our investigation of which AP-1 family members mediate the differential BCC- and SCC-associated chromatin accessibility switch. Importantly, our previous work indicates that JUN family members are highly expressed in human BCCs (Maglic et al., 2018), and we have shown that JUN cooperativity in the absence of c-FOS drives BCC resistance through non-canonical HH pathway activation (Yao et al., 2020). Interestingly, with their oncogenic roles reported in epidermal tumors (Angel et al., 2001; Briso et al., 2013; Saez et al., 1995), FOS family members are induced upon BST in BCC tumor cells (Figures S3C and S3D). c-FOS, in particular, has been reported as a master regulator of SCC identity and progression (Guinea-Viniegra et al., 2012), supporting its potential role in driving BST. Indeed, c-FOS accumulates in the nucleus upon experimental BST, where it interacts with both c-JUN and JUNB (Figures 3D, 3E, and S3E). To elucidate direct c-FOS transcriptional targets, we performed c-FOS chromatin immunoprecipitation sequencing (ChIP-seq) upon Ras/MAPK activation and TGF-β inhibition and found enrichment for the AP-1 motif in promoter, intergenic, and intronic regions (Figures S3F and S3G). Overlapping c-FOS ChIP-seq peaks with chromatin-accessible peaks enriched for AP-1 motif allowed the identification of 2412 high-confidence c-FOS binding sites and 1589 genomic regions enrichment of annotations tool (GREAT)-associated genes (Figure 3F; Table S5). We found a higher representation of the identified c-FOS target genes in UP than in DN genes upon BST induction (Figures 3B–3G). Among the AP-1-associated genes identified by ATAC-seq (8819 genes in total), 1589 genes annotated as c-FOS direct targets demonstrate enrichment in SCC tumors, while the remaining 7230 genes with no c-FOS binding in the proximity are enriched in BCC tumors (Figure 3H). Altogether, we conclude that BST involves the distinct regulation of AP-1-dependent BCC and SCC transcriptomic programs and that c-FOS appears to play a direct role in toggling those programs.

c-FOS drives BST in vitro and in vivo

We interrogated the role of c-FOS in BST in vitro and in vivo. First, we used transient silencing (siRNA) for c-Fos and achieved a partial rescue in both the BST-associated reduction in Gli1 and the BST-associated increase in Cd44 (Figure S4A). Next, we assayed c-FOS sufficiency in driving BST in vitro by stably transfected BCC cells with either empty or doxycycline-inducible c-FOS constructs (Figures S4B and S4C) and observed reduced nuclear GLI1 and primary cilia along with increased CD44 and MUC1 expression upon c-FOS induction (Figures S4D and S4E). In vivo, c-FOS induction in BCC xenografted tumors drove BST as well, as shown by paler and larger nuclei, eosinophilic cytoplasmic changes, keratohyalin granules, keratin pearls, reduced nuclear GLI1, primary cilia, and increased CD44 and MUC1 expression, closely resembling Ras/MAPK activation and/or abrogation of TGF-β signaling (Figures 4A and 4B). We further validated c-FOS sufficiency in primary, unselected BCCs by enforcing c-FOS expression in Ptch+/− ;K14-Cre-ER;p53;RFPfl-stop-fl RFPpos BCCs (Wang et al., 2011a). We performed intratumoral injection of AAV6_GFP_cFos or AAV6_GFP viruses and identified infected tumor epithelial cells (RFPposGFPpos) by flow cytometry and immunofluorescence (Figures S4F–S4I). As with BCC cells, GFP_cFos, but not GFP_control-infected areas, showed squamous differentiation with reduced nuclear GLI1 and increased CD44 (Figures 4C–4E and S4J). Importantly, RNA-seq and ATAC-seq upon c-FOS induction confirmed reduced enrichment for the BCC signature and increased enrichment for the SCC signature (Figures S4K, S4L, and 4F), as well as significant chromatin remodeling in regions enriched with hair follicle and epidermal TF motifs, as observed upon Ras/MAPK activation and/or TGF-β inhibition (Figures 4G, S4M, and S4N). In particular, AP-1 appeared again among the top-enriched motifs in remodeled chromatin during c-FOS-mediated BST. Altogether, these in vitro and in vivo studies confirm c-FOS as a master regulator of the epigenetic and transcriptomic reprogramming associated with BST.

Figure 4.

The AP-1 family member c-FOS drives BST

(A) H&E, GLI1, ARL13B (for primary cilia), CD44, and MUC1 staining of ASZ xenografted tumors upon c-FOS induction.

(B) Quantification of GLI1, CD44, and MUC1 protein expression in (A) by pixel intensity measurements (n ≥ 40 cells, measured in n ≥ 3 fields, across n ≥ 3 tumors for each condition). Quantification of primary cilia (ARL13B) in (A) by the percentage of ciliated tumor cells (measured in n ≥ 6 fields, across n ≥ 3 tumors for each condition).

(C) Schematic representation of in vivo intratumoral adenoviral-based induction of c-FOS using Ptch+/− ;K14-Cre-ER;p53;RFP mice. Lower panels show H&E, GFP, GLI1, and CD44 staining of BCC tumors upon vehicle, AAV6_GFP, or AAV6_GFP_cFos injection. Infection rate was determined as shown in Figures S4F, S4G, and S4I.

(D) Quantification of SCC-like foci density shown in (C) (number of fields measured - n ≥ 10, across n = 2 tumors for each condition).

(E) Quantification of GLI1 and CD44 protein expression shown in (C) by pixel intensity measurements (n ≥ 40 cells, measured in n ≥ 3 fields, across n = 2 tumors for each condition).

(F) Relative enrichment for BCC signature (blue) and SCC signature (red) in RNA-seq data obtained from ASZ upon c-FOS induction (n = 2 biological replicates). RNA-seq profile and GSEA plots are shown in Figures S4K and S4L.

(G) Enrichment of TF motifs in the regions with decreased (left) or increased (right) chromatin accessibility in ASZ upon c-FOS induction. The y axis is −log10(p value) of a motif enrichment score. The x axis is the ranking number of sorted motifs. Motifs belonging to one TF family are labeled with the same color. ATAC-seq profile is shown in Figure S4M.

(H) scRNA-seq projection plot of tumor cells obtained from both nBCCs and rBCCs (shown in Figure 1A) analyzed for SCC (red) versus BCC (blue) and expBST UP (red) versus expBST DN (blue) signatures, respectively (upper panels). Violin plots show the signature enrichment in individual clusters. Clusters 6 (pink) and 7 (purple) are highlighted as representing opposed transcriptomic profiles. p values calculated using two-sided Wilcoxon rank sum test. Definition and composition of expBST-UP (red) and expBST-DN (blue) signatures are shown in Figure S5E and Table S6. Inflammatory response enrichment in individual clusters is shown in Figure S5G.

Scale bars in (A) and (C): 100 μm. Horizontal bars and error bars in (B) and (D)–(F) represent the mean ± SD. p values in (B) and (F) calculated using unpaired two-sided Student’s t test. p values in (D) and (E) calculated using one-way ANOVA with Tukey’s post-test. *p < 0.05, ***p < 0.001.

c-FOS and surrogates of c-FOS activity also correlated with BST in patient tumor samples. Patient tumors revealed increased expression of c-FOS in BSCs and well-differentiated SCCs (Figures S5A–S5C), correlating with reduced nuclear GLI1 and increased membranous CD44 (Figure S5D). Moreover, we established an experimental BST signature based on the concordantly UP and DN gene sets in our BST models (expBST UP; expBST DN; Figure S5E; Table S6) and validated it in bulk RNA-seq profiles of mouse and human tumors (Figure S5F). We then overlayed the expBST signatures on patient scRNA-seq data pooled from nBCC and rBCC tumor cells (shown in Figure 1A) and found that the UP and DN genes in the expBST signatures significantly overlapped with SCC and BCC signature enrichment, respectively (Figure 4H, upper panels). In particular, when comparing the opposed “SCC-like” cluster 7 to the “BCC-like” cluster 6, we observed a concomitant decrease in expBST DN and an increase in expBST UP signatures (Figure 4H, lower panels). Consistent with a previously reported role for epidermal c-FOS in mediating inflammation (Briso et al., 2013), we observed a concomitant increase in inflammatory (Hallmark_inflammatory_response) and BST (expBST_UP) signatures (Figure S5G). Altogether, we conclude that c-FOS is a key functional driver of BST in both human and mouse BCC tumors.

EGFR/Ras/MAPK inhibitors partially reverse the transcriptional reprogramming of c-FOS-mediated BST

Our data support BST as an epigenetic, rather than a genetic, reprogramming mechanism requiring continuous c-Fos activity. Intriguingly, RNA velocity analysis on BCC-RMpos clusters 3 and 7 reveals bidirectional transcriptional transitions (Figure 2B), supporting the notion of BST transcriptional reversibility (Guo et al., 2018; Haensel et al., 2020). To test this hypothesis, we first used our in vitro BST system and assessed BCC and SCC markers upon 6 days of c-FOS induction (Dox ON), followed by 6 days of c-FOS removal (Dox OFF) (Figure 5A). Remarkably, c-FOS removal reversed the marked transcriptional and morphological changes associated with BST, including restoration of the primary cilia, upregulation of GLI1, and reduced CD44 transcription (Figures 5A and S6A), supporting the bidirectional nature of BST.

Figure 5.

EGFR/Ras/MAPK inhibitors partially reverse the transcriptional reprogramming of c-FOS-mediated BST

(A) DAPI, GLI1, ARL13B (for primary cilia), and CD44 staining of ASZ cells upon 6 days of c-FOS induction (Dox ON), followed by 6 days of c-FOS removal (Dox ON/OFF). The right panels show the quantification of GLI1 and CD44 protein expression by pixel intensity measurements (n ≥ 40 cells, measured in n ≥ 3 fields for each condition) and primary cilia (ARL13B) by the percentage of ciliated tumor cells (measured in n ≥ 6 fields for each condition).

(B) Toplist of cancer-related pathways or biological processes identified using GSEA in expBST UP signature. The x axis is −log10(p value) of the enrichment score.

(C) Proportion (%) of genes having at least one ATAC peak enriched for the AP-1 or ETS motif, with decreased or increased accessibility in the proximity. p values calculated using proportional Z test. **p < 0.01.

(D) Schematic depicting the c-FOS-driven BST model.

(E) Dose-response curves of Gli1 and Cd44 expression in c-FOS-induced BST upon afatinib or UO126 treatment (n ≥ 4 biological replicates). EC50 (half maximal effective concentrations) are indicated in red dotted lines. The upper panel depicts the experimental design.

(F) Scatterplot of genes with differential expression upon c-FOS induction and afatinib or UO126 treatment (n = 2 biological replicates). Genes with expression fold change ≥ 1.5 and p < 0.05 are marked in red, while genes with expression fold change < −1.5 and p < 0.05 are marked in blue. The experimental design is shown in (E) (upper panel).

(G) Relative enrichment for expBST DN (blue) and expBST UP (red) signatures in RNA-seq data obtained from ASZ upon c-FOS induction and UO126 or afatinib treatment, as described in (E) and (F) (n = 2 biological replicates). Horizontal bars represent the mean.

(H) Heatmap and GSEA for expBST genes showing EGFR/MAPK dependency. Genes in the expBST DN signature are depicted in blue, while genes in the expBST UP signature are depicted in red and shown in Table S6.

(I) Scatterplot of peaks with differential accessibility upon c-FOS induction and afatinib or UO126 treatment (n = 2 biological replicates). Peaks with differential accessibility log2(fold change) ≥ 1 and p < 0.05 are marked in red, while peaks with differential accessibility log2(fold change) < −1 and p < 0.05 are marked in blue. The experimental design is shown in (E) (upper panel).

(J) c-FOS protein expression upon c-FOS induction and afatinib or UO126 treatment by immunofluorescence staining and quantification on ASZ cells. Quantification of c-FOS protein expression by pixel intensity measurements (n ≥ 40 cells, measured in n ≥ 3 fields for each condition).

(K) Chromatin immunoprecipitation and quantitative PCR (ChIP-qPCR) analysis of the DNA binding activity of c-FOS in the proximity of its target gene loci Fosl1 and Epgn upon c-FOS induction and afatinib or UO126 treatment (n = 3 biological replicates). Lower panels show the FPKM levels of Fosl1 and Epgn in similar conditions (n = 2 biological replicates).

Scale bars in (A) and (J): 25 μm. Horizontal bars and error bars in (A), (E), (G), (J), and (K) represent the mean ± SD. p values in (A), (G), (J), and (K) calculated using one-way ANOVA with Tukey’s post-test. ns, non-significant. *p < 0.05, **p < 0.01, ***p < 0.001.

We thus sought to target c-FOS-mediated BST in the absence of clinically available c-FOS inhibitors. Indeed, we previously showed the toxicity of the AP-1 inhibitor T-5224 for advanced BCC, despite the absence of c-FOS expression (Yao et al., 2020). To circumvent this issue, we sought to interrogate transcriptomic and chromatin remodeling data for clinically targetable pathways. Interestingly, GSEA for cancer-related pathways revealed significant enrichment of the epidermal growth factor receptor (EGFR)/Ras/MAPK signaling in expBST UP transcriptomic signature (Figure 5B). Moreover, ATAC-seq profiles identified significant ETS motif and AP-1 motif enrichment in open chromatin in the proximity of expBST UP genes (Figures S6B and 5C), suggesting a concomitant activation of the EGFR/Ras/MAPK/ETS signaling along with c-FOS induction. Indeed, we observed both EGFR (Y1068) and ERK1/2 phosphorylation upon c-FOS induction in our BST model (Figure S6C). EGFR phosphorylation (Y1068) further correlated with nuclear c-FOS in human BSCs and well-differentiated SCCs (Figures S6D and S6E). Examination of transitional BSC patient tumors revealed cellular colocalization of nuclear c-FOS, p-EGFR, and phospho-mitogen-activated protein kinase kinase (p-MEK), suggesting their concomitant activation (Figures S6F and S6G).

Based on these correlative observations, we tested the hypothesis that EGFR or MAPK inhibitors may revert c-FOS-induced BST (Figure 5D). We conducted an in vitro EGFR (afatinib) and MAPK (UO126) inhibition assay and found that short-term EGFR or MAPK inhibition (6 h) altered c-FOS-mediated BST in a dose-dependent manner, as shown by nuclear Gli1 rescue and Cd44 decrease, with clinically acceptable EC50 ranges (Figure 5E). Importantly, RNA-seq profiles confirmed that both inhibitors partially reverted BST at the transcriptomic level (Figures 5F and 5G), particularly for cell-differentiation-related gene sets (Figure 5H; Table S6). The BST transcriptomic reversal, however, does not come with epigenetic reversal in our experimental system, as suggested by the few chromatin accessibility changes observed upon both inhibitor treatments (Figure 5I). Consistently, c-FOS mRNA (Figure S6H) and protein (Figures S6I and 5J) expression levels were unaffected by the inhibitors, as well as its DNA binding in the proximity of EGFR-dependent genes (Fosl1 and Epgn) (Figure 5K). Altogether, our experimental system suggests that EGFR/MAPK inhibition partially reverses BST transcriptional output without affecting c-FOS expression level and DNA binding, which may expose to fast BST recurrence upon drug cessation.

Targeting EGFR/Ras/MAPK signaling pathway blocks BST

Based on the partial reversibility demonstrated by EGFR and MAPK inhibitors in vitro, we then investigated whether clinically approved afatinib could prevent BST in vivo. Remarkably, while afatinib had a negligible effect on the BCC phenotype in vivo (Figures 6A and 6B, left columns), it prevented c-FOS-induced BST as shown by restored dense nuclei, and basophilic cytoplasm reduced keratohyalin granules and keratin pearls, rescued nuclear GLI1 and primary cilia, and reduced CD44 and MUC1 expression (Figures 6A and 6B, right columns). Similar results were obtained with erlotinib, an additional clinically approved EGFR inhibitor (EGFRi) (data not shown). Encouraged by these promising results, we assessed prevention of BST in patient tumors. Remarkably, tumors that harbor BST histological features and intermediate nuclear c-FOS expression specifically achieved both GLI1 rescue and CD44 decrease upon short-term ex vivo afatinib treatment, while classical BCCs or well-differentiated SCCs did not or did only partially (Figure 6C). Altogether, our data indicate that c-FOS-mediated BST may be partially reverted using EGFR inhibitors.

Figure 6.

Targeting EGFR/Ras/MAPK signaling pathway blocks BST

(A) H&E, GLI1, ARL13B (for primary cilia), CD44, and MUC1 staining of ASZ xenografted tumors upon c-FOS induction and treatment with afatinib.

(B) Quantification of GLI1, CD44, and MUC1 protein expression shown in (A) by pixel intensity measurements (n ≥ 40 cells, measured in n ≥ 3 fields, across n ≥ 3 tumors for each condition). Quantification of primary cilia (ARL13B) shown in (A) by the percentage of ciliated tumor cells (number of fields measured: n ≥ 6, across n ≥ 3 tumors for each condition). Horizontal bars and error bars represent the mean ± SD. p values calculated using one-way ANOVA with Tukey’s post-test. ns, non-significant. *p < 0.05, ***p < 0.001.

(C) GLI1 and CD44 expression fold change upon afatinib compared to DMSO treatment in human tumor explant (lower panels) and related to histopathological diagnosis (BCCs in blue, n = 6; BSCs in green, n = 6; and well-differentiated SCCs in orange, n = 4) and nuclear c-FOS expression (number of fields measured: n ≥ 2 per tumor). Upper panels show representative H&E, DAPI, and c-FOS staining of the human tumor explants used for the experiment. qRT-PCR error bars represent the SD from two technical replicates per sample. Immunofluorescence error bars represent the SD from individual measurements across one tumor sample. p values calculated using unpaired two-sided Student’s t test. ns, non-significant. *p < 0.05.

Scale bars in (A) and (C): 100 μm.

DISCUSSION

In contrast to classical binary “steady states” of BCC and SCC, single-cell and multi-dimensional analysis reveal BST’s continuous and reversible nature. The present study and recent work from our group (Yao et al., 2020) support BCC evolution from SMOi-sensitive BCCs to rBCCs and from rBCCs to BSCs and well-differentiated SCCs. By integrating chromatin accessibility profiles with transcriptomic profiles in our BST models, we identify the pivotal role of c-FOS in toggling AP-1 chromatin accessibility and transcriptomic reprogramming from non-overlapping BCC-specific AP-1 sites to SCC-specific AP-1 sites. Previous studies identified a high frequency of epidermal genes with AP-1 motifs near their promoter regions (Eckert and Welter, 1996; Hu and Gudas, 1994; Navarro et al., 1995), suggesting a dynamic and sequential regulation of AP-1-dependent effector genes along epidermal differentiation trajectories (Angel et al., 2001) that have been repurposed in skin tumor resistance. In accordance with this notion, we recently demonstrated that JUN family members induce a chromatin accessibility profile in response to p38 JUN kinase/TGF-β signaling that drives SRF-MRTF-GLI chromatin association toward non-canonical HH pathway activation and SMOi resistance in BCCs (Blatti and Scott, 1992; Chung et al., 1996; Virolle et al., 1998; Yamamura et al., 2000; Yao et al., 2020; Zhang et al., 1998). By contrast, c-FOS-induced BST downstream of TGF-β inhibition and Ras activation induces the SRF-AP-1-ETS-associated transcriptomic repertoire and the switch to epidermal (Fischer et al., 1999; Sark et al., 1998) and SCC superenhancer-associated AP-1 and ETS binding sites (Yang et al., 2015).

AP-1 members are recognized as pioneer TFs establishing cell-type-specific enhancers and chromatin architecture (Heinz et al., 2010; Iwafuchi-Doi and Zaret, 2014; Madrigal and Alasoo, 2018; Phanstiel et al., 2017). Our work highlights how AP-1 isoform-specific toggling dictates cell-type-specific chromatin remodeling and transcriptional repertoires necessary for tumor lineage reprogramming. We initially identified TACSTD2, LYPD3, and LY6D BCC-RMs as defining the progression of SMOi-sensitive to JUN-dependent rBCCs (Yao et al., 2020). Our study demonstrates that the different cell states along AP-1 isoform-specific tumor lineage reprogramming share the expression of the previously identified BCC-RMs. Furthermore, the epigenetic proximity of BCC-RMpos resistant cells with suprabasal cells of the hair follicle (Yao et al., 2020), as well as their gradual increase along SCC differentiation in human tumors, suggests that they encompass heterogeneous morphological and transcriptional cell states of keratinocyte commitment, whose regulation is currently under investigation. From a clinical point of view, our work further supports that BCC-RMs, as surrogate markers for c-FOS-mediated BST propensity, may serve as useful markers to predict resistance to HH pathway inhibitors and poor prognosis (Tan et al., 2017).

Because of the central role c-FOS plays in driving BST, both genetic and epigenetic regulators of c-FOS activity will define a tumor’s propensity to undergo BST. Arguing for a genetic basis is the finding that BCC-RMpos clones rarely occur in Gorlin’s syndrome BCCs, suggesting that accumulated mutational burden and/or genomic instability may predispose patients to BST (Bonilla et al., 2016; Yao et al., 2020). Indeed, recurrent mutations in ciliome genes that reduce HH signaling increase c-FOS expression (Kuonen et al., 2019). Additionally, spontaneous BSC and pancreatic tumors with squamous features contain mutations in chromatin regulators like ARID1A (Chiang et al., 2019; Hayashi, 2020), which may favor c-FOS-induced chromatin remodeling (Sun et al., 2016). The overall genetic similarity between BSC and BCC (Chiang et al., 2019) and the reversibility of BST suggest additional non-genetic contributors to BST. The tumor microenvironment is a critical contributor to tumor cell plasticity in different cancer models (Hirata et al., 2015; Kaur et al., 2016; Lee et al., 2011; Roswall et al., 2018; Smith et al., 2014; Straussman et al., 2012), and inflammation, in particular, has been shown to regulate AP-1 chromatin accessibility (Naik et al., 2017). In the skin, chronic inflammation drives squamous metaplasia (Na et al., 2011; Trikudanathan and Dasanu, 2010) and preferential SCC, rather than BCC, development (Xiang et al., 2019). In the era of immunotherapies, and given the reported c-FOS-mediated inflammation in the skin (Briso et al., 2013), more sophisticated immunocompetent BST experimental models will be required to address the exact role of the immune system during c-FOS-mediated BST. Our results suggest that further insights into the genetic background and environmental influences that enforce c-FOS expression through epithelial EGFR/Ras/MAPK or TGF-β signaling pathways will help predict BST and therapy efficacy for BCC and BSC (Bhowmick et al., 2004; Pastore et al., 2008). Intriguingly, a significant proportion of SCCs harbor inactivating mutations in Ptch1 (Chiang et al., 2019; Ping et al., 2001), suggesting their genetic derivation from BCC. However, additional studies are needed to better characterize the genetic evolution of SCC, its heterogeneity, and how it may dictate progression and therapeutic sensitivities.

Motivated by BST’s apparent reversibility from our single-cell analysis, we identified and provided initial evidence for the clinically approved EGFRi afatinib in reversing BST as defined by GLI1, CD44, and cilia restoration. Human primary BCC tumor explants with SCC features (BSCs) demonstrated a better response to afatinib treatment as compared to nBCCs or SCCs, highlighting the EGFR/Ras/MAPK/ETS axis as a potential therapeutic target to revert c-FOS-driven BST. Our data predict that by restoring GLI1 and abrogating CD44, afatinib may resensitize tumors to HH inhibition, with SMO blockade in the correct genetic background sensitizing tumors to EGF-Ris. Indeed, our data provide an explanation for monotherapy failure as an alternative population arises and provides a mechanism for the observed synergistic effect of SMOi and EGFRi combination therapy in pancreatic (Hu et al., 2007), cerebellar (Götschel et al., 2013), lung (Della Corte et al., 2015), and head and neck tumors (Liebig et al., 2017). While we observed reversal of key BCC markers upon afatinib or UO126 treatment, supporting their role in arresting BST, detailed transcriptomic analyses suggest these drugs induce a partial transcriptional reversal without reversing chromatin accessibility or c-FOS binding to DNA. Therefore, we may expect BST to easily recur upon drug cessation. The data raise interesting questions regarding the role of tyrosine kinase-independent EGFR signaling (Eldredge et al., 1994; Han and Lo, 2012) or EGFR-independent signaling in c-FOS-mediated BST that may limit the efficiency of EGFR tyrosine kinase inhibitors, highlighting the need for the development of c-FOS and AP-1 isoform-specific inhibitors in the future. Further, the current study focusing on the transition of advanced rBCCs to well-differentiated SCCs by c-FOS complements our previous work demon-strating that JUN inhibitors can restore SMOi sensitivity (Yao et al., 2020) and suggests that additional in vitro and in vivo elucidation of BST signaling pathways will uncover a panoply of additive combination therapies that will further improve cancer therapy.

STAR★METHODS

RESOURCE AVAILABILITY

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Anthony Oro (oro@stanford.edu).

Materials availability

All materials generated in this study are available upon request.

Data and code availability

Sequencing data generated for this manuscript have been deposited at the Gene Expression Omnibus and are publicly available as of date of publication. The accession number is listed in the Key Resources Table. The accession numbers of existing, publicly available sequencing data analyzed in this paper are listed in the Key Resources Table.

KEY RESOURCES TABLE

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
ARL13B	Abcam	Cat#ab136648
GLI1	Novus	Cat#NBP1-78259; RRID:AB_11030198
CD44	Invitrogen	Cat#14-0441-82; RRID:AB_467246
MUC1	Abcam	Cat#ab15481; RRID:AB_301891
Phospho-MEK1	Abcam	Cat#ab96379; RRID:AB_10678572
Phospho-SMAD3	Abcam	Cat#ab52903; RRID:AB_882596
SMAD3	Abcam	Cat#ab40854; RRID:AB_777979
Phospho-EGFR	Abcam	Cat#ab40815; RRID:AB_732110
TACSTD2	Abcam	Cat#ab214488; RRID:AB_2811182
GFP	Abcam	Cat#ab13970; RRID:AB_300798
RFP	Abcam	Cat#ab124754; RRID:AB_10971665
LYPD3	Sigma-Aldrich	Cat#HPA041529; RRID:AB_2677534
LY6D	Proteintech	Cat#17361-1-AP; RRID:AB_2878397
C-FOS	LS-Bio	Cat#LS-B14369
MYC-tag	CST	Cat#2276; RRID:AB_2148465
Phospho-ERK1/2	CST	Cat#4377; RRID:AB_331775
ERK1/2	CST	Cat#9102; RRID:AB_330744
C-FOS	CST	Cat#2250; RRID:AB_2247211
B-TUBULIN	DSHB	Cat#AB_2315513; RRID:AB_528499
IRDye 680LT Donkey anti-rabbit IgG	LI-COR	Cat#926-68023; RRID:AB_10706167
IRDye 800CW Donkey anti-mouse IgG	LI-COR	Cat#926-32212; RRID:AB_621847
Donkey anti-rabbit IgG (Alexa Fluor 488)	ThermoFisher	Cat#A21206; RRID:AB_2535792
Donkey anti-mouse IgG (Alexa Fluor 555)	ThermoFisher	Cat#A31570; RRID:AB_2536180
Donkey anti-rat IgG (Alexa Fluor 555)	ThermoFisher	Cat#A48270
Donkey anti-goat IgG (Alexa Fluor 680)	ThermoFisher	Cat#A21084; RRID:AB_2535741
Rabbit IgG	CST	Cat#2729; RRID:AB_1031062
C-JUN	CST	Cat#9165; RRID:AB_2130165
Fluorescent TrueBlot®: Anti-Rabbit IgG DyLight 680	ROCKLAND	Cat#18-4416-32
FRA1	Abcam	Cat#ab252421
FOSB	Abcam	Cat#ab184938; RRID:AB_2721123
JUNB	Abcam	Cat#ab128878; RRID:AB_11144334
JUND	Abcam	Cat#ab181615; RRID:AB_2864350
FRA2	Abcam	Cat#ab124830; RRID:AB_11128294
Bacterial and virus strains
Virus: AAV6-CMV-CFOS-IRES-GFP	Vector Biosystems	Cat#AAV-229259
Virus: AAV6-CMV-GFP	Neurosciences Institute, Stanford	Cat#GVVC-AAV-5
Biological samples
Human BCC, BSC and SCC sections	Department of Dermatology, Stanford University	https://www.stanford.edu/
Human BCC, BSC and SCC tumor explants	Department of Dermatology, UNIL-CHUV	https://chuv.ch/
Lipofectamine LTX with Plus Reagent	ThermoFisher	Cat#A12621
Lipofectamine RNAiMax Reagent	ThermoFisher	Cat#13778100
Matrigel Matrix, Basement Membrane, Growth Factor reduced	FisherScientific	Cat#356238
Chemicals, peptides, and recombinant proteins
Afatinib	Selleckchem	Cat#S1011
Erlotinib	Selleckchem	Cat#S7786
UO126	Selleckchem	Cat#1102
Vismodegib	Selleckchem	Cat#1082
XAV-939	Selleckchem	Cat#1180
RO4929097	Selleckchem	Cat#1575
SB431542	Selleckchem	Cat#1067
Avanafil	Selleckchem	Cat#4019
Mouse recombinant TNFα	Biolegend	Cat#575202
Phorbol 12-myristate 13-acetate (PMA, TPA)	Sigma-Aldrich	Cat#P8139
4-OHT	Sigma-Aldrich	Cat#H7268
Tamoxifen	Sigma-Aldrich	Cat#T5648
Doxycycline	Sigma-Aldrich	Cat#D9891
Hoechst 33342	ThermoFisher	Cat#H3570
Critical commercial assays
MinElute PCR Purification kit	QIAGEN	Cat#28004
RNeasy Plus Mini kit	QIAGEN	Cat#74136
Nextera DNA Library Prep Kit	Illumina	Cat#FC-121-1030
Illumina Tagment DNA TDE1 Enzyme and Buffer Kit	Illumina	Cat#20034197
CUT&RUN Assay Kit	CST	Cat#86652
Pierce c-Myc-Tag Magnetic IP/Co-IP Kit	ThermoFisher	Cat#88844
Deposited data
Raw RNaseq and ATACseq	This paper	GEO:GSE168376
Naive BCC scRNaseq	Yao et al., 2020	GEO:GSE141526
Resistant BCC scRNaseq	Yost et al., 2019	GEO:GSE123813
Bulk human SCC RNaseq	Chitsazzadeh et al., 2016	GEO: GSE84194
Bulk human normal skin and BCC RNaseq	Bonilla et al., 2016	EGA: EGAS00001001540
Bulk mouse normal skin and SCC RNaseq	Nassar et al., 2015	ArrayExpress: E-MTAB-2889
Bulk mouse BCC RNaseq	Whitson et al., 2018	GEO: GSE78497
Experimental models: Cell lines
ASZ_001	Gift of Dr. E. Epstein	N/A
Experimental models: Organisms/strains
Mouse: NOD-SCID	Jackson Laboratory	https://www.jax.org/
Mouse: Ptch1+/−K14-Cre-ER2 p53fl/fl RFPfl/stop/fl	Gift of Dr. E. Epstein	N/A
Oligonucleotides
6-FAM non silencing siRNA	Sigma-Aldrich	Cat#SIC007
Mouse c-FOS silencing si RNA: mFOS_siRNA_01	Sigma-Aldrich	Cat#SASI_Mm01-00192758
Mouse c-FOS silencing si RNA: mFOS_siRNA_02	Sigma-Aldrich	Cat#SASI_Mm01_00192759
GLI1-FAM: Hs00171790 (TaqMan Gene Expression Assays)	ThermoFisher	N/A
CD44-FAM: Hs05662929 (TaqMan Gene Expression Assays)	ThermoFisher	N/A
KRT14-FAM: Hs00265033 (TaqMan Gene Expression Assays)	ThermoFisher	N/A
GAPDH-FAM: Hs02786624 (TaqMan Gene Expression Assays)	ThermoFisher	N/A
GLI1-FAM: Mm00494654 (TaqMan Gene Expression Assays)	ThermoFisher	N/A
CD44-FAM: Mm01277165 (TaqMan Gene Expression Assays)	ThermoFisher	N/A
MUC1-FAM: Mm00449604 (TaqMan Gene Expression Assays)	ThermoFisher	N/A
CFOS-FAM: Mm00487425 (TaqMan Gene Expression Assays)	ThermoFisher	N/A
GAPDH-FAM: Mm99999915 (TaqMan Gene Expression Assays)	ThermoFisher	N/A
Recombinant DNA
Piggybac vector	Gift of Dr. S. Yamanaka	N/A
ERRasV12	Gift of Dr. P. Khavari	N/A
Mouse myc-DYK-tagged c-OFS ORF	Origene	Cat#MR205933
MISSION Mouse TGFβRI_shRNA-pLKO.1-puro	Sigma-Aldrich	Cat#SHCLNG-NM_009370
MISSION Non-mammalian shRNA control-pLKO.1-puro	Sigma-Aldrich	Cat#SHC002
pMD2G	Addgene	Cat#12259
psPAX2	Addgene	Cat#12260
Software and algorithms
Prism	Graphpad	Version 8.0
ImageJ	National Institue of Health	Version 2.0.0
Seurat	Satija et al., 2015	Version 3.2.1
VelocytoR	La Manno et al., 2018	Version 0.6
TopHat	https://ccb.jhu.edu	Version 2.1.1
Homer	Univeristiy of California San Diego	Version 4.10
DESeq2	Bioconductor	N/A
EdgeR	Bioconductor	N/A
EnrichR	Maayan laboratory	maayanlab.cloud/Enrichr
GSEA	Broad Institute	Version 3.0
ssGSEA	GenePattern	Version 14; https://www.genepattern.org:443/
Bowtie	http://bowtie-bio.sourceforge.net	Version 1.1.2; Version 2.3.4.1
Samtools	http://www.htslib.org	Version 0.1.18
MACS2	Zhang et al., 2008	Version 2.2.5
GREAT	McLean et al., 2010	Version 4.0.4
FlowJo	BD Biosciences	Version 10.6
Image Studio Lite	LI-COR	Version 4.0
CHIPseekerR	Yu et al., 2015	N/A
Integrative Genomics Viewer	Broad Institute	Version 2.6.2
limma	Bioconductor	N/A

This paper does not report original code.

Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS

Animals

Mice were housed under standard conditions, and animal care was in compliance with the protocols approved by the Institutional Animal Care and Use Committee (IACUC) at Stanford University. Ptch1+/−K14-Cre-ER2 p53fl/fl RFPfl/stop/fl mice (C57BL/6J background) were generated passaged and utilized to develop BCC tumors as described previously (Wang et al., 2011a, 2011b). Here we irradiated mice (5Gy) using an X-ray irradiator. Adult male (5–7 weeks of age) NOD-SCID mice were used as host animals for grafted tumors. Littermate of the same sex were randomly assigned to experimental groups.

Human subjects

Human studies were approved by Stanford University (Stanford, CA, IRB #18325) and Lausanne University (Lausanne, Switzerland, CIR #102995) institutional review boards and were performed in accordance with the guidelines of the Declaration of Helsinki. Histological diagnoses of BCC, BSC, well-differentiated and poorly-differentiated SCC samples used in this study were confirmed by an independent dermatopathologist. Freshly resected tumors were obtained from patients undergoing dermatological surgery, with previous obtention of written informed consent for all patient samples. Samples were deidentified and there was no selection for age, gender, or other characteristics.

Cell lines

ASZ_001 cells (mouse BCC cell line, referred to as ASZ) were grown in M154CF (Thermo Fisher Scientific, Waltham, MA) supplemented with 2% chelated FBS, 1% Penicillin-Streptomycin and 0.05mM CaCl2 (So et al., 2006). Experiments were carried out in serum-starved conditions. Cell lines were assessed for mycoplasma using the MycoSeq Mycoplasma real-time PCR detection kit (ThermoFisher Scientific, Waltham, MA).

METHOD DETAILS

DNA constructs and transfection

The ERRasV12 sequence was inserted into the piggyback (pb) vector (courtesy of Yamanaka laboratory) using EcoRI and PacI restriction enzyme sites to generate a doxycycline-inducible expression of ERRasV12 as previously described (Kuonen et al., 2019). The mouse myc-DYK-tagged c-FOS ORF clone was obtained from Origene (Rockville, MD) and subcloned into the piggyback (pb) vector using EcoRI and NotI restriction enzyme sites to generate a doxycycline-inducible expression of c-FOS. The pb-ERRasV12 and pb-c-FOS plasmids were both transfected with a transposase-expressing plasmid (System Biosciences) into the indicated cell lines using Lipofectamine LTX and Plus reagents obtained from Thermo Fisher Scientific (Waltham, MA). Clones with stable insertion of the transposon were then selected in the medium containing the appropriate antibiotic for at least 48 hours.

Lentiviral constructs

Murine TGFβR1 silencing and non-silencing shRNA sequences cloned into the pLKO.1-puro vector were purchased from Sigma-Aldrich (St-Louis, MO). All sequences are available online (https://www.sigmaaldrich.com/). Lentiviruses were produced in 293T cells using Lipofectamine LTX and Plus protocol (Thermo Fisher Scientific, Waltham, MA) by co-transfecting the pLKO.1-puro constructs with pMD2G and psPAX2 plasmids. Infection of targeted cells was done by 1 hour of centrifugation in the presence of 8 μg/ml poly-brene. Selection was performed 48 hours after infection using 5 μg/ml puromycin (Sigma-Aldrich, St-Louis, MO).

Mouse model and drug treatment

Adult male (5–7 weeks of age) NOD-SCID mice were used as host animals for grafted tumors. Primary tumors were initiated by injection of ASZ (5×105 cells/mouse) tumor cells together with mouse dermal fibroblasts (1.0×106 cells/mouse) subcutaneously in 100 μL of PBS:Matrigel at a 1:1 ratio as previously described (Kuonen et al., 2018). Doxycycline 2mg/ml was provided in the drinking water with 2% sucrose. 2mg of tamoxifen (Sigma-Aldrich, St-Louis, MO) diluted in corn oil at 10mg/ml was administered by oral gavage. 25mg/kg of afatinib (Selleckchem, Houston, TX) diluted in 10% DMSO was administered by intraperitoneal injections. Drug treatments (doxycycline, tamoxifen, afatinib) were administered on a daily basis from day 3 after tumor inoculation until sacrifice. BST was assessed on HE staining.

siRNA knockdown

Indicated cell lines were transfected with 30 pmol of either a 6-FAM non-silencing siRNA or mouse c-FOS-targeting siRNA (Sigma-Aldrich, St-Louis, MO) using Lipofectamine RNAiMAX transfection reagents (Thermo Fisher Scientific, Waltham, MA) with standard protocol. siRNAs references are listed as follows:

6-FAM non-silencing siRNA SIC007

mFos_siRNA_01 SASI_Mm01_00192758

mFos_siRNA_02 SASI_Mm01_00192759

Transfection efficiency was evaluated by 6-FAM fluorescence. 24h after transfection, cells were starved in medium containing vehicle or doxycycline and 4-OHT for 24h before mRNA extraction. Knockdown efficiency was measured by qRT-PCR.

Immunofluorescence

GLI1 (NBP1-78259, 1:200), ARL13B (ab136648, 1:500), CD44 (14-0441-82, 1:100), MUC1 (ab15481, 1:100), p-MEK (ab96379, 1:250), p-SMAD3 (ab52903, 1:100), c-FOS (LS-B14369, 1:100), p-EGFR (ab40815, 1:100), TACSTD2 (ab214488, 1:150), LYPD3 (HPA041529, 1:100), LY6D (17361-1-AP, 1:100), GFP (ab13970, 1:500), RFP (ab124754, 1:100) and nuclei (Hoechst 33342, 1:2000) were stained using a standard immunofluorescence protocol for FFPE. Antigen retrieval was performed in pH 6.0 citrate buffer (Vector Laboratories, Burlingame, CA). For in vitro immunofluorescence, cells plated in 8-well chamber slides (Millipore) were fixed with 4% formaldehyde before standard immunofluorescence protocol. The fluorescent-labeled secondary antibodies were used as follows: anti-rabbit Alexa Fluor 488 (1:500, Life Technologies, A-21206), anti-mouse Alexa Fluor 555 (1:500, Life Technologies, A-31570), anti-rat Alexa Fluor 555 (1:500, Life Technologies, A-48270) and anti-goat Alexa Fluor 680 (1:500, Life Technologies, A-21084). Confocal imaging was carried out using a Leica SP8 microscope equipped with an adjustable white light laser and hybrid detectors (Leica, Allendale, NJ). To quantify nuclear GLI1, membraneous CD44, cytoplasmic MUC1, cytoplasmic p-MEK, nuclear p-SMAD3, nuclear c-FOS and membraneous p-EGFR, pixel intensity was measured using ImageJ software (NIH; Bethesda, Maryland). Primary cilia (ARL13B) were quantified as previously described in Kuonen et al. (2019).

Real-time quantitative PCR

RNA samples were obtained from adherent cells using RNeasy Plus Mini Kit from QIAGEN according to manufacturer’s instructions. Quantitative real-time PCR was performed using TaqMan assay mouse-specific primers obtained from Thermo Fisher Scientific (Waltham, MA). Each reaction was performed in replicate and values were normalized to GAPDH.

Western blotting

Total cell lysates were resolved by SDS_PAGE and blotted onto a nitrocellulose membrane (Whatman, 0.45μm, Sigma-Aldrich, St-Louis, MO). Membranes were blocked for 30 min at room temperature (RT) with 5% BSA blocking buffer, and incubated with anti-MYC-Tag (CST 2276, 1:1000), anti-phospho-ERK1/2 (CST 4377, 1:750), anti-ERK1/2 (CST 9102, 1:1000), anti-phospho-SMAD3 (ab52903, 1:1000), anti-SMAD3 (ab40854, 1:1000), anti-c-FOS (CST 2250, 1:1000), anti-phospho-EGFR (ab40815, 1:1000), anti-c-JUN (CST 9165, 1:1000), anti-JUNB (ab128878, 1:1000), anti-JUND (ab181615, 1:1000), anti-FOSB (ab184938, 1:10000), anti-FRA1 (ab252421, 1:1000), anti-FRA2 (ab124830, 1:1000) and anti-β-tubulin (AB_2315513, 1:2000) antibodies overnight at 4°C. Membranes were visualized using LI-COR secondary antibodies and Image Studio Lite software (LI-COR Biosciences, Lincoln, NE).

Single cell RNA sequencing analysis

KRT14 expressing tumor epithelial cells were extracted from naive BCCs (nBCC; GSE141526; Yao et al., 2020) and resistant BCCs (rBCCs; GSE123813; Yost et al., 2019). For the nBCCs, we filtered out non-epithelial cells/clusters such as fibroblasts, immune cells, and endothelial cells through the markers COL1A2, PTPRC, and PECAM1, respectively. For analysis of the individual nBCC3 and nBCC8 samples, the top 5 PCs with a resolution of 0.3 were used. For the rBCCs, we extracted tumor epithelial cells by using the metadata included by Yost et al. which labeled the tumor epithelial cells (Yost et al., 2019). Tumor epithelial cells from the nBCC and rBCC datasets were merged together using the Multiple Dataset Integration and Label Transfer (anchoring) standard workflow (Stuart et al., 2019) and the top 15 PCs with a resolution of 0.4 were used. For gene scoring, the AddModuleScore function in Seurat R package was used (Seurat v3.2.1) (Satija et al., 2015). The two-sided Wilcoxon rank sum test was used to determine if there were significant differences between specified groups of cells.

For RNA Velocity analysis (La Manno et al., 2018), we utilized the SeuratWrappers R package for Velocity (SeuratWrappers v0.2.0; Velocyto.R v0.6). We superimposed the velocity results onto BCC3, using the cellular embeddings and parameters previously described. The neighborhood value for velocity projection was n = 1000.

RNA sequencing and analysis

RNA samples were obtained from adherent cells using RNeasy Plus Mini Kit from QIAGEN according to manufacturer’s instructions. Library preparation and sequencing were performed as described previously (Atwood et al., 2015). The RNaseq libraries were constructed by TruSeq Stranded mRNA Library Prep kit (Illumina) and sequenced on an Illumina NextSeq sequencer. Read alignment was performed using TopHat with mm9 as a reference genome. Raw counts and FPKM values were called using the HOMER function “analyzeRepeats.pl.” The DESeq2 R package was used to generate the differential gene expression in:

ASZ_NSshRNA_RasV12/Dox-4OHT versus ASZ_NSshRNA_RasV12

ASZ_TGFβRIshRNA_RasV12 versus ASZ_NSshRNA_RasV12

ASZ_TGFβRIshRNA_RasV12/Dox-4OHT versus ASZ_NSshRNA_RasV12

ASZ_cFos/Dox versus ASZ_empty/Dox

ASZ_cFos/Dox/DMSO versus ASZ_empty/Dox/DMSO

ASZ_cFos/Dox/UO126 versus ASZ_cFos/Dox/DMSO

ASZ_cFos/Dox/Afatinib versus ASZ_cFos/Dox/DMSO

Differential gene expression analysis was performed with cutoffs at log2 fold change > 1 or < −1 and P-value < 0.05 unless otherwise specified.

Mouse RNA-seq data generated in this manuscript were submitted to GEO (GSE168376). Publicly available RNA-seq data used from human (Bonilla et al., 2016; Chitsazzadeh et al., 2016) and mouse (Nassar et al., 2015; Whitson et al., 2018) skin, BCCs and SCCs were similarly aligned using TopHat with hg19 and mm9 as reference genomes respectively. These data are accessible using the GEO accession codes GSE84194 and GSE78497, the European Genome-phenome Archive identifier EGAS00001001540 and the ArrayExpress identifier E-MTAB-2889. The top 300 differentially upregulated and downregulated genes were used to define the human and mouse SCC and BCC signatures respectively. Gene signature enrichment in each sample was calculated based on normalized gene counts using ssGSEA from GenePattern (Broad Institute, Cambridge, MA) and expressed as the ranking difference between categories or ROAST gene set testing (Wu et al., 2010) and logFCs visualized using a barcode enrichment plot with limma package (Bioconductor). Gene lists were submitted to Enrichr (Ma’ayan laboratory, Icahn School of Medicine, New York, NY) for cancer- and development-related pathways enrichment. The read counts of differential gene lists were normalized using quantile normalization and the heatmaps were generated by clustering genes using the “heatmap.2” function from gplots R package.

ATAC sequencing library preparation, sequencing and analysis

ATAC-seq was performed as described previously (Buenrostro et al., 2013; Corces et al., 2017) with minor modifications for on-plate transposition in ASZ cells (Yao et al., 2020). In brief, 100,000 ASZ001 cells were plated per replicate in a 48-well plate, then serum starved for 24h prior to transposition. Wells were washed with cold PBS and lysed using 0.1% NP40 in resuspension buffer for 10 min at 25°C. Wells were washed again, then 100 μL of transposition mixture was added per well and incubated at 37°C for 30 min, shaking at 1000 RPM using a thermoshaker. DNA was purified using QIAGEN MinElute PCR Purification Kit and then amplified for 8–15 cycles to produce libraries for sequencing. The libraries were initially sequenced on an Illumina MiSeq sequencer and analyzed using a custom script to determine the enrichment score by calculating the ratio of signal over background at TSS over a 2-kb window. Only libraries that had the highest score above the threshold (> 5) were chosen for deeper sequencing. Two independent, biological replicates were sequenced on an Illumina NextSeq sequencer or a NovaSeq sequencer. Paired-end reads were trimmed for Illumina adaptor sequences and transposase sequences using a customized script and mapped to mm9 using Bowtie v1.1.2 (Langmead et al., 2009) with parameters -S -X2000 -m1. Duplicate reads were discarded with Samtools v0.1.18. Narrow peaks were called using MACS2 (Zhang et al., 2008) with parameters–nomodel–extsize 200–shift 104 and FDR threshold 0.05. Background removal was carried out via submitting replicates to irreproducible discovery rate (IDR) filtering (Li et al., 2011). Overlapping peaks from all samples were merged into a unique peak list, and raw read counts mapped to each peak [using bedtools multicov (Quinlan laboratory, University of Utah, Salt Lake City, UT)] for each individual sample were quantified. Differentially accessible peaks from the merged union peak list were selected with the edgeR package from Bioconductor. Cutoffs were set at log2 fold change > 1 or < −1 and P-value < 0.05. Differential peak analysis was performed in:

ASZ_NSshRNA_RasV12/Dox-4OHT versus ASZ_NSshRNA_RasV12

ASZ_TGFβRIshRNA_RasV12 versus ASZ_NSshRNA_RasV12

ASZ_TGFβRIshRNA_RasV12/Dox-4OHT versus ASZ_NSshRNA_RasV12

ASZ_cFos/Dox versus ASZ_empty/Dox

ASZ_cFos/Dox/DMSO versus ASZ_empty/Dox/DMSO

ASZ_cFos/Dox/UO126 versus ASZ_cFos/Dox/DMSO

ASZ_cFos/Dox/Afatinib versus ASZ_cFos/Dox/DMSO

Read pileups at genomic loci were imaged using Integrative Genomics Viewer (Broad Institute). Differentially accessible peaks were annotated for gene associations using the Genomic Regions Enrichment of Annotations Tool (GREAT) (Stanford University, Stanford, CA) (McLean et al., 2010), with the parameter: “single nearest gene within 50kb.”

ChIP-sequencing library preparation, sequencing and analysis

To test c-FOS chromatin occupancy, chromatin immunoprecipitation for c-FOS was carried out as described previously (Whitson et al., 2018) with minor modifications. Six days after doxycycline/4OHT induction, ASZ_RasV12_TGFβR1kd cells were crosslinked with 1% formaldehyde for 10 min and then lysed in modified RIPA buffer (50mM Tris, 150mM NaCl2, 1% Triton X-100, 0.75% SDS, 0.5% Sodium Deoxycholate), which was supplemented with protease and phosphatase inhibitor cocktail (Roche, CA). Cellular extracts were sonicated using a Covaris B208 ultrasonicator to produce chromatin fragments between 100 to 400bp. Cleared extract was incubated with 3 μg of anti-c-FOS (9F6) (CST, Beverly, MA), or non-specific IgG control antibody (CST, Beverly, MA) overnight and precipitated using protein A/G Sepharose beads. Beads were washed with ChIP wash buffer (100mM Tris pH 9.0, 500mM LiCl, 1% Igepal, 1% Deoxycholic Acid, protease and phosphatase inhibitor cocktail) and protein/DNA complexes were eluted with IP elution buffer (1%SDS, 50mM NaHCO3). Crosslinks were reversed by incubation at 67°C overnight while shaking at 1400 rpm on a thermoshaker. RNA was digested with 0.2ug/ml RNase A at 37°C for 30 min. DNA was isolated using QIAGEN MinElute PCR Purification Kit according to manufacturer’s instructions. Relative fold enrichment of c-FOS was determined by adding DNA to Brilliant II SYBR Green qPCR Master Mix Kit (Agilent Technologies). ChIP with non-specific IgG control antibody was used as a control to calculate fold enrichment. ChIP libraries from two biological replicates were sequenced using the Illumina NextSeq (400M) platform. Sequencing reads were mapped to mm9 using Bowtie2 (2.3.4.1)(Langmead et al., 2009) with parameters -p 4–very-sensitive. Duplicates are then removed using Samtools rmdup (Li et al., 2009). Peaks were identified using MACS2 (Zhang et al., 2008) with input controls and -p 0.01. Read pileups at genomic loci were imaged using Integrative Genomics Viewer (Broad Institute). Background removal was carried out using bedtools (Quinlan laboratory, University of Utah, Salt Lake City, UT) to identify overlapped peaks between replicates. High-confidence peaks were annotated for gene associations using GREAT (Stanford University, Stanford, CA) (McLean et al., 2010), with the parameter: “single nearest gene within 50kb.” Genomic annotation of ChIP peaks was carried out using the ChIPseeker R package (Yu et al., 2015).

Motif analysis of peaks from ATAC-seq and ChIP-seq

Motif analysis on peak regions was performed using HOMER function (http://homer.ucsd.edu/homer/motif/) “findMotifsGenome.pl” with default parameters to calculate the occurrence of a TF motif in peak regions compared to that in background regions. We used −log10(P-value) to rank the enrichment level of TF motifs. For high confidence c-FOS genomic binding sites, significant ATAC peaks with AP-1 motif enrichment and c-FOS ChIP-seq peaks were overlapped using bedtools (Quinlan laboratory, University of Utah, Salt Lake City, UT).

Adenovirus inoculation

AAV6-CMV-cFos-IRES-GFP (AAV-229259) was obtained from Vector Biosystems (Malvern, PA). Control AAV6-CMV-GFP (GVVC-AAV-5) was obtained from the Neurosciences Institute (Stanford, CA). 50 μL of vehicle (PBS), AAV6-CMV-GFP (control, 1×1012 particles) or AAV6-CMV-cFos-IRES-GFP (1×1012 particles) were inoculated intratumorally. Three days after inoculation, tumors were harvested and processed as previously described for HE and immunofluorescence staining or mechanically and enzymatically digested for flow cytometry cell sorting and mRNA extraction.

Flow cytometry and FACS

For flow cytometry analyses, tumors were excised, mechanically disrupted, enzymatically digested in 0.5% collagenase for up to 1 h, and then strained through a 70 μm filter to obtain single-cell suspensions. Cells were washed twice using FACS buffer (2% BSA/PBS) before sorting. Samples were acquired with a FACS Aria II (Becton Dickinson, Franklin Lakes, NJ) and data analyzed using FlowJo (Ashland, OR). For mRNA extraction, tumor cells were sorted into RFPposGFPpos and RFPposGFPneg populations. The dead cell proportion (< 10%) and purity of the sorted samples were assessed using trypan blue and fluorescence microscopy respectively. A minimum of 104 viable sorted cells was used for mRNA extraction. Independent biological replicates from two tumors were used per condition. All samples were FACS analyzed with the same parameters.

Human primary tumor explant culture and drug treatment

Freshly resected tumors were obtained from patients undergoing dermatological surgery. The tumor subtype and nuclear c-FOS expression level were verified through immediate histological examination and immunostaining. Specimens were minced and cultured in EpiLife medium (Life Technologies) supplemented with 0.05 mM CaCl2 with DMSO or 5 μM of afatinib for 24 h. Drug-treated tissues were homogenized using a tissue homogenizer (Thermo Fisher Scientific, Waltham, MA) before RNA extraction as previously described. RNA extracts were used to carry out qRT-PCR with TaqMan probes for human GLI1, CD44, KRT14 (Thermo Fisher Scientific, Waltham, MA). KRT14 was used to normalize gene mRNA expression to the keratinocytic content of tissue samples.

Co-immunoprecipitation

ASZ cells with dox-inducible Myc-tagged c-FOS construct were treated with doxycycline in serum-free medium for 24h prior to lysis. Cells were lysed in Mag c-Myc IP/Co-IP Buffer-1 from the Pierce Magnetic c-Myc-Tag IP/Co-IP Kit (Pierce Biotechnology, Rockford, IL, USA), supplemented with protease and phosphatase inhibitors, on ice for 30 minutes. Lysates were cleared (10 min, 13,000 g) and 10% of each supernatant was set aside as input. The rest of the supernatants were incubated with pre-washed Pierce Anti-c-Myc Magnetic Beads (Pierce Biotechnology, Rockford, IL, USA) overnight at 4°C. Beads were then collected with a magnetic stand, washed and protein samples eluted following manufacturer’s protocol.

In vitro EGFR (afatinib) and MAPK (UO126) inhibition assay

ASZ cells with dox-inducible c-FOS construct were treated with doxycycline in serum-free medium for 24h. Cells were then treated with either 5μM of afatinib, 10μM of UO126 or DMSO for 6h. ASZ cells with an empty construct were treated with doxycycline in serum-free medium for 24h as the control. Following the inhibition assay, cells were harvested for subsequent analysis using RNA-seq, ATAC-seq and CUT&RUN-qPCR.

CUT&RUN-qPCR

CUT&RUN was performed with the CUT&RUN Assay Kit (86652, CST) according to manufacturer’s protocol. Briefly, cells were harvested, washed, and bound to activated Concanavalin A-coated magnetic beads and permeabilized. The bead-cell complex was incubated overnight with either an anti-c-Fos antibody (CST) or an IgG antibody control at 4°C. Cells were then washed with digitonin buffer, resuspended in 50 μl of pAG/Mnase and incubated for 2h at 4°C. Enriched DNA fragments were purified using QIAGEN MinElute PCR Purification Kit.

Quantitative PCR was performed with SimpleChIP® Universal qPCR Master Mix (#88989, CST). Each reaction was performed in replicate and Ct values were normalized to IgG controls. Primer sequences for qPCR analysis are as follows:

Fosl1 forward primer: 5′-AAGTCGGTCGCTTTCTGTCTGTA-3′

Fosl1 reverse primer: 5′-GAACTTCACGACCCTCTGCTC-3′

Epgn forward primer: 5′-CTTGCATCCTCCAAAGCTACCG-3′

Epgn reverse primer: 5′-AGTTGGCAGATTTAAAGGCTCCTA-3′

QUANTIFICATION AND STATISTICAL ANALYSIS

In general, data represent results from three or more independent biological samples, unless otherwise described. Deep sequencing data are from two biological samples per condition. Immunofluorescence analyses of human tumors were performed on n ≥ 3 tumors per condition. Immunofluorescence analyses of cultured cells were performed on n ≥ 40 cells per condition. Scale bar annotation and quantification of pixel intensity was performed using ImageJ software (NIH; Bethesda, Maryland) (Schindelin et al., 2012). Bar and line graph results reflect the mean with standard deviation (SD). Statistical comparisons were performed using unpaired two-sided Student’s t test or one-way ANOVA with Tukey’s post-test for multiple comparisons. Statistical comparisons of proportions were performed using z-test. The software used for statistical analysis is GraphPad Prism (La Jolla, CA), with the annotations: ns, non-significant, *p < 0.05, **p < 0.01, and ***p < 0.001. Principal component analysis (PCA) was used for to compare transcriptomic profiles of normal skin, BCCs, BSCs, well-differentiated and poorly-differentiated SCCs. Central distribution and distances between tumor types were calculated from PCA. Details of statistical methods for specific analysis are described in corresponding methods and figure legends. A normal distribution was observed for all data.