Phenotypic dissection of epithelial lineages and therapeutic manipulation of differentiation programs in human Adenoid Cystic Carcinomas (ACCs)

Adenoid Cystic Carcinoma (ACC) is a rare and aggressive form of salivary gland cancer, characterized by the co-existence within tumor tissues of two distinct populations of malignant cells, phenotypically similar to the myoepithelial and ductal lineages of normal salivary glands. Using a novel computational approach for single-cell RNA-seq analysis, we identified two cell-surface markers (CD49f, KIT) that enable the differential purification of myoepithelial-like (CD49fhigh/KITneg) and ductal-like (CD49flow/KIT+) cells from ACC patient derived xenografts (PDX). Using prospective xeno-transplantation experiments, we demonstrate that myoepithelial-like cells act as progenitors of ductal-like cells. Using three-dimensional (3D) organoid cultures, we demonstrate that agonists of retinoic acid (RA) signaling promote differentiation of myoepithelial-like cells into ductal-like cells, while inhibitors of RA signaling selectively kill ductal-like cells. Finally, we demonstrate that BMS493, an inverse agonist of RA signaling, can be successfully leveraged for the in vivo treatment of human ACCs.


INTRODUCTION
Adenoid cystic carcinomas (ACCs) are malignant adenocarcinomas which arise in secretory glands throughout the body, most commonly in the salivary glands (SGs) of the head and neck 1 . ACCs usually display an indolent growth, but their slow proliferation kinetics belie their biologically aggressive and clinically relentless nature, characterized by high rates of local invasion, peri-neural infiltration, and high frequency of distant metastases, often 10-15 years after initial diagnosis 2 .
Current treatment options for ACCs are limited to surgery and radiotherapy, which are often destructive and, in approximately 60% of cases, unable to prevent metastatic relapse and patient death.
Due to their slow proliferation kinetics and low rates of somatic mutation (i.e., low mutational burdens), ACCs are usually refractory to chemotherapy and immunotherapy regiments 3,4 and lack meaningful targets for the development of pharmacological agents with tumor-specificity (i.e., kinase inhibitors) 5,6 .
From a histological point of view, ACCs are characterized by a defining feature: the concurrent presence of two phenotypically distinct populations of cancer cells, termed "ductal-like" and "myoepithelial-like" because of their similarity to the ductal and myoepithelial cell lineages of the normal SG epithelium, in terms of both morphology and expression of lineage-specific biomarkers 1, [7][8][9] .
Although this feature has long been appreciated, the molecular mechanisms that underpin it remain poorly understood. It remains unclear, for example, whether the two cell types originate as a result of clonal events (i.e., the accumulation of divergent somatic mutations) or of multi-lineage differentiation processes akin to those that enable stem/progenitor cells to sustain the homeostatic turnover of adult tissues 10 , as in the case of colon cancer 11 . In the case of ACC, this important aspect of the tumor's biology has remained difficult to investigate, due to the lack of experimental means that would enable the isolation of the two cell populations from primary tumors, and allow their comparison in terms of transcriptional identities, functional properties, developmental relationships, and differential sensitivity to anti-tumor therapies.
In this study, we utilize single cell RNA sequencing (scRNA-seq) in concert with a novel computational approach derived from random matrix theory (RMT) 12 , to identify cell surface markers that enable, for the first time, the differential purification by fluorescence activated cell sorting (FACS) of the two malignant cell populations found in patient-derived xenograft (PDX) lines established from human ACCs 13 . We then proceed to test their developmental relationships using prospective transplantation experiments, to identify a signaling pathway that controls the differentiation of one into the other, and to develop a new therapeutic approach for the in vivo treatment of human ACCs.

Identification of cell-surface markers for the differential purification of myoepithelial-like and
ductal-like cell types from human ACCs. To identify cell-surface markers differentially expressed between the myoepithelial-like and ductal-like populations of human ACCs, we performed scRNA-seq on the PDX line ACCX22. We selected this model because it was established from a parotid gland ACC lesion, representing one of the most frequent sites of origin for ACCs 8 , and because it displayed two classical morphological features of human ACCs: a well-differentiated "cribriform" histology and a biphenotypic cell composition. Utilizing our previously established protocols for dissociation and singlecell transcriptional analyses of PDX tumors 11,14 , we dissociated ACCX22 into a single-cell suspension, purified human epithelial (EpCAM + ) tumor cells by FACS, and used scRNA-seq to capture their individual transcriptional profiles (Extended Data Fig. 1). We then implemented a novel computational pipeline (Randomly) 12 designed to leverage the statistical principles of RMT to remove gene-expression signals that follow a stochastic distribution and identify genes most responsible for biological signal. We clustered cells in different subgroups based on transcriptional patterns (Extended Data Fig. 2), with the optimal clustering solution (i.e., the one displaying the highest silhouette score) consisting of three subpopulations (Figure 1a). The two most abundant populations displayed mutually exclusive expression of known myoepithelial-like (e.g., ACTA2, CNN1, TP63) and ductal-like (KRT7, KRT18, ELF5) markers (Extended Data Fig. 3a-b), while a third population appeared to represent a highly proliferating (MKI67 high , TOP2A high , CDK1 high , PCNA high ) subset of the population expressing ductal-like markers ( Figure 1e, Extended Data Fig. 3c). Among the genes displaying the highest differential expression between the two main clusters, we identified cell-surface markers CD49f (ITGA6), and KIT/CD117 (KIT), which were preferentially detected in cells associated with myoepithelial-like and ductal-like markers, respectively (Figure 1b-d). We then tested whether the expression levels of CD49f and KIT proteins could be leveraged to visualize the two subsets of malignant cells using FACS. Indeed, the combined use of the two markers allowed us to discriminate two clearly distinct cell populations (CD49f high /KIT neg vs. CD49f low /KIT + ), across 5 independent PDX lines representative of bi-phenotypic ACCs (ACCX5M1, ACCX6, ACCX14, ACCX22, SGTX6) (Figure 1f). Analysis of the same tumors by immunohistochemistry (IHC) confirmed that KIT expression was restricted to cells with ductal-like morphology, and mutually exclusive to expression of TP63, a marker of myoepithelial-like cells, in agreement with previous reports 9 (Figure 1g). The transcriptional profiles of CD49f high /KIT neg and CD49f low /KIT + cells are highly conserved across ACCs from different patients. To understand whether CD49f high /KIT neg and CD49f low /KIT + cells isolated from ACCs of different patients displayed similar gene-expression patterns, we used FACS to sort autologous pairs of cells from all 5 bi-phenotypic PDX lines (ACCX5M1, ACCX6, ACCX14, ACCX22, SGTX6), and then analyzed them by conventional bulk RNA-seq. This allowed us to obtain robust transcriptional profiles of the two populations, in which the representation of genes would not be compromised by technical limitations (i.e., low capture efficiency, gene dropouts) intrinsic to scRNAseq techniques.
As a first step, we performed a Principal Component Analysis (PCA) of the 10 samples included in the RNA-seq dataset (Figure 1h). The analysis revealed that, indeed, the 10 samples grouped tightly into two separate and distinct clusters (5 samples/cluster), exactly corresponding to the two cell phenotypes.
Importantly, the analysis also showed that the two clusters separated along the first principal component (PC1), which accounted for a dominant fraction (58%) of the variability in gene-expression levels observed within the dataset. Taken together, these observations indicated that CD49f high /KIT neg and CD49f low /KIT + cells are defined by systematic differences in transcriptional profiles, which are strongly conserved across tumors originated from different patients. Furthermore, transcriptional differences between the two cell types are greater in magnitude than those attributable to variables which are unique to the individual tumors (e.g., site of origin, repertoire of somatic mutations, genetic background of patients) (Supplementary Table 1).
To understand which genes were differentially expressed between the two cell phenotypes in the most robust and systematic manner, we performed a differential expression (DE) analysis using the DESeq2 package 15 . The analysis identified 643 genes that displayed more than a 2-fold difference in expression levels between the two populations, with high statistical significance (p<0.05, after Benjamini-Hochberg adjustment, Supplementary Table 2). We then ranked all differentially expressed genes based on their adjusted p-value, and selected the top 100 for visualization by hierarchical clustering (Figure 1i). Among the genes preferentially expressed by CD49f high /KIT neg cells were several wellestablished markers of the myoepithelial lineage (e.g., ACTA2, MYH11, PDPN, TP63) 16, 17 . Among the genes preferentially expressed by CD49f low /KIT + cells were several well-established markers of the ductal and luminal lineage in ectodermal secretory glands (e.g., ELF5, KIT, SLPI, ANXA8) [18][19][20][21] . Taken together, our observations confirmed that CD49f and KIT represented robust markers for the differential purification of myoepithelial-like (CD49f high /KIT neg ) and ductal-like (CD49f low /KIT + ) populations from human ACCs.
In parallel, we also tested whether CD49f high /KIT neg and CD49f low /KIT + cells differed in the expression levels of the MYB and NFIB genes, or MYB-NFIB chimeric transcripts, detected using the STAR-Fusion software 22 . A majority (50-60%) of human ACCs contain t(6;9) MYB-NFIB chromosomal translocations, which cause constitutive upregulation of MYB expression, and are recognized as primary oncogenic drivers in ACC cancerogenesis 5,6,[23][24][25][26] . Although it is established that, when present, the MYB-NFIB translocation is found in both myoepithelial-like and ductal-like cells 25 , it remains unknown whether the two cell types actually differ in terms of expression of the corresponding chimeric transcript.
Our analysis found no significant difference in the expression level of MYB or NFIB between the two populations across the 5 analyzed tumors (Extended Data Fig. 4a-b). Out of the 5 PDX models with a bi-phenotypic histology, 3 (SGTX6, ACCX5M1, ACCX14) contained MYB-NFIB fusion transcripts which gave rise to chimeric transcripts, in agreement with previous results from whole genome and Sanger sequencing studies 27,28 . In these models, fusion transcripts were representative of multiple isoforms, originated from the alternative splicing of the 3' end (NFIB) of the chimeric gene 23,24 , and displayed no statistical difference in expression levels between the two populations (Extended Data   Fig. 4c, e).
CD49f high /KIT neg cells display a higher tumor-initiating capacity as compared to CD49f low /KIT + cells. To understand whether the two populations display functional differences in terms of their capacity to initiate and sustain tumor growth, we compared their tumor-initiating capacity upon prospective xenotransplantations by Extreme Limiting Dilution Analysis (ELDA) 29 . We used FACS to purify autologous pairs of CD49f high /KIT neg and CD49f low /KIT + cells from bi-phenotypic PDX lines and injected them, side-by-side, in the subcutaneous tissue of NOD/SCID/IL2Rγ -/-(NSG) mice, at progressively decreasing doses (10,000-250 cells/mouse) (Figure 2a). The results showed that the frequency of tumorigenic cells was substantially higher in CD49f high /KIT neg as compared to CD49f low /KIT + cells. As the dose of injected cells progressively decreased, the frequency of tumors formed by CD49f low /KIT + cells was rapidly reduced as compared to the frequency of tumors formed by the CD49f high /KIT neg population, until it was eventually lost (Figure 2d-e, Extended Data Fig. 5). To understand whether the higher tumor-initiating frequency observed among CD49f high /KIT neg cells could be attributed to higher proliferation rates, we analyzed their cell-cycle distribution and compared it to that of autologous CD49f low /KIT + cells, across all 5 bi-phenotypic PDX lines (Extended Data Fig. 5a-d). We found that CD49f high /KIT neg cells had a lower proportion of cells in the G2/M phase of the cell cycle compared to CD49f low /KIT + cells (Extended Data Fig. 5g-h), in agreement with our scRNA-seq data, which had revealed a sub-group of cells expressing proliferation markers within the ductal-like (KIT + ) population (Figure 1e, Extended Data Fig. 3). Taken together, our results indicated that myoepithelial-like cells represent an aggressive component of the malignant tissue, capable of initiating and sustaining tumor growth, despite their low proliferation kinetics.
CD49f high /KIT neg cells and CD49f low /KIT + cells do not represent distinct genetic clones, but emerge as a result of epigenetic differentiation. We next wanted to understand whether the two cell populations represented different genetic clones, originated from the accumulation of distinct somatic mutations, or different developmental lineages, originated through differentiation processes akin to those that enable stem cells to sustains the homeostasis and regeneration of multi-cellular tissues. We hypothesized that, in a "genetic/clonal" scenario, distinct cell types would maintain their phenotypic identity upon transplantation (Figure 2b). In contrast, in an "epigenetic/differentiation" scenario, one or more cell types could serve as a progenitor of the others, in a plastic and dynamic fashion (Figure 2c).
When we analyzed the cell composition of tumors originated from the ELDA experiments, we found that tumors derived from CD49f high /KIT neg cells recapitulated the cribriform histology and bi-phenotypic cell composition, with both cell populations (CD49f high /KIT neg and CD49f low /KIT + ) present at identical ratios to those observed in the parent lines (Figure 2f-g, i-j). This finding indicated that CD49f high /KIT neg cells can differentiate into CD49f low /KIT + cells, thus excluding the "genetic/clonal" hypothesis. When we analyzed the limited number of tumors originated from high doses of CD49f low /KIT + cells, we also found them to be indistinguishable from parent lines (Figure 2h, k). In the specific case of CD49f low /KIT + cells, however, given the higher number of cells required to initiate tumor growth, we could not exclude the possibility that tumors might have arisen from a small cross-contaminations of CD49f high /KIT neg cells, despite the high levels of purity achieved by double-sorting.
CD49f high /KIT neg and CD49f low /KIT + cells are characterized by differential activation of the retinoic acid (RA) signaling pathway. To shed light on the molecular mechanisms involved in regulating the differentiation of CD49f high /KIT neg cells into CD49f low /KIT + cells, we decided to investigate the RA signaling pathway. RA signaling is critical for proper morphogenesis and differentiation of salivary gland tissues during development [30][31][32] . Furthermore, recent studies have demonstrated that stimulation of RA signaling can antagonize MYB signaling and slow tumor growth in ACC PDX lines 33,34 . We first examined whether genes encoding effectors and modulators of RA signaling, including enzymes necessary for RA biosynthesis 30 , RA binding proteins 35,36 , and RA receptors 37 , were differentially expressed between CD49f high /KIT neg cells and CD49f low /KIT + cells isolated from PDX lines (Figure 3a). Our analysis revealed that, indeed, most of these genes displayed statistically significant differences in expression levels (Figure 3b). More importantly, when we categorized such genes as potentiators 30,32,36 or attenuators 38,39 of RA signaling, we found that potentiators were expressed at higher levels in CD49f low /KIT + cells, while attenuators were expressed at higher levels in CD49f high /KIT neg cells, in a coordinated and systematic fashion (Figure 3c).

Activation of the RA signaling pathway induces differentiation of CD49f high /KIT neg cells into
CD49f low /KIT + cells, while inhibition of the RA signaling pathway causes selective death of CD49f low /KIT + cells. To formally test whether RA signaling controls cell differentiation in ACCs, we utilized a three-dimensional (3D) in vitro organoid culture system. ACC organoid cultures retained the bi-phenotypic composition of PDX models, and thus allowed us to quantify the effects on cell composition caused by pharmacological perturbations of RA signaling (Figure 3d-g) 40,41 . We tested whether activation or inhibition of RA signaling altered the cell composition of ACC organoids, using direct and inverse agonists. We found that agonism of RA signaling with increasing doses of either all- To formally dissect the effects of activators and suppressors of RA signaling on the two cell populations, we purified CD49f high /KIT neg and CD49f low /KIT + cells by FACS and then treated them individually with either ATRA (10 µM) or BMS493 (10 µM) using in vitro monolayer cultures 42 (Figure   4m, q). The experiment revealed that agonism of RA signaling (i.e., treatment with ATRA) did not cause cell death of either CD49f high /KIT neg or CD49f low /KIT + cells (Figure 4n), but induced a dramatic change in phenotype in CD49f high /KIT neg cells, which became almost completely CD49f low /KIT + (Figure 4o).
This result indicated that, the effects observed in organoid cultures stimulated with RA agonists resulted from differentiation of CD49f high /KIT neg cells into CD49f low /KIT + cells (Figure 4p). Conversely, the experiment also demonstrated that suppression of RA signaling (i.e., treatment with BMS493) resulted in fragmentation and death of the majority of CD49f low /KIT + cells (Figure 4r-s), indicating that the effects observed in organoid cultures following treatment with RA inverse agonists were likely caused by selective cell death of CD49f low /KIT + cells (Figure 4t).

Inverse agonists of RA signaling can be leveraged as anti-tumor agents against human ACCs,
with preferential activity against "solid" variants of the disease, which are exclusively composed of CD49f low /KIT + cells. We next aimed to elucidate whether the effects observed in vitro following treatment with inverse agonists of RA signaling could be leveraged for in vivo therapy targeting human ACCs. Because we observed that inverse agonists of RA signaling had selective/preferential cytotoxicity against CD49f low /KIT + cells, we hypothesized that tumors containing a large percentage of KIT + cells would represent the most susceptible targets. Although ACCs usually present with a bi-phenotypic histology, over the natural history of the disease a subgroup of ACCs progresses to a "solid" histological pattern, predominantly composed of ductal-like (KIT + ) cells. Recent studies have shown that ACC progression to solid histology associates with activating mutations in NOTCH1 5, 6, 27, 43, 44 , higher proliferation kinetics and a more aggressive clinical course 9,45,46 .
To understand whether ACCs with solid histology indeed represented mono-phenotypic expansions of ductal-like cells akin to those found in bi-phenotypic ACC tumors, we analyzed two PDX lines with solid histology and NOTCH1 activating mutations (ACCX9, ACCX11) 27 . Indeed, analysis of the two PDX lines by FACS revealed that, in both cases, the bulk of the malignant cells consisted of a single population characterized by high expression levels of KIT (Figure 5a-b). Analysis by IHC confirmed that both lines consisted almost exclusively of KIT + cells, and lacked cells expressing the myoepithelial marker TP63, in agreement with previous reports 27 (Figure 5c-f). We then performed bulk RNAsequencing of KIT + cells sorted by FACS from both models, in order to obtain a better understanding of their transcriptional profiles. When we performed a PCA analysis, this time combining the two KIT + populations purified from solid ACCs with the 5 paired sets of CD49f high /KIT neg and CD49f low /KIT + cells isolated from bi-phenotypic ACCs, we observed that, indeed, the KIT + cells from the two solid ACCs clustered with the CD49f low /KIT + cells from the bi-phenotypic tumors, indicating that they retained a transcriptional profile characteristic of ductal-like cells (Figure 5g). A similar result was also obtained when we performed a hierarchical clustering of the same 12 populations, using the top 100 genes identified as the most differentially expressed between CD49f high /KIT neg and CD49f low /KIT + cells ( Figure 5h).
We then wanted to understand whether, in addition to sharing transcriptional similarities with the CD49f low /KIT + cells of bi-phenotypic ACCs, KIT + cells from solid ACCs also shared a vulnerability to pharmacological suppression of RA signaling. We generated organoid cultures from both PDX lines with solid histology, and then treated them in vitro with BMS493 (10 µM). After one week of culture, organoids treated with BMS493 began to lose structural integrity, appearing fragmented and apoptotic As a final step, we wanted to test whether an inverse agonist of RA signaling (BMS493) could be leveraged for the in vivo treatment of human ACCs with either solid (mono-phenotypic) or cribriform (bi-phenotypic) histology. We engrafted three PDX lines, two representative of solid ACCs (ACCX9, ACCX11) and one representative of a cribriform ACC (ACCX5M1), in NOD/SCID/IL2Rγ -/-(NSG) mice. We then let tumors grow up to a volume of approximately 0.4 cm 3 , and then treated tumor-bearing mice with BMS493 (40 mg/kg, i.p.) using a slightly more intense regimen for the cribriform models (4 times/week x 3 weeks) as compared to the solid models (3 times/week x 3 weeks), under the assumption that the cribriform model would be more resistant to the treatment. As expected, treatment with BMS493 was associated with side-effects reminiscent of those previously described as secondary to vitamin A deficiency in rodents (e.g., encrusted eyelids, rough coat, scaling of skin) 47 , which began to manifest during the second week of treatment and worsened during the third. Overall, out of 18 tumor-bearing animals treated with BMS493, 33% (n=6/18) experienced reductions in tumor volume (Extended Data   Fig. 8). Four animals (4/18 22%) had to be prematurely euthanized due to a sudden deterioration in general health; in three of these animals, euthanasia followed tumor shrinkage, suggesting a possible acute complication due to widespread tumor lysis (Extended Data Fig. 8). Overall, treatment with BMS493 was associated with a statistically significant reduction in the growth kinetics of engrafted tumors across all three tested models, even after removal of the four mice undergoing premature euthanasia (Figure 6c, e, h, Extended Data Fig. 8). Taken together, our data revealed that BMS493 has robust anti-tumor activity against human ACCs, both in vitro and in vivo, and indicate that inverse agonist of RA signaling are deserving of further study as a new class of anti-tumor agents for the clinical treatment of salivary gland malignancies.

DISCUSSION
The results of our study reveal that the two major sub-types of malignant cells (myoepithelial-like and ductal-like), which co-exist within human ACCs, can be differentially purified by FACS using two surface markers (CD49f, KIT). We find that these two markers are differentially expressed between ACC cell populations in a robust and reproducible manner, and across multiple tumor models. This technical advancement has important implications, as it enables, for the first time, experiments that are aimed at comparing the functional properties of the two populations in "side-by-side" prospective assays, including experiments designed to elucidate their developmental relationships.
The results of our gene-expression studies, conducted on autologous pairs of cell populations sorted in parallel from the same tumors, indicate that the two populations are defined by transcriptional signatures that are closely recapitulated across different patients. The two cell-specific transcriptional programs appear to represent one of the most dominant sources of intra-tumor heterogeneity, irrespective of each tumor's mutational repertoire, its site of primary origin, or the patient's sex. The results of our xeno-transplantation studies in immune-deficient mice, designed to compare the tumorigenic capacity of the two cell populations, show that the two populations do not represent distinct genetic clones (i.e., they do not originate from the divergent accumulation of independent somatic mutations), but distinct epigenetic lineages (i.e., different cell types originated through the dynamic and systematic activation of "hard-wired" developmental programs, akin to those that enable stem cell populations to support the homeostatic turnover and/or regeneration of normal tissues) 10 . In agreement with these observations, we also find that the relative proportion of myoepithelial-like and ductal-like cells within ACC tumors represents a stable, intrinsic feature of each PDX line, as each PDX model maintains the same relative frequencies over serial passaging, and most notably, in tumors originated from purified cell populations.
From a translational point of view, one the most important findings of our study is the observation that, despite their apparent low proliferation rates, myoepithelial-like cells (CD49f high /KIT neg ), represent a highly tumorigenic component of human ACCs, capable of supporting the formation of tumors when transplanted in immune-deficient mice, and acting as a progenitors of the more proliferative ductal-like cells (CD49f low /KIT + ). Traditionally, in tumors originated from secretory glands of ectodermal origin (e.g. breast cancer), myoepithelial-like components are considered to represent a relatively indolent subset of malignant cells, often used as a biomarker of tumors with a more benign clinical behavior, and postulated to exert a suppressive effect on metastatic spread 48,49 . Our findings caution against this interpretation in the case of the majority of human ACCs, especially those originated from salivary glands. Indeed, our findings indicate that, in order to be curative, treatment strategies for ACCs will likely need to eradicate both myoepithelial-like and ductal-like cell populations.
From a more mechanistic point of view, our study also demonstrates that, in human ACCs, the capacity of myoepithelial-like cells to differentiate into ductal-like cells is controlled by RA signaling. Recent studies have identified RA signaling as a critical pathway in the biology of human ACCs, demonstrating that RA agonists, such as ATRA, have powerful anti-proliferative activity in PDX models 33,34 . A recent clinical trial, however, showed that treatment with ATRA was associated with limited clinical benefit in ACC patients 50 . Our in vitro studies, using 3D organoids and purified cell preparations of both myoepithelial-like and ductal-like cells, provide novel mechanistic insights into the anti-tumor activity of retinoids in ACCs, and a possible explanation for the conflicting results between pre-clinical and clinical data. We find that activation of RA signaling induces differentiation of myoepithelial-like cells into ductal-like cells, without compromising their viability, in agreement with the well-known function of ATRA as a differentiation agent across multiple tissues 51 . We therefore hypothesize that, in a clinical setting, the therapeutic benefit of ATRA might be short-lived, because of the transient and dynamic nature of the perturbations that it is predicted to induce in the cell composition of ACC tumors.
Perhaps more importantly, we found that inhibition of RA signaling induces selective cell death of ductal-like cells. This finding provides an opportunity for the development of novel therapeutic agents against ACCs, especially in tumors with solid histology, which are characterized by a mono-phenotypic expansion of ductal-like cells. These tumors usually originate during the progression of human ACCs, with kinetics that are reminiscent of the "blast crisis" phase of human chronic myelogenous leukemias (CMLs), in which a population of more differentiated, yet highly proliferating myeloid progenitors becomes dominant, as a result of mutations that aberrantly activate the pathways that confer self-renewal to normal hematopoietic stem cells 10,52 . Our pre-clinical data show that, in solid ACCs, inhibition of RA signaling using inverse agonists (e.g., BMS493) can have profound and robust anti-tumor activity.
Unlike agonists of RA signaling, which have been extensively explored as anti-tumor agents in a number of human malignancies [53][54][55] and have known and manageable associated toxicities 56 , inverse agonists of RA signaling, to our knowledge, have not yet been evaluated for clinical use. Because RA signaling is critical to the biology of many tissues, the anti-tumor benefit of RA antagonists will have to be weighted carefully against any potential toxicities, based on the results of future toxicology and pharmacodynamics studies.  Fig. 2b-c). The final output produced by Randomly corresponds to the latent space, which is then analyzed to identify cell populations, using synergistic machine-learning techniques for cell clustering.
In our experiment, clustering was performed using the Leiden algorithm as implemented in the code by Wolf et al. 59 (https://scanpy.readthedocs.io/en/stable). The optimal number of clusters was selected based on their mean silhouette score 60 . More precisely, we performed a set of clusterings using different Leiden resolutions, then computed the mean silhouette score for each case, and finally selected the clustering solution with the highest mean silhouette score (Extended Data Fig. 2d).

Visualization of scRNA-seq results and differential expression (DE) analysis. Dimensionality reduction and Uniform Manifold Approximation and Projection (UMAP) representation of the scRNA-
seq data were performed using the visualization functions of the Randomly package, with default parameters. Genes with differential expression between the myoepithelial-like cell-population (Cluster 1) and the ductal-like cell population (Cluster 2) were identified based on the following two-step approach: 1) of all the genes identified as informative based on their projection in the 47 eigenvectors corresponding to biological signal, those displaying preferential expression in one of the two cell populations (i.e., expressed in at least 60% of cells in one population, and no more than 20% of cells the other) were selected for analysis; 2) genes selected for analysis were tested for statistical significance of their difference in mean expression levels between the two populations (Student's t-test, two-tailed). PCA was performed using the plotPCA function with default parameter for the number of top genes to use for principal components (i.e., 500 genes). Differentially expressed genes were defined as those displaying a more than two-fold difference in mean expression levels between the two populations (log2 fold-change > 1) that was considered to be statistically robust after ranking based on the Wald statistic and correction for multiple comparisons (FDR<0.05; Benjamini-Hochberg method). The variance in gene-expression levels across independent samples was visualized using heatmaps, generated using the pheatmap function, with scaling performed by mean-centering the gene expression values for each gene.
Heatmaps were organized by hierarchical clustering of both genes and samples.  Fig. 2b-c). Double-sorting consisted in two sequential rounds of sorting, whereby, after the first sort, cells were spun down, resuspended in 0.5 mL of fresh FCB with DAPI, and then sorted a second time, using identical gates. Cells were assessed for purity and viability after the second sort, resuspended in fresh FCB and counted using a hemocytometer. Cells were then aliquoted at various doses (250- Tumors originated from the injection of purified preparations of either CD49f high /KIT neg or CD49f low /KIT + cells were analyzed by flow-cytometry and IHC, to evaluate their cell composition.

Identification and quantification of MYB-NFIB
Differences in the ratio of CD49f high /KIT neg cells and CD49f low /KIT + cells observed between parent tumors and tumors generated from purified populations were tested for statistical significance using an ordinary one-way ANOVA with Dunnett's multiple comparisons test.  , 10 mM), were prepared in DMSO and stored at -20ºC, protected from light. On the day of use, stock solutions were thawed and added to complete organoid medium, at appropriate concentrations (100 nM -10 µM). Due to the short half-life of retinoids, medium with retinoid compounds was kept for a maximum of 3 days at 4ºC and changed daily for the duration of culture (7 days). cell composition (histology, FACS, brightfield microscopy). Organoid cultures established from human ACCs were dissociated from Matrigel by incubation with a solution of 2 mg/mL Dispase-II (Thermo Fisher, 17105041) and 200 U/mL Collagenase-III (Worthington, LS004183) in DPBS at 37ºC for 15 minutes. Organoid suspensions were then transferred to 1.5 mL plastic tubes and pelleted by centrifugation (10,000 RPM, 2 minutes). Excess Matrigel and disaggregation solution were carefully aspirated. To dissolve remaining Matrigel, organoid pellets were briefly (3 minutes) resuspended in 0.25% Trypsin at 37ºC, followed by a wash step with cell culture medium containing 10% FBS. To prepare organoids for immunohistochemistry, organoids were pelleted by centrifugation and fixed in 10% formalin for 4-12 hours. Fixed organoids were embedded in paraffin blocks, from which 4 µm tissue-sections were cut and stained following protocols identical to those used for tumor tissues (described above). To prepare organoids for FACS, organoid pellets were resuspended in disaggregation medium containing DNase-I (100 U/mL), Collagenase-III (200 U/ml), and Hyaluronidase (100 U/ml) and incubated at 37 ºC for 20-30 minutes. Disaggregated organoids were then pelleted and incubated in 0.25% Trypsin for 10-15 minutes to generate a single cell suspension.

Analysis of organoid
Dissociated cells were washed with cell culture medium containing 10% FBS to inhibit Trypsin activity, followed by blocking with human IgGs (5mg/ml) and staining with antibodies. Differences in the percentage of CD49f high /KIT neg and CD49f low /KIT + cells in ACC organoids treated with retinoid compounds were tested for statistical significance using a one-way ANOVA with Dunnett's multiple comparisons test. Brightfield images of organoid cultures were acquired using a Cytation 5 Cell Imaging Reader (BioTek) at 4x magnification. Images of hematoxylin and eosin (H&E) or IHC-stained organoids were acquired using a Nikon Eclipse E600 microscope with NIS-Elements Software (version 5.21).
Brightness and contrast were adjusted uniformly throughout whole images using Adobe Photoshop at individual time points were tested for statistical significance using a two-tailed Student's t-test, while differences between growth rates (i.e., Log10(fold-increase in tumor volume)/Time(days)) were tested for statistical significance using a two-tailed Welch's t-test (i.e., assuming unequal variance).        (g-h) Comparison of tumor growth kinetics, expressed as either fold-increases in tumor volumes or growth rates, between mice treated with BMS493 and mice treated with the drug's vehicle alone (DMSO), following subcutaneous engraftment with a bi-phenotypic ACC model (ACCX5M1). Differences between mean fold-increases in tumor volume were tested for statistical significance using a two-tailed Student's t-test (ns = non-significant, *p<0.05, **p<0.01, ***p<0.001). Differences between mean growth rates were tested for statistical significance using a two-tailed Welch's t-test. Growth rates were calculated assuming exponential kinetics. Error bars: mean +/-standard deviation. scRNA-seq. Single-cell libraries were prepared using the 10x Chromium system (Single Cell 3' v3 chemistry), and sequenced using the NovaSeq 6000 platform (Illumina). Single-cell data was filtered to retain cells with >500 genes/cell and normalized (log2(1+TPM)). Gene-expression data was processed to remove signals attributable to stochastic variation (noise) using the Randomly package, then used to sub-group cells into distinct clusters, which were visualized using Uniform Manifold Approximation and Projection (UMAP). Schematic in panel (c) was created with BioRender.com tumors treated with DMSO (gray) or BMS493 (magenta). Cross symbols ( †) denote two animals who were sacrificed early, due to deterioration of their health condition. (b) Growth rates of ACCX9 tumors treated with either DMSO (n=6) or BMS493 (n=7). Two treatment cohorts are identified by different symbols (circles = cohort 1; triangles = cohort 2). Differences in mean growth rates were statistically significant (Welch's t-test; **p<0.01). Black symbols denote the two animals who were sacrificed because of a deterioration of their health condition. (c) Individual animal weights over the course of in vivo treatment with BMS493 ( †: animals sacrificed due to deterioration of health conditions). (d) Individual growth curves of ACCX11 tumors treated with DMSO (gray) or BMS493 (magenta). (e) Growth rates of ACCX11 tumors treated with either DMSO (n=4) or BMS493 (n=5). Differences in mean growth rates were statistically significant (Welch's t-test; **p<0.01).
(f) Individual animal weights over the course of in vivo treatment with BMS493. (g) Individual growth curves of ACCX5M1 tumors treated with DMSO (gray) or BMS493 (magenta). Cross symbols ( †) denote two animals who either were sacrificed early due to deterioration of their general health or who died at the final time point. (h) Growth rates of ACCX5M1 tumors treated with either DMSO (n=6) or BMS493 (n=6). Two treatment cohorts are identified by different symbols (circles = cohort 1; triangles = cohort 2). Differences in mean growth rates were statistically significant (Welch's t-test; *p<0.05). Black symbols denote the two animals who either were sacrificed early due to deterioration of their general health or died at the final time point. (i) Individual animal weights over the course of treatment ( †: animals sacrificed or found dead).