LY6D marks pre-existing resistant basosquamous tumor subpopulations

Targeted therapies have proven beneficial in treating cancers, but dynamic tumor heterogeneity and an array of pre-existing and induced resistance mechanisms act to circumvent treatment. Growing evidence points to inherent cellular plasticity or phenotype switching as one of the most challenging resistance mechanisms. While epithelial stem cell plasticity confers a remarkable ability of cells to toggle their behaviors and cellular identities, drivers of changing phenotypes and markers for populations that correlate with poor clinical outcomes are currently lacking1.

Non-melanoma skin cancers or keratocarcinomas such as basal cell carcinomas (BCCs) and squamous cell carcinomas (SCCs) represent ideal model systems to study dynamic tumor heterogeneity and plasticity because of their accessibility for serial biopsies and the well-defined lineage relationships within the skin2,3. BCC growth depends on the canonical activation of the Hedgehog (HH) signaling pathway, primarily driven through various mutations in genes such as Patched (PTCH1) or the G-protein coupled receptor Smoothened (SMO)4. While HH pathway inhibitors targeting SMO (SMOi) are the primary therapy in advanced non-surgically resectable cases, 60% are already resistant2. Furthermore, acquired resistance occurs in ~20% of patients who initially respond to SMOi, reinforcing the need to delineate tumor evolution mechanisms2.

Although distinct in their molecular drivers, clinical observations point to a fluid relationship between BCCs and SCCs, which have clear implications for treatment strategies. First, it is routinely observed that BCCs treated with SMOi can develop resistance in the form of a transition towards SCCs after being re-biopsied5,6. In contrast to BCCs, which exhibit features of the hair follicle (HF) epithelium progenitors, SCCs are not HH-driven and resemble the interfollicular epidermis (IFE)2,3. In line with this clinical observation, some BCC mouse models point to a SMOi-dependent lineage transition resistance mechanism, producing a persister population of cells characterized by loss of HH dependence and toggling towards a differentiated IFE-like cell fate7. Second, the existence of sporadic basosquamous carcinomas (BSCs), which show features of both BCCs and SCCs, further points to the presence of a basosquamous transition (BST) in keratocarcinomas8,9. Further investigation found that the SCC portions of BSCs are BCC-derived but lose their dependence on HH signaling9. With both sporadic and induced observations, BST etiology remains complex, with various drivers and inducers such as SMOi, transcriptional inducers like c-FOS, or subsequent mutations such as ARID1A5,6,9,10.

A critical unresolved question is whether populations that undergo BST originate from the selection of small pre-existing resistant populations or are induced in response to the therapy. It is unclear what the relationship is between these observed populations, the mechanisms that drive their accumulation, the key markers that define and track them, and whether they exist in other epithelial tumors.

Here we deconvolute BST mechanisms utilizing the identified stem cell surface marker Lymphocyte Antigen 6 Family Member D (LY6D), which marks and tracks basosquamous populations. LY6D marks the SMOi-induced IFE-like resistant population, but importantly, it also marks a pre-existing population in naive human sporadic BCCs. We show that LY6D+ tumor cells lie on a differentiation spectrum between BCC and SCCs, and resemble basosquamous cell carcinoma (BSC), highlighting pathway switching as a cellular plasticity-mediated resistance mechanism. Further, we show how SMOi expedites the transition towards the LY6D+ state and demonstrate that the LY6D+ population marker can identify other basosquamous states in epithelial tumors with poor outcomes.

LY6D marks a pre-existing spatially localized basosquamous BCC population

Our approach to deconvolute BST-associated plasticity aimed to define populations and markers within BCCs, identify potential signs of BST, and then interrogate factors that promote interconversion. We generated a comprehensive dataset by merging seven existing and new (BCC19, BCC20, and BCC21) datasets from naive human BCC tumors with varying histological features and then used appropriate canonical correlation analysis (CCA) to mitigate patient-to-patient variability (Fig. 1a, b, Supplementary Fig. 1a–c, e)11. CopyKAT, which specifically differentiates tumor and normal tissue via copy number variations (CNVs), validated our ability to discriminate between tumor and normal tissue (Supplementary Fig. 1d)12. As a result, we found seven distinct clusters, each with its unique marker genes and corresponding Gene Ontology (GO) associated terms (Fig. 1b, Supplementary Fig. 1e–h, Supplementary Data 1).

Surprisingly, the naive BCC dataset already contained SCC-like resistant populations arising from BST. We used BCC and SCC gene signatures to score the different cellular populations and found a spectrum of BCC- and SCC-like states among the different clusters where C6 from our merged dataset (Merged-Hu-C6) had significant enrichment for SCC-associated gene signatures, identifying the cluster of SCC-like tumors’ pre-existence in naive BCCs (Fig. 1b, c)10. Consistent with the SCC-like state, Merged-Hu-C6 was enriched for epidermis-related GO terms, keratinocyte differentiation gene scoring, and differentiated related cell types such as spinous and granular IFE but not the HH responsive HF epithelium (Fig. 1d, e, Supplementary Fig. 1h–i)13. This observation suggests that merged-Hu-C6 represents an intermediate state between BCC and SCC, resembling a BSC/basosquamous state, which may be resistant to SMOi9. To assess the incidence of this population in naive BCCs, we interrogated the presence of the Merged-Hu-C6 population in all seven naive BCCs. We found the basosquamous population present in all, albeit with varied frequency, supporting the commonality of this pre-existing population (Fig. 1f). Furthermore, we noted that other clusters representing intermediate SCC-like clusters with an enriched for (Cluster 2) and non-SCC-like (Cluster 3) were also distributed in all seven naive BCCs, suggesting a spectrum of tumor cell states (Fig. 1f).

To prospectively isolate the Merged-Hu-C6 basosquamous population, we used highly expressed marker genes specific to Merged-Hu-C6 with subsequent verification by COMET to probe our scRNA-Seq for associated surface markers and identified LY6D, CLDN4, LYPD3, and TACSTD2 (Fig. 1g, Supplementary Fig. 1j)14. We have found the surface markers interchangeable, especially with subsequent staining experiments, but found that LY6D had the greatest Merged-Hu-C6 specificity (Fig. 1g). LY6D is part of the Ly6 family of lymphocyte lineage and stem cell markers that have roles in marking solid organ lineages as well15. Notably, LY6D expression is also enriched in the well-differentiated human SCC cluster and marks the overlap between human BCC and SCC scRNA-Seq data and the IFE suprabasal cells of human foreskin scRNA-Seq data, further drawing a parallel to the basosquamous cellular state observed in naive human BCCs (Supplementary Fig. 2a–g)16,17. We then examined naive murine BCCs from the Ptch1+/−;K14-creER;p53fl/fl strain to determine whether the same marker could be used in mouse tumor biology (Fig. 1a, h). Highlighting its utility, Ly6d also marked the basosquamous population in naive murine BCC tumors (Fig. 1j). Although multiple clusters expressed Ly6d, perhaps suggesting increased Ly6d correlates with larger and more advanced tumors (a point we will return to), only C5 (Mouse-C5) expressed all basosquamous markers (Ly6d, Cldn4, Lypd3, and Tacstd2) and was the most enriched for IFE/keratinocyte differentiation-related genes (Fig. 1h–j, Supplementary Fig. 2h, i, Supplementary Data 2). Additionally, the intersection between the human Merged-Hu-C6 and mouse Mouse-C5 clusters reveals a 77 gene set that includes LY6D, suggesting a similar conserved basosquamous state across species (Fig. 1k; Supplementary Data 3). We conclude that LY6D/Ly6d marks the basosquamous population in human and mouse tumors.

Using our surface markers as basosquamous readouts, we found that the populations occupy unique spatial localizations within the tumor. RNAScope delineated the population, with LY6D marking the cells in the most central portions of tumor nodules (Fig. 1l, Supplementary Fig. 2j). Like the suprabasal cells in the normal human epidermis, these clusters of cells were noticeably separated from the basement membrane zone (BMZ) and did not interface with the stroma (Fig. 1l, Supplementary Fig. 2j). Mouse BCCs had similar patterns with Ly6d expression found centrally localized in nodules (Fig. 1m). Extensive expression was also observed in well-differentiated human SCCs and the more SCC-like portions of BSCs compared to the BCC-like regions (Fig. 1m).

We next interrogated the status of the HH signaling pathway in LY6D+ basosquamous cells by checking whether LY6D+ cells overlapped with the pathway target gene GLI1. In line with the previous finding that Merged-Hu-C6 was not enriched for HH-related genes, LY6D and the associated surface markers failed to colocalize with GLI1 (Fig. 1n, Supplementary Figs. 1i, 2k). Additionally, scATAC-Seq (which we will discuss at length later) lineage trajectories suggested an inverse relationship between GLI1 and LY6D (Supplementary Fig. 2l). Finally, we confirmed by RNAScope that LY6D+ clusters in both naive BCCs and BSCs do not colocalize with GLI1 (Fig. 1o).

Overall, we conclude that LY6D marks a previously unidentified spatially localized basosquamous-like population lacking GLI signaling, which is present before SMOi treatment. This population lacks GLI1 expression, is spatially localized, and is demarcated by TACSTD2, LYPD3, and most specifically by LY6D (Fig. 1p). Our analysis of multiple naive BCCs, BSCs, and SCCs suggests levels of LY6D can vary from tumor to tumor but dramatically increase along the BST spectrum towards SCC (Fig. 1q).

BCCs transition towards an LY6D+ state while maintaining proliferative capacity

The spatial localization of LY6D− and LY6D+ tumor populations were of great interest, prompting further investigation into their potential relationship. We utilized several computational approaches to support the lineage progression from LY6D− to LY6D+ tumor states. First, principal component analysis (PCA) of the merged human BCC scRNA-seq dataset indicates the LY6D+ Merged-Hu-C6 is a BCC lineage endpoint (Fig. 2a). Furthermore, PCA and Monocle analysis of BCC21 (where BCC21 C5 expresses high levels of LY6D, referred to as BCC21-Hu-C5) reinforce a transition from an LY6D− state to an LY6D+ state (Fig. 2b, c, Supplementary Fig. 2l, m, Supplementary Data 4). These observations were further validated with Monocle pseudotime analysis of our mouse scRNA-Seq data, revealing the Ly6d− to Ly6d+ transition (Supplementary Fig. 2n).

Second, single-cell ATAC-Seq (scATAC-Seq) supports sequential chromatin accessibility changes between LY6D− and LY6D+ populations. We analyzed three human BCC tumors with accompanying scRNA-Seq data and used the ArchR algorithm to define cluster relationships (Fig. 2d, Supplementary Fig. 3a–d, Supplementary Data 5)18. An example analysis of the tumor epithelial subset of BCC21 found at least four distinct clusters, which varied in their chromatin accessibilities near various marker genes (Fig. 2e, Supplementary Fig 3e–k, Supplementary Data 6–8). Integration of BCC21 scRNA-Seq data identified scATAC-Seq C1 as highly similar to BCC21-Hu-C5 (Fig. 2e, Supplementary Fig. 3l). scATAC-Seq C1 showed enhanced accessibility for LY6D along with reduction for GLI1, and like the lineage predictions from the scRNA-Seq, the scATAC-Seq lineage prediction tool via ArchR predicts scATAC-Seq C1 as the endpoint (Fig. 2f, g). Looking at the scATAC-Seq tracks, we note distinct peaks associated with scATAC-Seq C1 in the promoter regions of LY6D (Fig. 2h). Lastly, consistent with sequential cluster transition, we note distinct transcription factor motifs appear and close (Fig. 2i, Supplementary Fig. 3m, n). Previous work has implicated the AP-1 family motifs in the transition from BCCs to SCCs10. We found that AP-1 family motifs were enriched at the early stages of the transition (scATAC-Seq C3 to C4 and C4 to C2), but that enriched motifs such as RUNX1 and various forkhead box family of transcription factors were involved in the regulation of the LY6D+ state (Fig. 2i). Our chromatin dynamics analysis strongly supports a spontaneous transition to an LY6D+ population.

Third, to functionally confirm that the transition from an LY6D− to LY6D+ state is tumor cell-autonomous and spontaneous, we generated an SMOi-sensitive mouse BCC organoid system to culture bulk or sorted populations of tumor epithelial cells based on specific surface markers (Fig. 2j, Supplementary Fig. 5a–c). We sorted Ly6d− cells, cultured them as organoids for 14 days, and found that after 14 days in culture, expression levels of Ly6d in the organoids were significantly higher compared to levels at initial plating as well as comparable to organoids from Ly6d+ cells grown for 14 days (Fig. 2j, k, m Supplementary Fig 4d). Additionally, when we cultured Ly6d− organoids, we detected Ly6d+ cells by flow analysis at D14, suggesting the ability of Ly6d− to transition Ly6d− to the Ly6d+ state (Fig. 2n). Overall, this data reinforces the cell’s autonomous nature of the Ly6d− to Ly6d+ transition. Additionally, we noted a significant growth advantage of the Ly6d+ organoids compared to their Ly6d− counterparts (Fig. 2k, l). Importantly, this growth advantage was specific to tumor Ly6d+ cells as Ly6d+ cells from the normal epidermis were growth-arrested compared to normal Ly6d− cells (Fig. 2k, l, Supplementary Fig. 5g). We conclude that despite the spontaneous suprabasal conversion from LY6D− to LY6D+ tumor cells, LY6D+ BCC cells do not exit the cell cycle and maintain their proliferative potential.

Distinct epigenetic regulation and biological processes associated with tumor epithelial cells

Because LY6D+ cells in the suprabasal layers of the normal IFE normally uniformly lead to proliferative arrest and terminal differentiation, we investigated the difference between the interconversion from LY6D− to LY6D+ in both normal and BCC contexts. We took advantage of the ability to discriminate between normal and tumor cells in our human scATAC-Seq and scRNA-Seq data via chromatin accessibility and expression of GLI1, respectively (Fig. 3a–c, Supplementary Data 9). Using four scRNA-Seq clusters defined by both normal and BCC cells as well as LY6D− and LY6D+ fractions (Normal: LY6DLow, Normal: LY6DHi, Tumor: LY6DLow, and Tumor: LY6DHi), we were able to impute and integrate these clusters onto the scATAC-Seq clustering for subsequent analysis (Fig. 3d). We would presume that the basis of these key cellular differences stem from the large difference in mutations between the tumor and normal epithelium. We had previously used CopyKat to differentiate tumor and normal epithelial tissue through CNV inferences (Supplementary Fig. 1d). Visualization of the altered CNVs between Normal: LY6DHi and Tumor: LY6DHi indicate extensive alterations in CNVs, likely reflecting the altered mutational status (Supplementary Fig. 3o)12.

a UMAP plot of scATAC-Seq for both the tumor epithelia and normal epithelia regions within BCC21. b Feature plot of the GeneScoreMatrix for GLI1 for human BCC21 scATAC-Seq data as shown in Fig. 3a. c The differential marker features from the scATAC-Seq data of human BCC21 between tumor epithelia and normal epithelia populations. d UMAP plots of scRNA-Seq (left) and scATAC-Seq (right) data integration for the tumor epithelia and normal epithelia populations of human BCC21. e GO analysis for marker genes from (1) scRNA-Seq (top) and (2) scATAC-Seq (bottom) of C1 vs. C2. f Edu% of Ly6d+ cells from either normal mouse skin organoids or tumor organoids. Each point represents Edu% from distinct organoids (n = 4) derived from BCCs or normal mouse skin from different mice. g Heatmap indicating the enriched motifs in the various clusters from normal and BCC cells. h Feature plot of the MotifMatrix for various motifs identified in Fig. 3g superimposed onto the different clusters from the normal and BCC cells to indicate the specificity of enrichment. i Feature plot highlighting expression of the various genes (JUNB, SMARCC1, FOXN3, and RUNX1), whose motifs were identified in h. j Brightfield images of organoids from mouse BCC tumors that were treated with SMOi, Runx1i, or SMOi with Runx1i. Inhibitor treatments were at 1 μM. k Quantification of the organoid area after growth with the various drugs from distinct primary tumors (n = 4). l Diagram of the CRISPR strategy to target Foxn3 and Runx1 in ASZs. m Gene expression of Gli1, Tacstd2, and Lypd3 in wildtype and different clones of Foxn3-KO and Runx1-KO cells. Each point represents technical replicates (n = 4) of gene expression values of the specified gene normalized to Gapdh. Length of each scale bar is noted in figure. f, k, and m error bars represent mean ±SD. f, k p values were calculated using unpaired, two-tailed t test. Source data are provided as a Source Data file.

Full size image

To further understand the underlying biological processes enriched in the LY6D+ BCC cells compared to LY6D+ normal cells, we did GO analysis of differential marker genes from scRNA-Seq and peaks (in proximity to genes) from scATAC-Seq between LY6D+ (C1) and LY6D− (C2). We identified similar terms associated with differentiation, proliferation, locomotion, and adhesion (Fig. 3e). Consistent with the observed increased proliferative differences, we noted prominent Mki67+/Ly6d+ cells in our mouse scRNA-Seq data and subsequently found Edu incorporation of Ly6d+ cells from the tumor are more significant than from the normal epidermis (Fig. 3f, Supplementary Fig. 5e, f).

Motif analysis of these four clusters allowed found significant enrichment for forkhead box family of transcription factors (such as FOXN3) within the BCC compared to normal cells (Fig. 3g). Interestingly, we found that although enriched in the LY6DHi BCC cells compared to the LY6DLow BCC cells, the AP-1 motifs (such as JUNB) and SMARCC1 were largely specific to normal cells compared to BCC cells (Fig. 3g, h). Closer examination identified that the RUNX1 motif was uniquely enriched in the LY6DHi portion of the BCC, suggesting that this tumor-specific motif might regulate the basosquamous state transition (Fig. 3g, h). Although there was motif enrichment specificity, the expression of the identified transcription factors was largely consistent throughout the tumor and normal tissue (Fig. 3i). To functionally test whether FOXN3 and RUNX1 promoted or inhibited distinct tumor states, we performed two separate experiments. First, we tested the effects of the Runx1 molecular inhibitor (Runx1i) with and without Smoi in our mouse tumor organoid model (Fig. 3j, Supplementary Fig. 4a). We found that both Smoi and Runx1i had significant effects on organoid growth, and together, the inhibitors did not display synergism, suggesting that Runx1 acts positively on BCC growth through the same downstream mechanism (Fig. 3j, k, Supplementary Fig. 4a). Second, we generated Foxn3 (Foxn3-KO) and Runx1 (Runx1-KO) mutant murine BCC cell lines (Fig. 3l, Supplementary Fig. 4b). Compared to wildtype controls, both Foxn3-KO and Runx1-KO cells had reduced levels of Gli1 (Fig. 3m). Furthermore, we noted that along with the reduction in Gli1, we saw an increase in the expression levels of BST-associated surface markers Lypd3 and Tacstd2 (Fig. 3m). Together, these data support a role for Runx1 and Foxn3 in maintaining the LY6D− BCC state through repression of the BST transition.

Physical features and genetics influence LY6D+ cellular accumulation

While we have shown that LY6D+ BCCs lack HH/GLI signaling, form spontaneously without treatment, and remain proliferative compared to normal tissue, we further investigated the determinants affecting LY6D+ tumor cell accumulation. As previously mentioned, we noted that the bulk expression of Ly6d was higher in our mouse model compared to our human tumors. Upon quantitative analysis, we noted higher Ly6d bulk expression in larger mouse tumors when compared to paired smaller tumors (Fig. 4a). These results would suggest that the amount of Ly6d increases as tumors progress and get larger in size. Additionally, within BCC21, the expression of LY6D correlated with large nodules (Fig. 4b). Furthermore, basosquamous regions were typically found 2–3 cells away from the BMZ, suggesting that overall tumor size, intratumor nodular size, and distance from the BMZ contribute to LY6D/Ly6d levels (Fig. 4c).

a RNAScope and quantification of Ly6d+ (yellow) clone size from small or large tumors from the same mouse. Each point represents the average clone size from small and large tumors from the same primary mouse (n = 4). Lines between dots indicate the small and large tumors from the same primary tumor. b Tumor nodule areas from human BCC21 that are either LY6D+ or LY6D− from Fig. 1l. The median area is indicated. c The number of cells between basosquamous regions and the BMZ from human BCC21 in Fig. 1l. d H/E and corresponding RNAScope of the spatial localization of LY6D (yellow) and GLI1 (red) in resistant human BCCs, which is co-stained for K14 protein (green). The upper tumor displays an example of an rBCC with basosquamous/SCC-like features, while the lower tumor does not. Zoomed-in regions are indicated with yellow dashed boxes. Regions enriched for LY6D are indicated with a yellow arrow. Regions enriched for GLI1 are indicated with a red arrow. n = 2. e Feature plot of two resistant human BCCs showing dual expression of GLI1 and LY6D. R8 (left) has basosquamous markers, where R6 (right) does not. f Persister score (from Biehns et al., 2018) for the various BCC datasets with the Merged-Hu-C6 (left), BCC21-Hu-C5 (center), and Mouse-C5 (right) all indicated. g Violin plot for LY6D expression (top), Persister scoring (middle), and mutational status (bottom) for various resistant human BCCs (from Yost et al., 2019). h H/E examples of baseline (left) and SMOi (right) treated tumors from Gorlin patients. Below is RNAScope indicating the spatial localization of LY6D (yellow) and GLI1 (red), which is co-stained with K14 protein (green) in baseline (left) and on SMOi (right) Gorlin patients. The SMOi (right) treated RNAScope image is matched with the H/E above. Zoomed-in regions are indicated with yellow dashed boxes. Yellow arrows point to regions enriched for LY6D. n = 1. Length of each scale bar is noted in figure. a p values were calculated using unpaired, two-tailed t test. Source data are provided as a Source Data file.

Full size image

Along with tumor size, drug therapy also affects the proportion of LY6D+ cells. When resistant BCCs (rBCCs) are re-biopsied after SMOi, many present extensive SCC-like phenotypes (Fig. 4d)5. Examining rBCCs with SCC-like features by RNAScope reveals LY6D expression that does not colocalize with GLI1 (Fig. 4d). Conversely, rBCCs without SCC-like features lack LY6D expression but maintain GLI1 levels throughout the tumor (Fig. 4d). Similar observations can be seen in rBCC scRNA-Seq data, with some tumors homogeneously expressing LY6D or GLI1 (Fig. 4e, Supplementary Fig. 5h)19. Overall, this confirms that within naive BCC tumors, targeted therapy via SMOi selects for LY6D+ populations.

Previous work has identified a persister population of cells that maintain resistance to SMOi, which have IFE-like characteristics20. We found that our LY6D+/Ly6d+ cells shared considerable overlap with this population of cells in our human and mouse naive datasets (Fig. 4f). Furthermore, rBCCs with high levels of LY6D subsequently had a high persister score (Fig. 4g)19. In contrast to previous observations, this data would suggest that the persister population is spontaneously generated in naive sporadic tumors, which offers a mode of BCC resistance to SMOi.

Not all rBCCs express LY6D, so we postulated that variation could be based on differential mutational landscapes that yield differential resistance. Consistent with this idea, we noted a lack of LY6D expression in those tumors with SMO activating mutations, irrespective of size (Fig. 4g). Additionally, we noted correlations of LY6D expression in rBCCs with NOTCH1 and NOTCH2 mutations, typically seen in SCCs (Fig. 4g)9. Patients with basal cell nevus syndrome, also known as Gorlin syndrome, have a rare heritable inactivation of PTCH1, predisposing them to BCCs with relatively few background mutations21. Clinical trial data demonstrate that while SMOi therapy shrinks syndromic BCCs, they typically reform in the same location after SMOi is discontinued from the persister cell population22. Examining naive and SMOi-treated tumors, we find that naive syndromic tumors express little LY6D but high levels of GLI1 (Fig. 4h). By contrast, persister cell clusters in SMOi-treated tumors enrich for LY6D and lack GLI1 expression, reinforcing the therapy dependent LY6D basosquamous state in these tumors (Fig. 4h). Altogether these observations suggest that tumors without an SMO activating mutation employ both spontaneous and drug-induced conversion strategies towards an LY6D+ basosquamous persister state.

Therapeutic intervention can alter SMOi Resistant LY6D+ cellular accumulation

To better understand how SMOi and other drug interventions regulate the transition towards an LY6D+ state, we utilized our organoid system, plating bulk tumors and subjecting them to various inhibitors to examine how they would influence the levels of Ly6d (Fig. 5a). To test SMOi selectively, we treated organoids with SMOi and then checked the apoptotic levels in both Ly6d− and Ly6d+ cells by flow cytometry (Fig. 5b, Supplementary Fig. 5i). We found enhanced cell death in the Ly6d− cells suggesting that, as one might predict, SMOi targeted the cells enhanced with HH signaling (Fig. 5c). Treatment with the MRTF inhibitor CCG, which has been shown to associate with the LY6D+ nMRTF state switching, did not show similar acute selectivity (Fig. 5c)10,11.

a Schematic diagram for tumor dissociation and subsequent experimental designs with in vitro organoid or in vivo allograft drug treatments. b Schematic diagram of organoid treatments designed to examine acute and long-term drug treatments. c Fold change in Apotracker Green in organoids treated with either SMOi or CCG. Each point represents Apotracker Green% of Ly6d− or Ly6d+ cells from distinct drug treated organoids derived from different primary tumors (n = 3). d Fold change in Ly6d% in organoids treated with either SMOi or CCG. Each point represents the fold change in Ly6d% from distinct drug treated organoids derived from different primary tumors (n = 3). e UMAP plot of merged scRNA-Seq data from control and SMOi-treated mouse tumors by clusters (left) and treatment groups (right). f Feature plot for Ly6d expression data for Fig. 5c. g The proportion of cells from control or SMOi treated tumors with a Ly6d expression level >2. h Keratinocyte Differentiation Score for the control or SMOi treated tumors in Fig. 5c. i Mouse-C5 marker gene score from the mouse scRNA-Seq data for the control or SMOi treated tumors in Fig. 5c. j Persister score for the control or SMOi treated tumors in Fig. 5c. c, d Error bars represent mean +/−SD. c, d p values were calculated using unpaired, two-tailed t test. Source data are provided as a Source Data file.

Full size image

To test the long-term effect of drug intervention, we treated organoids for 14 days immediately after plating, and to test the acute effect, we treated already-formed organoids for 48 h (Fig. 5b). As predicted from our previous observations, treatment with SMOi led to a modest increase in Ly6d, suggesting that long-term treatment would select for a Ly6d+ state (Fig. 5d; Supplementary Fig. 5j). In contrast, when organoids are treated with CCG, we see a dramatic decrease in Ly6d levels (Fig. 5d; Supplementary Fig. 5k).

To examine the global long-term in vivo effects of SMOi on BCCs, we transplanted the same primary tumor and then treated the allograft tumors with either DMSO or SMOi (Fig. 5a). scRNA-Seq of the treated tumors revealed distinct clustering, with enhanced Ly6d levels in the SMOi treated tumor (Fig. 5e–f, Supplementary Data 10). In addition, the global analysis revealed enhanced differentiation and enrichment for Mu-C5 marker genes from the primary mouse scRNA-Seq (Fig. 5h, i). Finally, we noted an enrichment for the resistant persister cell state in the SMOi-treated tumor (Fig. 5j). This data supports the notion that targeted therapy through SMOi can induce the LY6D+ persister basosquamous state.

LY6D+ populations are targetable and exist in other epithelial tumors, correlating with poor clinical outcomes

Our studies indicate that LY6D distinguishes between mouse and human skin BCC- and SCC-like states, raising the possibility that it serves as a general diagnosis and treatment biomarker for other basosquamous cancers. Examining LY6D levels across all major cancers, we noted high levels in adenoid cystic carcinoma (ACC), bladder urothelial carcinoma (BLCA), cervical squamous cell carcinoma, and endocervical adenocarcinoma (CESC), esophageal carcinoma (ESCA), head and neck squamous cell carcinoma (HNSC), lung squamous cell carcinoma (LUSC), and pancreatic adenocarcinoma (PAAD) (Fig. 6a). Compared to normal tissue, LY6D levels were higher in the tumor tissue of CESC, LUSC, and PAAD (Fig. 6b, Supplementary Fig. 6a). Poor clinical outcomes were associated with high levels of LY6D in CESC and LUSC but dramatically higher in PAAD (Fig. 6c, Supplementary Fig. 6b).

To understand the range of LY6D expressed in different pancreatic tumors, we examined the ductal cells from a sizeable scRNA-Seq dataset containing both normal and tumor samples, noting a spectrum of LY6D expression (Fig. 6d, Supplementary Fig. 6c, d)23. Furthermore, we noted that tumor ductal cells with higher levels of LY6D had a higher persister score, raising the possibility of a similar squamous-associated state observed in BCCs (Fig. 6e).

Squamous features are sometimes noted in PAAD, referred to as adenosquamous carcinoma of the pancreas (ASCP). Reported as a generally rare pancreatic cancer event (with frequencies ranging between 0.3–10% of cases), ASCP typically has dramatically poor clinical outcomes relative to PAAD24,25. The frequency range stems from the difficulties in histologically identifying ASCP, indicating a clear unmet clinical need to have markers that could distinguish PAAD and ASCP. To determine whether LY6D is a readout for ASCP, we examined histologically defined patient tumors from conventional PAAD as well as ASCP (Fig. 6f). We found no or very low levels of LY6D in the K19+ ducts of PAAD samples or the inflamed pancreatitis state (Fig. 6f, g). However, we did find high LY6D expression uniformly associated with the ASCP state (Fig. 6f, g). These observations suggest that LY6D defines this notoriously challenging to identify ASCP tumor state.

To determine whether LY6D has a functional role beyond its clear utility as a biomarker of squamous states, we utilized our tumor organoid model and introduced anti-Ly6d antibodies (Fig. 6h). We found that the addition of anti-Ly6d antibodies led to reduced organoid size, suggesting that targeting Ly6d may, in fact, have therapeutic significance (Fig. 6h). Overall, this work highlights that LY6D is expressed in multiple cancers and can be targeted, suggesting that LY6D has functional relevance beyond a biomarker.

Using BCCs as a model system, we have identified a spontaneously occurring resistant basosquamous state marked by LY6D surface expression and show that accumulation can be enhanced or reduced by tumor mutations or drug treatment strategies. Highlighting its utility, we show LY6D expression in other cancers, such as the ASCP subtype of pancreatic cancer, which correlates with extremely poor clinical outcomes.

While BCC tumors can develop resistance through cell identity switching between BCC and SCC, whether the transition occurs spontaneously or after therapy remained unclear2,9,20,26. Identifying the basosquamous state marked by LY6D and lack of GLI1 allowed us to track and characterize this resistant population, demonstrating its appearance even before SMOi therapy, which may explain the high failure of monotherapy. Moreover, we unify previous observations and show the LY6D marks persister cells that previously were found only after SMOi treatment20. We identify this population before SMOi but show that SMOi can further promote LY6D+ tumor accumulation. Consistent with the low SMO activity driving LY6D+ accumulation, human rBCCs containing high levels of SMO activity through activating mutations lacked spontaneous LY6D+ tumor accumulation. By contrast, Gorlin’s syndrome patients have a low base level of spontaneous LY6D+ cells, respond well to SMOi, but induce residual LY6D+ persister tumors that recur upon SMOi cessation21.

It is interesting to draw parallels between the LY6D+ cell in BCCs and whether they can potentially serve as SCC cancer stem cells (CSCs), giving rise to the SCC-associated tumor epithelium observed in multiple contexts such as rBCCs and BSCs27. Data from this work and others would suggest that persister populations from Gorlin patients as well as mouse models of BCC tumor persistence enter an LY6D+/IFE-like state with SMOi but then return to GLI1+/HH-dependent state after stopping therapy and not an SCC-like state20. This evidence would suggest that additional factors may drive cells from a BCC-associated LY6D+ state to an SCC-associated LY6D+ state with the potential to grow out an SCC-like tumor epithelium. Evidence from BSCs would suggest that additional sporadic mutations such as ARID1A are needed to make the transition to an LY6D+ state capable of giving rise to an SCC-associated tumor epithelium9. In the case of rBCCs, it has been made clear from this work that SMOi can lead to an accumulation of LY6D+ cells suggesting that drug intervention may drive LY6D+ cells into a similar state capable of giving rise to an SCC-associated tumor epithelium. These data demonstrate the key role of SMO activity and background genetic mutations in the accumulation kinetics of LY6D+ resistant tumors and suggest kinetics may be sensitive to other therapeutic, genetic, or environmental determinants.

Our work also provides additional insights into signaling pathways that may promote LY6D+ BCC tumors. The central localization within tumor nodules of the basosquamous state and distance from the BMZ provides insight and speculation as to what factors may drive this intrinsic resistance. As the tumor and tumor nodules grow, the distance from the characteristic BCC tumor stroma becomes larger, suggesting that lack of stroma (its contact and/or secreted factors) may drive towards an LY6D+ state. Deconvoluting the contribution of the BCC stroma to general tumor growth and potential resistance has become of clear clinical importance as immunotherapy approaches are now being used to treat advanced BCCs19. In addition, our studies rely on previous work to identify signaling pathways that can also influence LY6D+ tumor accumulation10,11. LY6D reduction through blocking transition pathways therapeutics targeting MRTF or EGFR reduce LY6D accumulation, raising the possibility of using these drugs in concert with SMOi therapy to minimize the accumulation of LY6D-associated resistance in both syndromic and sporadic BCCs. Furthermore, our work highlights the therapeutic implications of targeting LY6D itself may have functional roles, but additional considerations are needed to target tumor-specific LY6D+ cells.

The existence of LY6D surface expression in other epithelial cancers suggests that the work on BCCs will apply to more lethal cancers such as ASCP incidence in pancreatic adenocarcinoma. With varied suggestions on ASCP prevalence, LY6D expression can now be used to quickly identify ASCP and allow for more functional experimentation to better understand differences with PAAD. Roughly 11 of the 24 (45.8%) scRNA-Seq tumor datasets analyzed here have higher levels of LY6D, suggesting that ASCP might be more widespread than initially thought. Additionally, it is possible that, like the BCC-SCC epithelial plasticity spectrum, pancreatic cancer may have a more plastic nature that toggles between a non-squamous and squamous state. EGFR inhibitors, which reduced LY6D+ accumulation in our BCC system, are approved to treat PDAC, suggesting that they might be more effective in LY6D-expressing ACSP10,28. Overall, this work demonstrates the utility of LY6D as a generalizable marker of a basosquamous state with poor clinical outcomes.