Spatially distinct epithelial and mesenchymal cell subsets along progressive lineage restriction in the branching embryonic mammary gland

Abstract

How cells coordinate morphogenetic cues and fate specification during development remains a fundamental question in organogenesis. The mammary gland arises from multipotent stem cells (MaSCs), which are progressively replaced by unipotent progenitors by birth. However, the lack of specific markers for early fate specification has prevented the delineation of the features and spatial localization of MaSC-derived lineage-committed progenitors. Here, using single-cell RNA sequencing from E13.5 to birth, we produced an atlas of matched mouse mammary epithelium and mesenchyme and reconstructed the differentiation trajectories of MaSCs toward basal and luminal fate. We show that murine MaSCs exhibit lineage commitment just prior to the first sprouting events of mammary branching morphogenesis at E15.5. We identify early molecular markers for committed and multipotent MaSCs and define their spatial distribution within the developing tissue. Furthermore, we show that the mammary embryonic mesenchyme is composed of two spatially restricted cell populations, and that dermal mesenchyme-produced FGF10 is essential for embryonic mammary branching morphogenesis. Altogether, our data elucidate the spatiotemporal signals underlying lineage specification of multipotent MaSCs, and uncover the signals from mesenchymal cells that guide mammary branching morphogenesis.

Keywords: Branching Morphogenesis; Embryonic Explants; Fate Specification; Mammary Stem Cells; scRNA-seq.

Figure 1. Developmental atlas of the transcriptional signatures and 3D trajectory analysis of luminal and basal differentiation of single mammary epithelial cells from E13.5 until birth.

(A) Scheme showing the isolation and sequencing strategy of mammary embryonic cells at four developmental stages spanning embryonic MG development. (B) UMAP plot of embryonic MECs isolated at E15.5 after subset analysis of non-proliferative MG epithelial cells. Cells are color-coded by cluster. (C) UMAP plots from (B) color-coded according to the expression of the single-cell ID scores in MECs: basal score (left) and luminal score (right). (D) UMAP plot of MECs isolated at P0 after subset analysis of MG epithelial cells. (E) UMAP plots from (D) color-coded according to the expression of luminal progenitors (LP), mature luminal (ML) and basal cell (BC) scores. (F) Violin plots showing the expression levels of the basal and luminal scores in each cluster. n = 22 cells at E13.5; n = 28 cells at E14.5; n = 98 basal-like cells, n = 199 hybrid cells and n = 86 luminal-like cells at E15.5; n = 19 basal cells, n = 140 LP cells and n = 39 ML cells at P0. Statistical significance was assessed with Wilcoxon test. (G) 3D trajectory of MECs from E13.5 at the origin of the mammary cellular hierarchy to P0 MECs positioned at the end of two divergent differentiation routes. Basal and luminal trajectories were calculated using a minimum spanning tree (MST) connecting the center (in the 3D space) of each cell cluster.

Figure 2. Luminal and basal progenitors segregate into different distribution patterns at E15.5.

(A–C) Characteristic examples of the selected genes presenting increased expression in pseudotime towards luminal differentiation (Anxa1, A), basal differentiation (Cxcl14, B) or in the hybrid cells at E15.5 (Cd74, C). Cells are color-coded by cluster as indicated in the legend. (D) Representative sections of embryonic mammary buds at E13.5, E14.5 and E15.5 showing the expression of Cxcl14 (in green) and Anxa1 (in magenta) detected by RNAscope and of K5 by IF (in white). A white dotted line delineates the BM (n = 3). (E) Segmentation of individual mammary cells from (D) that were color-coded according to the ratio between the number of transcripts for Cxcl14 and Anxa1, as indicated on the color key at the bottom. The color key goes from dark purple (100% of transcripts in a cell are Cxcl14) to bright yellow (100% of transcripts in a cell are Anxa1). (F) Quantification of the proportion of Anxa1 and Cxcl14 transcripts in each ring (see Fig. EV3D) at each developmental stage (n = 5). (G, I) Representative sections of embryonic mammary buds at E13.5 and E15.5, showing the expression of Cd74 (in green) and Anxa1 (in magenta) in (G) or Cxcl14 (in magenta) in (I), detected by RNAscope, and of K5 by IF (in white). Dotted lines delineate the BM (n = 2). (H, J) Segmentation of individual mammary cells from (G) and from (I) that were color-coded according to the ratio between the number of transcripts for Cd74 and Anxa1 in (H) or for Cd74 and Cxcl14 in (J). The color key goes from dark purple (100% of transcripts in a cell are Cd74 in (H) or Cxcl14 in (J)) to bright yellow (100% of transcripts in a cell are Anxa1 in (H) or Cd74 in (J)). (K, L) Quantification of the proportion of Cd74 and Anxa1 transcripts in (K) and of Cd74 and Cxcl14 transcripts in (L) in each ring at each developmental stage (n = 2). Data information: graphs in (F, K, L) show mean ± SEM. Statistical significance was assessed with a two-way ANOVA test. ****P < 0.0001, **P < 0.01, ns non-significant. Scale bars: 50 μm in (D, G, I). Source data are available online for this figure.

Figure 3. The early embryonic mammary mesenchyme contains two spatially distinct cell populations.

(A) UMAP plots of embryonic mammary mesenchymal cells isolated at E13.5, E15.5, and P0. Cells are color-coded by cluster. (B) Violin plots representing the expression levels of Esr1, Lef1, Plagl1, Vcan, Cdkn1c, and Crapb1 in sub-epithelial mesenchyme (SE-M) and dermal mesenchyme (D-M) at E15.5 (n = 208 cells). (C, D) Representative sections of embryonic mammary buds at E15.5 immunostained for ERα (in red) and K5 (in white) (C) or LEF1 (in red) and K8 (in white) (C); PLAGL1 (in red) and DAPI (in blue) (D). (E) Representative sections of embryonic mammary buds at E15.5 showing the expression of Vcan (in magenta), Cdkn1c (in green) or Crabp1 (in magenta), detected by RNAscope and immunostained for K8 (in white). Dotted lines delineate the BM (in white) and the two mesenchymal compartments (in orange). Scale bars: 100 μm (n = 3).

Figure 4. Changes in ligand–receptor interaction pairs between both mesenchymal populations and the epithelial cell clusters before and during cell fate switch.

(A) CellPhoneDB analysis with the predicted ligand–receptor interactions between the two mesenchymal populations, sub-epithelial (SE-M) or dermal mesenchyme (D-M), and mammary epithelial cells at E13.5 and E15.5 (P value < 0.01). The arrowhead highlights the ligand–receptor interaction between FGF10 and FGFR2 that was functionally investigated in embryonic ex vivo cultures. Permutation test was used for statistical analysis. (B) UMAP plots from Fig. 3A illustrating the expression of selected Wnt-related cluster-specific genes in mesenchymal cells at E15.5 (left side) and from Fig. 1B illustrating the expression of some cluster-specific Wnt ligands in epithelial cells at E15.5 (right side). (C) Violin plots representing the expression levels of Fgf10 and Fgfr2 in mesenchymal and epithelial clusters, respectively, at E15.5 (n = 383 cells).

Figure 5. FGF signaling is required for embryonic mammary branching.

(A) Time-lapse images of a mammary explant grown in control medium (top) or in the presence of FGF10 (bottom) for 24 h. T = 0 h refers to 4 days in culture. The rendered surface of the mammary epithelium is outlined in blue (in the control bud) and in magenta (in the FGF10 condition). (B) Quantification of the velocity of branch growth in control conditions and in the presence of FGF10 in the medium (n = 43 branches analyzed from 8 independent explants in control conditions and n = 56 branches analyzed from 12 independent explants in FGF10 conditions). (C) Fold change increase in area in control and FGF10 conditions. In both cases, the area is doubled within 16 h in culture. (n = 8 independent explants in control conditions; n = 7 independent explants in FGF10 conditions). (D) Representative whole-mount immunostaining of an embryonic mammary gland cultured in control or FGF10 conditions showing EdU⁺ cells (in green), membrane tdTomato (in white) and DAPI (in blue). Mammary buds were dissected at day E13.5 and cultured ex vivo for 7 days. Orange outlined insets show a duct region and blue outlined insets show a tip region (n = 3). (E–G) Quantification of EdU⁺ cells (E), number of branches (F) and branch diameter (G) in control and FGF10 conditions. n = 6 regions analyzed from at least two independent explants per condition in (E); n = at least seven independent explants per condition in (F); n = 37 branches analyzed from eight independent explants in control conditions and n = 28 branches analyzed from seven independent explants in FGF10 conditions in (G). (H) Representative whole-mount immunostaining of an embryonic mammary gland cultured in control conditions (DMSO) or in the presence of a pan-FGFR inhibitor (BGJ398) starting from day 4 in culture, showing the correct positioning of K8⁺ (in magenta) and P63⁺ (in green) cells even when branching is abrogated by FGFR inhibition. Orange outlined insets show a duct region and blue outlined insets show a tip region (n = 3). (I) Area quantification in control and BGJ398 inhibitor conditions (n = 11 regions analyzed from three independent explants in control conditions and n = 7 regions analyzed from three independent explants in BGJ398 treated conditions). Data Information: scale bars: 100 μm. (A, D, H) Graphs show mean ± SEM. Statistical significance was assessed with two-tailed unpaired Welch’s t test. ****P < 0.0001, **P < 0.01, ns non-significant. Source data are available online for this figure.

Figure 6. Proposed model for lineage segregation of embryonic mammary epithelial cells during development.

(A) Proposed model of luminal and basal differentiation trajectories from E13.5 to P0. (B) Cartoon depicting the spatial distribution of the different cell types distinguishable in the embryonic mammary bud at E13.5, E15.5, and P0.

Figure EV1. Related to Fig. 1. Lineage-committed cells exist in early MG development.

(A) Representative FACS dot plots of the gating strategy used to sort E15.5 and P0 epithelial and mesenchymal cells. (B) UMAP plots of embryonic MECs and surrounding mesenchymal cells isolated by scRNA-seq at E13.5, E14.5, E15.5 and P0. Cells are color-coded by cluster. (C) UMAP plot of embryonic MECs isolated at E15.5 after subset analysis of all MECs (including proliferative cells shown in light blue). (D) Violin plot representation of the cell cycle score in each mammary epithelial cluster at E15.5 (n = 430 cells analyzed). (E) Heatmap showing the expression of genes specific for each cell cluster at E15.5. Each column is color-coded according to the cell cluster from (B). (F) UMAP plots from (C) showing the expression of specific luminal (Krt8 and Krt18) and basal (Krt5 and Trp63) genes commonly used to distinguish adult LCs and BCs but unable to discriminate distinct cell clusters at E15.5. (G) Box plots illustrating the log2 fold change of the luminal/basal score ratio in each cluster. n = 22 cells at E13.5; n = 28 cells at E14.5; n = 98 basal-like cells, n = 199 hybrid cells and n = 86 luminal-like cells at E15.5; n = 19 basal cells, n = 140 LP cells and n = 39 ML cells at P0. Statistical significance was assessed with Wilcoxon test. Lower and upper hinges correspond to the first and third quartiles. The upper whisker extends from the hinge to the largest value no further than 1.5 * IQR from the hinge (where IQR is the inter-quartile range, or distance between the first and third quartiles). The lower whisker extends from the hinge to the smallest value at most 1.5 * IQR of the hinge. Data beyond the end of the whiskers are called “outlying” points and are plotted individually.

Figure EV2. Related to Fig. 1. Pseudotime ordering identifies genes associated with early luminal and basal differentiation and changes in ligand-receptor interaction pairs between basal-like and luminal-like at early and late developmental timepoints.

(A, B) Heatmaps illustrating genes exhibiting a differential pattern of expression along the pseudotime (from E13.5 to P0) towards the basal lineage (A) or the luminal lineage (B). Genes (rows) are clustered based on the dendrogram on the left and color-coded by their expression levels (from blue to red). The gene expression levels were smoothed using the generalized additive model (GAM) and scaled by row. Genes of interest are indicated on the right. Each set of genes with a specific pattern is color-coded on the left: 5 distinct patterns in the basal lineage (A) and 7 unique patterns in the luminal lineage (B). (C) CellPhoneDB analysis showing predicted ligand-receptor interactions between the two epithelial populations at E15.5, basal-like and luminal-like cells (left side of the dotplot), and between the three epithelial populations at P0; BCs, LP and ML cells (right side of the dotplot) (P value < 0.01). Permutation test was used for statistical analysis.

Figure EV3. Related to Fig. 2. Identification of novel genes that distinguish lineage-biased embryonic mammary cells.

(A) UMAP plots showing the expression levels of selected basal (Cxcl14, Pthlh and Ndnf), hybrid (Cd74) and luminal (Lgals3, Anxa1 and Plet1) genes at E15.5. (B) Heatmap illustrating the expression of genes specific for each MEC cluster at E15.5. Each column is color-coded according to the cell cluster from Fig. 1B. Black arrowheads indicate genes used in RNAscope experiments. (C) Representative section of a mammary duct at P0 showing the expression of Cxcl14 (in green) and Anxa1 (in magenta) detected by RNAscope and of K8 by IF (in white). The white dotted line delineates the BM (n = 3). (D) Optical section of a mammary bud at E15.5 illustrating the ROIs: outer ring (in green), middle ring (in blue) and internal ring (in magenta) used for the quantitative analysis. (E, F) Representative sections of embryonic mammary buds at E15.5 showing the expression of Ndnf (basal gene, in green) and Plet1 (luminal gene, in magenta) (E) or Pthlh (basal gene, in green) and Lgals3 (luminal gene, in magenta) (F), detected by RNAscope and immunostained with antibodies anti-K5 (in white) (n = 2). (G) Single optical section showing the expression of the luminal epithelial marker K8 (in magenta), and the basal epithelial marker P63 (in green) in an embryonic mammary bud at E15.5. K8 and P63 are co-expressed by all MECs at E15.5. (H) Representative sections of the 5 different embryonic mammary buds (#1 to #5) at E15.5, showing the conserved expression of Anxa1 (in magenta) and Cxcl14 (in green) detected by RNAscope. Anti-K5 or anti-K8 immunostaining delineate the mammary bud epithelium (in white) (n = 2). Data information: scale bars: 100 μm in (C) and 50 μm in (E–H).

Figure EV4. Related to Fig. 3. The heterogeneity of mesenchymal cells increases at birth.

(A) Violin plots representing the cell cycle score in each mammary mesenchymal cluster at E13.5 and E15.5. n = 129 cells at E13.5; n = 252 cells at E15.5. (B) Representative sections of mammary bud at E15.5 and E16.5 and whole-mount staining at E17.5 and P0 showing PH3⁺ cells (in green), K8 (in white) and DAPI (in blue) (n = 2). Dotted lines delineate the BM (in white). (C) Heatmap illustrating the expression of genes specific for each mesenchymal cluster at E15.5. Each column is color-coded according to the cell cluster from Fig. 3A. The black arrowheads indicate the genes that were further investigated for their specific expression in sub-epithelial or dermal mesenchyme. (D) UMAP plots from Fig. 3A illustrating the expression of cluster-specific genes in mesenchymal cells at P0. Scale bars: 100 μm in (B).

Figure EV5. Related to Fig. 5. Mammary bud ex vivo cultures recapitulate embryonic mammary morphogenesis and epithelial lineage segregation.

(A–C) Representative images of mammary embryonic buds dissected at day E13.5 and cultured ex vivo for 1 day, immunostained for the following lineage markers: P63 (in green) and K14 (in magenta) (A), P63 (in green) and K8 (in magenta) (B), and K5 (in green) (C). (D–F) Representative images of mammary embryonic buds dissected at day E13.5 and cultured ex vivo for 8 days, immunostained for the following lineage and polarity markers: P63 (in green) and K8 (in magenta) (D), α-SMA (in green) and K8 (in magenta) (E), and P63 (in green) and ZO-1 (in magenta) (F). Data information: scale bars: 50 μm (in A–C), 100 μm (in D–F) (n = 3).