Single-cell RNA sequencing of salivary gland tumors

To identify and characterize the cellular heterogeneity that is specific to the solid tumor context, we first established controlled ways to dissociate tumor-bearing (double-mutant: β-cat^GOF; Bmpr1a^LOF) and control salivary glands into high-quality single-cell suspensions (Fig. 1a). After dissociation, dead cells and enucleated cellular debris were excluded and live intact cells obtained by fluorescence-activated cell sorting (FACS) (Supplementary Fig. 1a). Cells were directly sorted into methanol for fixation¹⁸, and further processed to profile their transcriptomes by a high-throughput droplet-based approach (Drop-Seq)¹⁹. In total, 26 single-cell RNA libraries were generated from 12 control and 14 double-mutant (tumor-bearing) salivary glands of either female or male mice from an early and a late tumor stage at postnatal days 40 (P40) and 90 (P90), respectively (Fig. 1a). To validate our experimental approach, we compared all single-cell samples, computationally pooled by disease status (control or double-mutant), to bulk mRNA-seq data that were generated from equivalent, freshly dissected but unprocessed, salivary glands (Supplementary Fig 1b). Although gene expression levels correlated better within experimental procedures and samples grouped by genotype, correlations between all samples were generally high (R ≥ 0.74). Additionally, comparison of global transcript counts (Supplementary Fig. 1c) show that individual single-cell RNA libraries correlated well to each other (R ≥ 0.8) with no apparent bias towards the disease-, sex-, stage-related status or the experimental batch in which a specific sample was processed. After computational cell selection and filtering, we obtained a total of ~23,000 cells from 26 individual salivary glands (Supplementary Fig. 1d) and typically detected a median of ~500 genes and ~1000 unique molecular identifiers (UMIs) per cell (Supplementary Fig. 1e, f). To quantify sample-to-sample variation and the extent of possible batch effects, we developed an entropy-based approach²⁰ to measure how evenly cells’ nearest neighbors are distributed among different samples. This analysis showed that although these distributions were not completely random, ~80% of the variance was explained by the three major biological variables (genotype, sex, and stage) (Supplementary Fig. 2a). Further, cell distributions were largely balanced except in cell types that are affected by one of these biological variables (Supplementary Fig. 2b). We thus conclude that our protocol is reproducible given biologically-relevant sample differences.

A comprehensive salivary gland cell atlas

In order to provide a comprehensive cell atlas of the salivary gland that could also serve as an appropriate reference to the tumor context, we first pooled and analyzed single-cell datasets from control salivary glands using Seurat²¹ (Fig. 1b). After inspection of marker genes, we assigned the cell type identity to clusters, for which the expression of specific genes has been reported in previous literature^{22,23,24,25,26,27,28,29,30,31}, and molecularly characterized and validated additional cell types for which no or only ambiguous information was available (Fig. 1b, c; Supplementary Figs. 3, 4). In agreement with previous studies^32,33,34, we observed that Egf and Smgc were expressed in cell populations with strongly sex-dependent representation (Supplementary Fig. 3) and identified several other marker genes with similar patterns (Supplementary Fig. 4a). Interestingly, we also noted that Dcpp1+ cells were more abundant in tissues from mice at the P40 than at the P90 stage and validated this finding by immunofluorescence in tissue sections (Supplementary Fig. 3b, c). Altogether, this indicates that we can reliably identify and characterize cell types at high resolution and provide evidence that the ductal composition of healthy salivary glands is sexually dimorphic and stage-dependent. Our approach combined with imaging-based validations enabled us to chart all epithelial cells onto a consolidated anatomical sketch depicting female- (Fig. 1c) and male-specific (Supplementary Fig. 3e) features.

Identification of cancer stem cells and other tumor-specific cell populations

To systematically uncover cells that were specific to the tumor samples, we analyzed and clustered cells pooled from all single-cell datasets of control and double-mutant samples (Fig. 2a). This combined analysis enabled us to recapitulate all previously identified cell types without the need to resort to advanced sample alignment methods²¹ as cells from both control and double-mutant mice were distributed evenly on the tSNE for many clusters as shown in the local density plot (Fig. 2b). However, in line with our entropy-based analysis (Supplementary Fig. 2), we noted that a number of cell types, including stromal and immune-related cells were significantly more abundant in double-mutant (Supplementary Fig. 5a, Supplementary Fig. 6), while acinar and some ductal cell types were accordingly more prevalent in control tissues (Supplementary Fig. 5b).

a Clustering of combined single-cell sequencing data from control and double-mutant (β-cat^GOF; Bmpr1a^LOF) tumor-bearing salivary glands shown in the tSNE. Tumor-specific epithelial cell clusters are shown in the lower right part of the tSNE in shades of pink: cancer stem, basal, luminal Clu+ cells. b Same tSNE plot showing the relative local densities of control and double-mutant cells. Color gradients reveal clusters that are predominantly represented by double-mutant (red), control (blue) or both double-mutant and control cells (gray). c Expression of Axin2, Ptn, Wif1, Clu, and Wfdc18 in the tumor-specific epithelial subset of the tSNE representation as indicated by a color scale ranging from gray (no expression) to maroon (high expression).

In particular, our analysis revealed clusters of epithelial cells with transcriptional profiles unique for the tumor context. This included luminal- and basal-like cells as well as a small but distinct cancer stem cell (CSC)-like population in which Wnt-specific genes were activated (Fig. 2b, c). Among others, the expression of several genes such as Axin2, Ptn, Wif1, Clu and Wfdc18 was particularly characteristic for these tumor-specific cell clusters (Fig. 2c). Using these genes as markers, we further confirmed that the tumor-specific epithelial cells identified in our data were truly located in tumor regions. Immunofluorescence analysis of submandibular gland tissue sections from double-mutant mice showed that antibodies against these gene products positively stained cells within evident tumor lesions (Supplementary Fig. 7). To further increase tumor-specific resolution, we included other epithelial markers where possible. In line with our transcriptome data, Clu and Wfdc18 distinctively stained K8-positive luminal-like cells within the tumor region. Nuclear β-catenin is considered to be the hallmark of active Wnt signaling which ultimately drives the expression of its target genes^35,36 and was previously described to be a key feature of CSCs in several cancers including our model^16,17,37. We, therefore, used this as an additional marker, and found that nuclear β-catenin-positive cells were generally K8 negative and greatly overlapped with high Axin2, but only partially with Ptn and Wif1 stainings. Taken together, we identified a small subset of tumor-specific epithelial cells and provide a marker set that can identify them both at the RNA and protein level.

To validate this trajectory, we made use of the different time points (P40 and P90) at which we collected samples and monitored the contributions of tumor-specific cells pooled by stage (Fig. 5d, e; Supplementary Fig. 12) instead of cell type. In agreement with our model, we found that the number of cells from the early P40 tumor stage was significantly increased at the initial phase of tumorigenesis compared to that from the late P90 one, while the reverse was true at a more advanced phase. Stage-specific quantification of tumor-specific cells showed that although relative proportions of both CSC-like populations were similar, basal tumor and luminal Clu + cells were, in fact, significantly more and significantly less abundant in P40 cells than in P90 cells, respectively (Fig. 5e).

To study potential changes in global expression dynamics along the trajectory, we took the most differentially expressed genes between the four subpopulations and ordered their expression by pseudotime (Fig. 5f, Supplementary Fig 13). The results replicated the differential expression analyses summarized in Fig. 4c, and showed that we could additionally capture more subtle differences for individual genes for which the expression changed as a function of the tumor progression. We found that EMT-related genes such as Ctgf⁵⁵ and Sparc^56,57 were first activated in basal tumor cells, and that their expression extended and decreased within the cell population (Ctgf) itself or throughout both CSC-like subpopulations (Sparc) as tumorigenesis proceeded (Fig. 5f, Supplementary Fig. 13). A similar pattern was observed for Meg3, a long non-coding RNA previously found to regulate EMT⁵⁸ as well as the Wnt/β-catenin and p53 pathways^59,60.

While several genes (e.g., Ptn, S100a6, S100a14) were expressed in a more ubiquitous manner in CSC-like cells, Wnt target genes such as Axin2 and Bmp4 were specifically switched on in the CSC 1 subpopulation and continued to be only mildly expressed and subsequently inactivated in CSC 2 (Supplementary Fig. 13). However, in the latter, increased expression of the transcriptional regulators Nupr1 and Elf3 could be detected, both genes with known functions in epithelial cell differentiation and tumorigenesis^61,62. Finally, Clu, and Wfdc18, which we identified as being specific to luminal Clu+ cells (Fig. 2), were expressed at a later stage towards the end of our inferred trajectory. Validating previously proposed associations between EMT, Wnt signaling and stemness⁶³, we found that gradual loss of the EMT expression signature correlated with activation of Wnt signaling and induction of known CSC marker genes in head and neck squamous cell carcinoma⁶⁴ (Fig. 5c). Moreover, by inspecting differential pathways with mean log2 fold changes below the cutoff imposed in Fig. 4c, we identified dysregulated metabolic signatures suggesting that activation of the EMT program in basal tumor cells is accompanied by metabolic reprogramming from oxidative phosphorylation towards glycolysis (Supplementary Fig. 14). Together, our results show that our approach can identify specific genes and expression patterns potentially regulating and driving tumorigenesis.

In this study, we created a high-resolution salivary gland cell atlas and systematically dissected the cellular heterogeneity in a genetically-controlled, Wnt-dependent mouse model of a solid tumor. We identified, molecularly characterized and validated cell clusters that were specific to the tumor and established their lineage relationship by uncovering the progenitor and progeny populations of CSCs.

Different from other studies involving single-cell transcriptome profiling of human biopsy samples^5,6,7, we could not computationally pinpoint the tumor cells by global copy-number variation (not present in our data), nor by identifying β-catenin and Bmpr1a mutations (not detectable) or EYFP-positive cells (mRNA transcripts not sufficiently captured by our sequencing method). Instead, we leveraged the reproducibility and large sample size of our data to systematically compare cells from control and tumor-bearing tissues to reliably and robustly identify those that were specific to the tumor samples, and confirmed these results by immunostainings in tissue sections.

Although only representing <1% of cells in our data, we were able to detect and extensively characterize CSC-like cells without relying on pre-defined surface markers, which often also target unrelated cells^9,13. We found that Wnt target genes such as Axin2 and Bmp4 were exclusively activated in this small population, but also other genes with highly tumor-specific functions such as Ptn (promoting tumor angiogenesis⁶⁵) and S100a14 (involved in tumorigenesis⁶⁶). Further investigation revealed that these CSC-like cells consist of a Wnt-high/K14+  and a luminal-like/K8+ subpopulation, also characterized by differential expression of the transcriptional regulators Nupr1 and Elf3, which are thought to promote metastasis and the induction of chemoresistance^61,62. Further experiments will be required to functionally characterize these CSC subpopulations and to address their stemness potential. Nevertheless, it has been recognized that different stem cell types with tissue renewal capacity can reside within the same tissue, that a continuum of stem cell states may provide a higher degree of flexibility⁶⁷, and that Wnt signaling is essential for stem and progenitor cell formation and function throughout development^10,68. Our results indicate that the stabilizing β-catenin mutations in our mouse salivary gland SCC model trigger similar processes, as uncontrolled Wnt signaling can lead to aberrant expansion of stem cells or confer stem cell behavior, paving the way for malignant proliferation⁹.

Since the tumor in our model was driven by a K14-cre, we argued that the CSCs most likely originated from the tumor-specific basal population, characterized by an EMT expression signature. The importance of the EMT program is well established as a major mechanism for the ‘invasion-metastasis cascade’ in cancer biology⁶³. With an ability to degrade and reorganize the ECM, it promotes the loss of cell-cell adhesion and the acquisition of migratory and invasive traits. In fact, EMT signatures have been observed in the tumor microenvironment^69,70, but also in malignant cells⁷, and its critical role for the induction of a stem-like phenotype has only been recognized over the last years¹². It has emerged that epigenetic and environmental cues can control EMT without introducing new genetic alterations, that this transition is therefore reversible, rarely fully executed under physiological or malignant conditions, and that an intermediate state between the poles of fully epithelial and mesenchymal cell identity is particularly favorable for the induction of the CSC phenotype^63,71. Specifically, Snai2/Slug in coorporation with Sox9 have been identified as stemness- and EMT-inducing factors in the basal layer of the mammary epithelium⁷². Moreover, the activity of these transcription factors can promote the existence of transient epithelial cell populations which can convert to stem cells with long-term tissue reconstituting ability⁷², and we propose that similar mechanisms play an essential role in our system.

Further, while connections between Wnt signaling both downstream and upstream of the EMT have long been drawn⁶³, we find here that genetic Wnt activation and Bmp deactivation in our model triggers EMT (presumably through Snai2) and metabolic reprogramming, before other Wnt targets, such as Bmp2/4 or Dkk2/4 can be detected.

Our data allowed us to put the different tumor-specific cell populations on a robust trajectory and to map the continuous shifts in cellular identities during tumorigenesis by means of diffusion maps, which have been widely used together with single-cell data to study differentiation dynamics during development^73,74. We note, however, that the directionality of the inferred trajectory is not conclusively fixed and requires further investigation. Nevertheless, based upon the K14-cre-dependence of mutations triggering the tumor in our model as well as the relative timing of early- and late-staged samples, we hypothesized that the tumorigenesis here is initiated in basal cells. This model further suggests that after EMT induction, Wnt signaling was upregulated in one CSC subpopulation in concert with the acquisition of a stemness phenotype that persisted throughout another CSC subpopulation already primed towards a luminal-like cell fate, towards which tumorigenesis finally transitioned. Comparable to epithelial differentiation processes, we noticed a general K14+ to K8+ gradient indicating that tumor-specific cells shifted from exhibiting basal-like to more luminal-like characteristics along the trajectory. This was accompanied by a gradual loss of EMT markers and therefore also suggests a reversal to an epithelial cell identity. Moreover, we found that the luminal-like cells exhibited high expression of Wfdc18, a gene that we identified to be specific to the excretory and intercalated ducts in control tissues. This suggests that in our model, Wfdc18-positive cells in the tumor may have undergone a similar differentiation path as normal cells in these ductal compartments.

Our work demonstrates the importance of using controlled models to robustly and reproducibly study essential mechanisms of carcinogenesis and tumor progression which would otherwise not be possible. Finally, our approach provides a blueprint to molecularly identify markers and characterize transcriptional events that are fundamental in the regulation of tumorigenesis, and ultimately facilitate further clinical studies in the design of appropriate treatment strategies.

Mouse strains

K14Cre(Δneo), β-catenin^flox, Bmpr1a^flox alleles, and Cre-inducible R26^EYFP reporter mice have been described, and mutant mice were analyzed for genotype and recombination by PCR^{16,75,76,77,78,79}. To obtain the double mutants, homozygous mice carrying the β-catenin^flox gain-of-function, the Bmpr1a^flox loss-of-function allele and R26^EYFP were crossed with K14-Cre mice that were homozygous for the Bmpr1a^flox allele. Control mice were obtained by crossing K14-Cre mice with mice carrying the Cre-inducible R26^EYFP reporter. All mice used in this study had a C57BL/6 background. Animal experiments were approved by LAGeSo Berlin and performed according to EU and national institutional regulations.

Tissue dissociation and single-cell sample preparation

Submandibular glands and primary tumor samples were collected, minced and dissociated with a GentleMacs Dissociator (Miltenyi Biotec) in digestion buffer (DMEM/F12 1:1 (Invitrogen), 1.67 mg/ml collagenase (Invitrogen), 1.33 mg/ml hyaluronidase (Sigma) and 1.67 mg/ml dispase (Invitrogen)). Cell suspensions were passed through a stainless filter (70 μm) and centrifuged at 900 × g for 5 min at 4 °C. Pellets were suspended in 10 ml Dulbecco’s modified Eagle/F12 1:1 medium and washed three times with PBS containing 10% fetal bovine serum (Invitrogen). Prior to sorting, cells were stained with DAPI, and then filtered through a 40 μm mesh. The FACS Aria (BD Biosciences) instrument was used for sorting and dead cells excluded by elimination of DAPI-positive cells and gates set to exclude cell clusters. Cells were directly sorted and fixed in ice-cold 80% methanol and stored at −80 °C until further processing.

Drop-seq procedure, single-cell and bulk library generation, and sequencing

Monodisperse droplets of about 1 nl in size were generated using microfluidic PDMS devices (Drop-SEQ chips, FlowJEM, Toronto, Canada; pre-coated with Aquapel). Barcoded microparticles (Barcoded Beads SeqB; ChemGenes Corp., Wilmington, MA, USA) were prepared and flowed in using a self-built Drop-seq set up¹⁹ (Online-Dropseq-Protocol-v.3.1: http://mccarrolllab.com/dropseq/) as previously described¹⁸. Cell preparations and reagents were kept on ice and handled in the cold. Methanol-fixed cells¹⁸ were centrifuged at 3000–5000 × g for 5 min, rehydrated in 1 ml PBS + 0.01% BSA supplemented with RNAse inhibitors (1 unit/μl RiboLock, ThermoFisher), pelleted and resuspended again in 0.5 ml PBS + 0.01% BSA in the presence of RNAse inhibitors. Cells were manually counted by means of a hemocytometer and diluted to a suspension of typically ~50–100 cells/μl in PBS + 0.01% BSA. Droplets were collected in 50 ml Falcon tubes for ~13 min, corresponding to ~1 ml of combined aqueous flow volume (1  ml cells and 1 ml of beads). Droplets were broken immediately after collection and barcoded beads with captured transcriptomes were reverse transcribed and exonuclease-treated. First strand cDNA was amplified by equally distributing beads from one run to 24 PCR reactions (50 μl volume; 4 + 9 to 11 cycles). 20 μl fractions of each PCR reaction were pooled (total = 480 μl), then double-purified with 0.6x volumes of AMPure XP beads (Beckman Coulter). Amplified cDNA libraries were assessed and quantified on a BioAnalyzer High Sensitivity Chip (Agilent) and the Qubit dsDNA HS Assay system (ThermoFisher). If necessary, more cDNA was purified from the PCR reactions. 600 pg of each cDNA library was fragmented, amplified (12 cycles) and indexed for sequencing with the Nextera XT v2 DNA sample preparation kit (Illumina) using custom primers enabling 3’-targeted amplification as described¹⁹. The libraries were double-purified with AMPure XP Beads (0.6 × , 1 × ), quantified and sequenced on Illumina NextSeq500 sequencers (library concentration 1.8 pM; NextSeq 500/550 High Output v2 kit (75 cycles) in paired-end mode; read 1 = 20 bp using the custom primer Read1CustSeqB¹⁹, read 2 = 64 bp).

For bulk sequencing, RNA was extracted with Trizol (Invitrogen) from freshly dissected submandibular glands. RNA integrity was assessed on a BioAnalyzer RNA Nano Chip (Agilent). Strand-specific cDNA libraries were generated from 500 ng total RNA according to the Illumina TruSeq protocol (TruSeq Stranded mRNA LT Sample Prep Kit, Illumina). Libraries were sequenced on an Illumina NextSeq 500 sequencer using the High Output v2 Kit (150 cycles), single read: 150 bp, index read: 6 bp.

CITE-seq experiments

Antibodies were covalently and irreversibly conjugated to DNA-barcoding oligonucleotides by iEDDA click chemistry as previously described⁸⁰. The antibody panel (see Supplementary Table 1) was prepared by mixing equal quantities of each DNA-barcoded antibody and concentrating the panel on an Amicon Ultra 0.5 ml 30 kDa MWCO centrifugal filter (Millipore).

Live intact cells from control and double-mutant submandibular glands of ~P70 mice were prepared as described above. Cells from each animal were processed separately and immediately prepared for simultaneous transcriptome and epitope profiling as outlined in the online CITE-seq protocol (https://cite-seq.com/) using the panel of 63 barcoded antibodies. We further supplemented the cells with a low amount (3%) of human HEK cells (cultured and prepared as described before¹⁸) used as a spike-in control. The cells were then processed for Drop-seq and cDNA libraries prepared as described earlier. Antibody-derived tag (ADT) libraries were amplified for 12 cycles together with TruSeq Small RNA primers for indexing and sequenced together with cDNA libraries on Illumina NextSeq500 sequencers using the same settings as those used previously for single-cell transcriptome libraries.

Immunostainings

Immunofluorescence analyses were performed on formalin-fixed paraffin-embedded tissue sections as described⁷⁶. Antigen retrieval was accomplished by Tris-EDTA (10 mM Tris, 1 mM EDTA, 0.05% Tween-20, pH 9.0) at 99–100 °C for 20 min. Following retrieval, sections were stained with one or several of the following primary antibodies for immunodetection: mouse-anti-β-catenin (BD Transduction Laboratories, 610153), rabbit-anti-EGF (Abcam, ab9695), goat-anti-Smgc (Sigma, SAB2501988), rabbit-anti-Kal1 (aka Wfdc18; Abcam, Ab115270), rabbit-anti-Axin2 (Cell signalling, 2151), rabbit-anti-Wif1 (Abcam, ab186845), rabbit-anti-SMA (aka Acta2; Abcam, ab5694), goat-anti-Prol1 (Abcam, Ab119999), mouse-anti-PGRP (aka Pglyrp1; ThermoFisher, MA1-41044), rabbit-anti-Hepacam2 (Abcam, ab189943), goat-anti-Muc19 (Abcam, ab121014), mouse-anti-AQP5 (Santa Cruz Biotechnology, sc-514022), rabbit-anti-Clusterin (Abcam, ab92548), guinea pig-anti-CK8 (aka K8 or Krt8; Progen, GP-K8), rabbit-anti-CK14 antibodies (aka K14 or Krt14; ThermoFisher, MA5-11599), rabbit-anti-Klk1 (Boster, PA1709), mouse-anti-Ptn (Santa Cruz Biotechnology, SC74443), goat-anti-Sparc (R&D Systems; AF942). Secondary antibodies were conjugated with Cy2, Cy3, or Cy5 fluorochromes (Jackson ImmunoResearch Laboratories). Images were captured using an Axio imager Z1m and AxioCam MRm (Carl Zeiss) and a Leica TCS SP8.

Processing and analysis of single-cell RNA-seq data

Drop-seq data were processed using Drop-seq tools v1.12¹⁹, based on the Gencode vM7 reference augmented by two pseudo-chromosomes containing the Cre and EYFP sequences, respectively. Selecting valid barcodes with the ‘knee’ method, resulting DGEs for each sample were combined and analyzed in R (version 3.4.4) using Seurat²¹ (version 2.3.4). We additionally filtered for at least 100 genes and less than 15% mitochondrial content. Clustering and t-SNE were performed based on a PCA with significant components chosen using JackStraw. For clustering, we used a resolution of 0.7 for the control samples and of 1.0 for the combination of all samples, and merged several acinar subclusters which were not distinguishable from their marker genes. We controlled replicate-to-replicate consistency (see below) across samples, and pooled and analyzed singe-cell datasets without any sample alignment method. For improved visualization, we additionally filtered out all cells located more than 3 standard deviations away from their cluster center and those which had a different identity than the majority of their 10 nearest neighbors in the t-SNE.

Bulk RNA-seq data were mapped to the same augmented Gencode vM7 reference using STAR⁸¹ (version 2.6.0c), quantified using featureCounts (version 1.6.0), and converted to TPM. For the comparison between bulk and single-cell RNA-seq, we used AverageExpression followed by a TPM transformation, and selected all genes found in all samples with TPM > 1.e−4.

For the CITE-seq samples, we used the same processing pipeline on the RNA fraction as for the other Drop-seq samples and filtered for at least 100 but less than 2500 genes and less than 15% mitochondrial content. We regressed out number of UMIs as well as the percentages of mitochondrial and ribosomal genes before performing clustering and t-SNE based on a PCA using 30 significant components. We also performed an additional analysis where we mapped against a combined hg19 and mm10 reference. From the second pipeline run we used reads mapping to the human genome to identify cells of human origin or doublets, and removed these cells from the standard pipeline output. We used CITEseq-count (v1.2) to quantify ADT counts for the respective cells found in both pipeline runs. Specifically binding antibodies were identified by inspecting ADT count distributions for human and mouse cells.

Batch effects and sample-to-sample variability were quantified using an entropy-based approach inspired by the one in Azizi et al.²⁰. Specifically, we aimed to measure how well the local distribution of cells among samples mirrors the global one, taking into account that the latter might not be perfectly uniform. Therefore we used the relative entropy (Kullback-Leibler divergence) \(D_{j}= \sum _{i\epsilon I}q_{i} \,log(q_{i}/q_{i}^{0})\), where q_i is the proportion of cells from group i in the local neighborhood N(j) of cell j, and \(q_i^0\) is the proportion of cells of group i in the entire data set. The set I of sample groups was defined by the biological factors sex, stage, and genotype, by their combination, by technical factors such as cell dissociation protocol, or simply by replicate, and the local neighborhood was defined by taking the 30 nearest neighbors in the kNN-Graph (calculated in PCA space). Controls were obtained by randomly shuffling the group assignment between cells.

The relative density of two sample groups on the t-SNE was plotted using the log2 ratio of two separate 2D kernel density estimators interpolated on the t-SNE coordinates of each cell. Differences in cluster proportions were analyzed using mixed-effects linear models (lme4 package;⁸² version 1.1) using a binomial model with sex, stage and genotype (if applicable) as fixed effects and sample identity as random effect.

Subclustering of cancer stem and basal cells was performed on these populations separately, removing one sample (B9T90R4F) with high content of ribosomal genes and the few cells from control samples from the cancer stem cell clusters. Differential genes between subclusters were detected by pooling counts over all cells from one sample and subcluster, using DESeq2⁸³ (version 1.18) on the pooled counts with the number of cells per group as covariate, and then filtering for genes expressed in at least 5% of cells. Pathway analysis was done with GAGE⁸⁴ (version 2.28) on the log2 fold changes estimated by DESeq2. We used HALLMARK, KEGG and REACTOME gene sets from the Molecular Signature Database (MsigDB: http://www.broad.mit.edu/gsea), kept pathways with a q-value <0.1 and absolute average log2 fold changes >0.5, and ignored those with more than 10% ribosomal genes. For the diffusion map, we combined the tumor-specific basal and CSC subclusters with the luminal Clu+ cluster and used destiny⁸⁵ (version 2.6.2) on the highly variable genes to create the diffusion map embedding and diffusion pseudotime. We also ran destiny on each sample with more than 10 cells in the relevant subpopulations separately and compared sample-specific to global pseudotime estimates. Differential genes along the trajectory were identified on pooled counts using DESeq2 with a likelihood ratio test on all 4 clusters simultaneously, and gene expression values plotted along the trajectory after imputation using SAVER⁸⁶ (version 1.1.1). Gene set scores were computed supplying imputed expression values to Seurat’s AddModuleScore function. An additional “stemness” gene set (Pou5f1, Nanog, Sox2, Prom1, Bmi1, Lgr5, Msi1, Tdgf1, Bmp4, Cspg4, Cxcr4, Alcam, Slc2a13, Aldh1a) was curated from the literature⁶⁴.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

The sequencing data generated is available in GEO under the accession GSE124425. Metadata for each sequenced cell, from which figures based on the single-cell RNA sequencing data can be reproduced, are provided as a source data file.