Differentiation dynamics of mammary epithelial cells revealed by single-cell RNA sequencing

Abstract

Characterising the hierarchy of mammary epithelial cells (MECs) and how they are regulated during adult development is important for understanding how breast cancer arises. Here we report the use of single-cell RNA sequencing to determine the gene expression profile of MECs across four developmental stages; nulliparous, mid gestation, lactation and post involution. Our analysis of 23,184 cells identifies 15 clusters, few of which could be fully characterised by a single marker gene. We argue instead that the epithelial cells-especially in the luminal compartment-should rather be conceptualised as being part of a continuous spectrum of differentiation. Furthermore, our data support the existence of a common luminal progenitor cell giving rise to intermediate, restricted alveolar and hormone-sensing progenitors. This luminal progenitor compartment undergoes transcriptional changes in response to a full pregnancy, lactation and involution. In summary, our results provide a global, unbiased view of adult mammary gland development.

Fig. 1

Single-cell RNA sequencing identifies 15 clusters of mammary epithelial cells. a Schematic diagram highlighting the experimental setup for isolating and sequencing the RNA of single cells using the 10× chromium system. b t-SNE plot of 23,184 cells visualises general structure in the data. Cells are coloured by the four developmental time points as follows: pink = NP, dark green = G, light green = L, purple = PI. c Same as b but colouring cells by clusters. d t-SNEs coloured by the normalised log-transformed expression of the basal marker Krt5 and the luminal marker Krt18

Fig. 2

Putative identities of mammary epithelial cell clusters. a Dendrogram of clusters based on the log-transformed mean expression values of the 15 clusters. The tree was computed based on Spearman’s rank correlation with Ward linkage. b t-SNEs with overlaid expression of cluster-specific genes. c Heatmap highlighting key marker genes that were used to infer putative identities. Colour scale represents log-transformed and normalised counts scaled to a maximum of 1 per row. Upper bars represent the cluster assignment and stages for the individual cells. For visualisation purposes only 100 randomly selected cells were shown for large clusters

Fig. 3

Computational reconstruction of differentiation processes in the mammary gland. a Diffusion map of epithelial cells from the NP and G time points, showing the first three diffusion components. b Differentiation trajectory of the luminal compartment based on the first two diffusion components. c Same plot as in b, coloured by the normalised and scaled expression values of various genes

Fig. 4

Pseudotime ordering identifies genes associated with luminal differentiation. a Definition of the hormone-sensing and secretory differentiation branch. Cells are coloured by their progression through pseudotime, where low values represent undifferentiated cells. b–d Examples of transcription factors with pseudotime-dependent expression with the same overall trend on both branches (b) or branch specific trends (c, d). e, f Heatmap of all genes with branch-specific, pseudotime-dependent expression for the hormone-sensing lineage (e) or the secretory lineage (f). Pseudotime and the cluster assignment are annotated above the heatmap. The values in the heatmap represent z-scaled, spline-smoothed expression values. Genes in the heatmaps were clustered using hierarchical clustering with a dynamic tree cut

Fig. 5

The effect of parity on the transcriptomic landscape of the luminal progenitor compartment. a Comparison of fold changes from C6 vs. luminal compartment and fold changes from C7 vs. the luminal cells. The genes represent the top 500 differentially expressed genes between C6 and luminal cells. b Volcano plot illustrating differential expression between C6 and C7, coloured dots represent significant genes with known function in lactation and immunity, dashed lines highlight the P-value threshold of 0.01 and a log fold change of 1. P-values are adjusted for multiple testing using Benjamini–Hochberg. c, d Visualisation of expression difference for genes related to lactation (c) or immunity (d) for the six luminal clusters in NP and PI. Expression values correspond to normalised UMI counts. e Cells from PI projected on the diffusion map from the NP and G time point. Cells from NP and G are coloured in grey