High Expression of CD200 and CD200R1 Distinguishes Stem and Progenitor Cell Populations within Mammary Repopulating Units.

Aiming to unravel the top of the mammary epithelial cell hierarchy, a subset of the CD49fhighCD24med mammary repopulating units (MRUs) was identified by flow cytometry, expressing high levels of CD200 and its receptor CD200R1. These MRUCD200/CD200R1 repopulated a larger area of de-epithelized mammary fat pads than the rest of the MRUs, termed MRUnot CD200/CD200R1. MRUCD200/CD200R1 maintained a much lower number of divergently defined, highly expressed genes and pathways that support better cell growth, development, differentiation, and progenitor activity than their MRUnot CD200/CD200R1 counterparts. A defined profile of hierarchically associated genes supporting a single-lineage hypothesis was confirmed by in vitro mammosphere analysis that assembled 114 genes with decreased expression from MRUCD200/CD200R1 via MRUnot CD200/CD200R1 toward CD200+CD200R1- and CD200R1+CD200- cells. About 40% of these genes were shared by a previously published database of upregulated genes in mammary/breast stem cells and may represent the core genes involved in mammary stemness.

Introduction

The mammary gland is a highly regenerative organ that experiences periods of development, function, and regression dictated by embryo and newborn development. The mammary gland develops as a rudimentary ductal network headed by terminal end buds in mice, or terminal ductal lobuloalveolar units in humans and bovines, which branch upon pregnancy (Cardiff and Wellings, 1999). Hormonally stimulated lobuloalveolar structures differentiate into milk-producing alveoli during lactation and involute upon neonate weaning. These changes mainly involve the epithelial layers. Indeed, the mammary gland is an epithelial organ. Luminal epithelium lines the apical part of the ducts and encompasses the milk-producing cells upon lactation. An outer contractile basal myoepithelium forces the milk toward the nipple (Hennighausen and Robinson, 2005).

A Lin−, CD24medCD49fhigh/CD24+CD29high cell population was identified as the entity that enables mammary gland development and regeneration (Shackleton et al., 2006, Stingl et al., 2006). This population is enriched in mammary repopulating units (MRUs) with stem-like properties that can generate an entire functional gland upon transplantation into de-epithelized mammary fat pad. Further studies reported multiple characteristics of a relatively heterogeneous entity: MRUs with relatively lower repopulation ability but rapid growth were identified during pregnancy (Asselin-Labat et al., 2010). Possibly different, actively proliferating s-SHIP-expressing stem cells with relatively high self-renewal ability were identified at the forefront of the virgin's developing terminal end buds and in the alveolar buds during pregnancy (Bai and Rohrschneider, 2010). In contrast, slow-cycling H2b-GFP+ mammary stem cells exhibiting Cd1d expression were also identified (dos Santos et al., 2013).

An additional layer of complexity at the top of the mammary cell hierarchy was added by the identification of distinct long-lived progenitors (Rios et al., 2014), which may resemble the multipotent progenitors that descend from the hematopoietic stem cells (Seita and Weissman, 2010). Further characterization of this defined population is needed to understand the initial steps of stem cell differentiation.

CD200 (OX-2), a transmembrane surface glycoprotein, is a member of the immunoglobulin superfamily. It is widely distributed across tissues, including lymphocytes, endothelial cells, keratinocytes, and neuronal cells (Ko et al., 2009). Its cognate receptor, CD200R1, is also an immunoglobulin transmembrane glycoprotein and is mainly expressed on myeloid cells. Binding of CD200 to its receptor attenuates immune activity (Minas and Liversidge, 2006). No information regarding CD200's role in the normal mammary gland is currently available. However, its expression has been detected in breast cancer (Leccia et al., 2012).

CD200 expression has been associated with stemness. In the presence of low expression levels of CD24, CD34, CD71, and CD146, high CD200 expression enabled successive enrichment of stem cells in human hair follicle bulge cells (Ohyama et al., 2006). Furthermore, high CD200 expression distinguished fetal- from maternal-originated placental mesenchymal stem cells (Zhu et al., 2014). In the mammary gland, however, high CD200 expression levels do not appear to mark stem cells (G.R., unpublished).

Seeking markers of mammary epithelial stem cells, a CD200highCD200R1high epithelial cell population was identified in the current study. Members of this population enriched the MRUs for cells with better repopulating activity compared with the rest of the MRUs. Gene-expression profiling combined with comparative dataset analysis, supported by further in vivo and in vitro studies, suggested the coexistence of two MRU subsets: a poorly metabolically active stem cell population with relatively high complement activity, and highly active differentiation-oriented progenitors.

Results

An MRU Subpopulation Expressing High Levels of CD200 and CD200R1 Exhibits Enhanced Repopulation Ability

A mouse mammary cell population expressing high levels of CD200 and CD200R1 was identified by flow cytometry (Figure 1A). This population included 82.8% ± 16.6% (n = 3) epithelial cells and represented 3.3% ± 0.8% of the mammary epithelial cells. Projecting the CD200highCD200R1high cells on the CD24/CD49f expression axes located 49.2% ± 18.7% of the cells within the CD24medCD49fhigh (MRU) boundaries (Stingl et al., 2006), representing 50.1% ± 11.9% (n = 3) of the MRUs (Figure 1B). The MRUs that expressed high CD200 and CD200R1 levels are termed here MRUCD200/CD200R1. To examine their repopulating potential, outgrowth development from these cells was compared with that developing from the rest of the MRUs, termed MRUnot CD200/CD200R1. As shown in Figure 1C, no difference was observed between the repopulating potential of the two subpopulations, and 40%–50% of the fat pads transplanted with 40 cells from each MRU subset were occupied by newly developed epithelium. Further analyses identified fat pads that were completely filled with outgrowths, whereas others were only partially occupied (Figures 1D and 1E, respectively). Transplantation of limiting numbers of MRUCD200/CD200R1 and MRUnot CD200/CD200R1 into cleared mammary fat pads revealed a 2.6-fold decrease in full repopulation frequency for the MRUnot CD200/CD200R1 (versus MRUCD200/CD200R1), which tended toward significance (p = 0.06, Figure 1F).

Figure 1

MRUs that Express High Levels of CD200 and CD200R1 Exhibit Better Repopulation Ability Than the Other MRUs

(A and B) Dot plots depicting the gating strategy for mouse mammary cells sorted for transplantation. Cells were analyzed simultaneously for CD200, CD200R1, CD24, and CD49f expression. (A) Identification of CD200highCD200R1high epithelial cell population. (B) Projection of the CD200highCD200R1high population on the CD49/CD24 axis identified two MRU (CD24medCD49fhigh) subpopulations: MRUCD200/CD200R1 that expresses high levels of both CD200 and CD200R1 (red) and MRUnot CD200/CD200R1 that represents the rest of the MRUs (green).

(C) Limiting dilution analysis of the repopulating frequency of MRUCD200/CD200R1 and MRUnot CD200/CD200R1 cells from 8-week-old virgin mice. CI, confidence interval. ∗Number of outgrowths per number of injected fat pads.

(D and E) Whole-mount Carmine alum staining of transplanted mouse fat pads depicting fully (D) and partially (E) reconstituted glands. Bar, 1 mm.

(F) Limiting dilution analysis of the potential for full fat pad occupancy by cells of the two MRU subpopulations. ∗Number of outgrowths per number of injected fat pads. For details, see (C).

(G and H) Analysis of the absolute and relative areas occupied by new epithelium after transplantation of the indicated number of cells from the two MRU subpopulations. Bars represent mean ± SEM of three replications. Asterisks mark statistically significant difference (p < 0.05).

See also Tables S1 and S5.

The absolute and relative areas of the reconstituted glands were also larger by 30% and 33%, respectively, for glands transplanted with MRUCD200/CD200R1 compared with those transplanted with MRUnot CD200/CD200R1 (Figures 1G and 1H). Significant differences (p ≤ 0.05) between the two subpopulations were obtained for transplantation of 100 cells. Combining data from transplanting 100 and 250 cells for each of the two MRU subpopulations (Figure 1G) also resulted in significantly (p ≤ 0.05) higher occupancy rates for the MRUCD200/CD200R1 developing cells (not shown). It is possible that transplantation of 100 cells per fat pad provided the most permissive environment, necessary for executing the outgrowth's full developmental potential.

In a serial transplantation analysis, mammary epithelial cells prepared from pooled outgrowths developed from MRUCD200/CD200R1 occupied 11% (2/18) of the fat pads transplanted. In contrast, MRUnot CD200/CD200R1 did not give rise to new epithelium (Table S1).

Differential Gene Expression Allocates a Larger Number of Highly Expressed Genes and Activated Pathways to the MRUnot CD200/CD200R1 Subpopulation

The observed difference in fat pad occupancy upon transplantation of the two MRU subpopulations reflected diversity within the MRUs. Therefore, the two MRU-composing populations, as well as the non-MRU CD200+ population (MRU−CD200R1−) and CD200R1+ population (MRU−CD200−), were sorted. Sets of 1,000 cells of each population were collected from five individual mice, lysed, and subjected to RNA sequencing (RNA-seq). On average, ∼21 million reads per sample were obtained. Of these, ∼95% were mapped to the mouse genome and ∼57% were found on known genes. Gene-expression analysis clearly clustered the samples into four cell populations, independent of the mouse effect (Figure 2A). To further distinguish the two MRU subpopulations, the expression level of each gene in MRUCD200/CD200R1 cells was related to its expression in the MRUnot CD200/CD200R1 cells in individual glands, and log2 values were calculated. Values of >1 indicated higher expression in MRUCD200/CD200R1 and <−1 marked higher expression in their counterpart. As shown in Figure 2B, a larger number of highly expressed genes was detected in the MRUnot CD200/CD200R1. The difference between this population and the MRUCD200/CD200R1 population ranged from 2.3-fold with significant differences of p ≤ 0.05 up to 10.3-fold for adjusted p ≤ 0.01. Within a defined list of annotated genes with 1.5-fold difference between the two populations and adjusted p < 0.07, the six genes with the highest relative expression in the MRUCD200/CD200R1 were CD5l, Siglec1, Kcnn1, Fcna, C6, and Ccl8 (Figure 2C). These genes were expressed at 5.5- to 8.6-fold higher levels than in the MRUnot CD200/CD200R1 cells, and were all associated with immune system regulation (Asano et al., 2015, Begenisich et al., 2004, Endo et al., 2012, Maggiani et al., 2011, Martinez-Picado et al., 2016, Mayilyan, 2012). Higher differences, 64- to 194-fold their counterparts' expression, characterized genes that were highly expressed in the MRUnot CD200/CD200R1; none of these genes had a direct association with immune regulation. This list included Igf2, coding a protein hormone known to regulate cell proliferation, growth, migration, differentiation, and survival (Bergman et al., 2013), but whose effect on mammary gland development has not been elucidated.

Figure 2

Unsupervised Hierarchical Clustering of Distinct Epithelial Cell Populations in This Study

The MRUnot CD200/CD200R1 subpopulation maintains a higher number of highly expressed genes and metabolic pathways compared with its MRUCD200/CD200R1 counterpart.

(A) Heatmap demonstrates clear separation of the MRUCD200/CD200R1, MRUnot CD200/CD200R1, CD200+CD200R1− and CD200−CD200R1+ populations from mammary glands of five individual mice.

(B) Number of genes that are highly expressed in each of the compared subpopulations according to the defined criteria.

(C) The most highly expressed genes in each of the two MRU subpopulations.

(D) Number of pathways, determined by IPA, that are highly activated in each of the compared populations according to the defined criteria.

(E) Highly activated pathways in each of the compared populations.

See also Figures S1 and S2.

Pathway activity supported by defined genes was determined using Ingenuity Pathway Analysis (IPA) software. The parameters defined for selected differences in gene expression between the two MRU populations, >1.5-fold difference, adjusted p < 0.05, provided a list of 1,278 genes. Importantly, only a few pathways were highly activated in the MRUCD200/CD200R1 compared with MRUnot CD200/CD200R1 cells (Figures 2D and S1). Among these, the highest Z-scores were associated with genes regulating the complement and coagulation systems. Many more pathways were highly activated in the MRUnot CD200/CD200R1 population. Higher Z-scores were calculated for nitric oxide signaling, interleukin (IL)-8 signaling, and endothelial nitric oxide synthase (eNOS) (Figure 2E). Pathway activities that characterized MRUCD200/CD200R1 cells involved inflammatory response, organismal injury, cell death, and cancer. In contrast, major activity in the MRUnot CD200/CD200R1 cells involved cell movement, tissue development, and cell development, growth, and differentiation (Figure S2).

Table 1 assembles the main molecular regulators and activators of the IPA-defined pathways. It shows that the molecular basis of the highly activated pathways in the MRUnot CD200/CD200R1 population consisted of an induced core of genes encoding multipotent kinases, GTPases, and G-protein subunits that support major cellular activities: proliferation, differentiation, and migration. Fewer encoded proteins: SRC, protein kinase C (PKC), and SHIP, maintained the opposite activity in these pathways and may serve as negative regulators. A much more modest repertoire of activated pathways characterized the MRUCD200/CD200R1 population. High activity of C proteins supported higher activity of the complement system. A chaperone of coagulation factor VII, von Willebrand factor (vWF), IL-6, and MAGI negatively regulate the coagulation system, and TREM1 and pTEN signaling, respectively, which were also upregulated.

Table 1

Key Encoded Proteins that Are Highly Involved in Mediating Pathway Activity in MRUCD200/CD200R1 and MRUnot CD200/CD200R1 Populations

Key Proteins	Pathways
Key proteins highly expressed in MRUCD200/CD200R1 populationa	highly active pathways in MRUCD200/CD200R1 populationa
Complement (C1q, C3, C3a, C3b, C4b, C6)	complement pathway
Key proteins highly expressed in MRUnot CD200/CD200R1 populationb	highly active pathways in MRUCD200/CD200R1 populationa
vWF	coagulation system
IL6	TREM1 signaling
MAGI	pTEN signaling
Key proteins highly expressed in MRUnot CD200/CD200R1 populationb	highly active pathways in MRUnot CD200/CD200R1 populationb
PI3 kinase	eNOS signaling; nitric oxide signaling; IL-8 signaling; Tec signaling; Gαq signaling; inhibition of angiogenesis by TSP1; endothelin-1 singling; anti-proliferative role of somatostatin receptor 2; thrombin signaling; HMGB1 signaling; p70s6k signaling; integrin signaling; hepatocyte growth factor (HGF) signaling
ERK 1/2	nitric oxide signaling; IL-8 signaling; antiproliferative role of somatostatin receptor 2; vascular endothelial growth factor (VEGF) family of ligand receptor interactions; HMGB1 signaling; prolactin signaling; agrin interactions; Stat3 signaling; nuclear factor (NF)-κB signaling; thrombin signaling; acute-phase response signaling; CAMP-mediated signaling; chemokine signaling; Rho signaling; role of NFAT in cardiac hypertrophy; renin-angiotensin signaling; HGF signaling; mouse embryonic stem cell pluripotency
p38 MAPK	inhibition of angiogenesis by TSP1; thrombin signaling; acute-phase response; type 1 diabetes; HMGB1 signaling; iNOS signaling; chemokine signaling; IL-1 signaling; VEGF signaling; TGFβ signaling; mouse embryonic stem cell pluripotency
PKC	nitric oxide signaling; IL-8 signaling; eNOS signaling; Inhibition of angiogenesis by TSP1; p70s6k signaling
Gα	Tec kinase signaling; endothelin 1 signaling; signaling by Rho family GTPases; IL-1 signaling
Gβ	Tec kinase signaling; IL-8 signaling; Gαq signaling; cardiac β-adrenergic signaling; antiproliferative role of somatostatin receptor 2; thrombin signaling; role of NA in cardiac hypertrophy; signaling by Rho family GTPases; α-adrenergic signaling; IL-1 signaling
GƳ	Tec kinase signaling; IL-8 signaling; Gαq signaling; cardiac β-adrenergic signaling; antiproliferative role of somatostatin receptor 2; thrombin signaling; role of NA in cardiac hypertrophy; signaling by Rho family GTPases; α-adrenergic signaling; IL-1 signaling
RAS	phospholipase signaling; thrombin signaling; acute-phase response; NF-κB signaling; integrin signaling; antiproliferative role of somatostatin receptor 2; thrombin signaling; glioma invasiveness phospholipase signaling; thrombin signaling; acute-phase response; NF-κB signaling; integrin signaling; HMGB1 signaling; agrin interactions; P70S6K signaling; α-adrenergic signaling; NRF-mediated stress response; prolactin signaling; p21 activated kinase (PAK) signaling
IκB	Gαq signaling; acute-phase response signaling; type 1 diabetes; iNOS signaling
JNK	IL-8 signaling; Tec signaling; inhibition of angiogenesis by TSP1; HMGB1 signaling; type 1 diabetes; Rho signaling; IL-1 signaling; PAK signaling; HGF signaling
FAK	glioma invasiveness; IL-8 signaling; Tec signaling; agrin interactions; Rho signaling; chemokine signaling; integrin signaling
RHO	glioma invasiveness; HMGB1 signaling; IL-8 signaling; Gαq signaling
VEGF and receptors 1/2/3/C	nitric oxide signaling; IL-8 signaling; VEGF signaling
CAVEOLIN-1	nitric oxide signaling; eNOS signaling; integrin signaling
eNOS	eNOS signaling; nitric oxide signaling; VEGF signaling
P21CIP1	HGF pathway
Key proteins highly expressed in MRU CD200/CD200R1 populationa	highly active pathways in MRUnot CD200/CD200R1 populationb
SRC	IL-8 pathway; antiproliferative role of somatostatin receptor 2; thrombin pathway; agrin interactions; Stat3 pathway; integrin pathway; cAMP-mediated signaling; endothelin 1 signaling; phospholipase C signaling; role of NA in cardiac hypertrophy; p70s6k signaling; integrin signaling; VEGF signaling; Gα12-13 signaling
BTK	Gαq signaling; phospholipase C signaling; NF-κB/PI3K signaling; p70s6k signaling; Gα12-13 signaling
PKC	VEGF signaling
SHIP	NF-κB/PI3K signaling

See also Figures S1 and S2.

Divergent pathway activity in the MRU subpopulations was determined by IPA.

Corresponds to MRU separation presented in Figure 1: MRUCD200/CD200R1 expresses high levels of both CD200 and CD200R1.

Corresponds to MRU separation presented in Figure 1: MRUnot CD200/CD200R1 represents the rest of the MRUs.

Expression of CD200 and CD200R1 Marks Distinct MRU Subpopulations with Stem and Progenitor Characteristics

The higher occupancy rate of the mammary fat pad by cells of the MRUCD200/CD200R1 population (Figure 1) suggested two possible hierarchical interactions with the MRUnot CD200/CD200R1 population. The latter may represent progenitors of the MRUCD200/CD200R1 population (Figure 3A) or an independent subgroup of stem cells (Figure 3B). The much higher metabolic activity and, in particular, expression levels of pro-differentiation molecules such as p38 mitogen-activated protein kinase (MAPK) and ERK1/2 in the MRUnot CD200/CD200R1 population supported a single hierarchy, initiated by stem cells via immediate progenitors toward differentiated CD200- and CD200R1-expressing cells.

Figure 3

Co-occupancy of the MRU Fraction by MRUCD200/CD200R1 and MRUnot CD200/CD200R1 May Involve Hierarchal Association or the Presence of Two Stem Cell Populations

(A and B) Demonstration of the two options.

(C) Genes with expression patterns that fit with stem cells being in a single hierarchy: expression of genes in MRUCD200/CD200R1 is greater than in MRUnot CD200/CD200R1 at p < 0.05; expression of genes in MRUnot CD200/CD200R1 is 4-fold greater than in CD200+CD200R1− and CD200R1+CD200−.

(D) Gene-expression pattern indicating progenitors: expression of genes in MRUnot CD200/CD200R1 is 2-fold greater than in MRUCD200/CD200R1, CD200+CD200R1− and CD200R1+CD200−.

(E) Gene-expression pattern indicating differentiated cells: expression of genes in CD200+CD200R1− and CD200R1+CD200− is at least 2-fold greater than in MRUnot CD200/CD200R1 with adjusted p ≤ 0.05; expression of genes in MRUnot CD200/CD200R1 is greater that in MRU CD200/CD200R1 at adjusted p ≤ 0.75.

See also Table S2.

To confirm this hypothesis, a gene subset was arrayed, with an expression profile that fit the pattern of a single lineage consisting of decreased expression from the MRUCD200/CD200R1 cells via the MRUnot CD200/CD200R1 cells toward mammary epithelial cells expressing either CD200 or CD200R1 (Figure 3C). Initial gene ontology analysis of the 114 identified genes (Table S2) showed a strong association between 19 of them and epithelium. The contribution of other tissues—epidermis, muscle, and mesenchyme—also produced significant false discovery rate values (0.001–0.0.4) (data not shown). Next, the 114 genes were analyzed by STRING (Jensen et al., 2009) to elucidate the interactions among them. Two complexes and a connecting branch were identified (Figure 4). Complex A included nine extracellular matrix-encoded proteins, mainly fibril-forming, basement-forming, FACIT and transmembrane collagens (Gelse et al., 2003). Complex A interacted with complex B via LAMA1 and its heterodimeric receptor ITGA6 that binds ACTN1. Complex B contained nine genes involved in muscle (mainly smooth muscle) architecture and function (Chereau et al., 2008, Lazar and Garcia, 1999, Rensen et al., 2000). Genes coding for proteins composing the connecting branch represented mesenchyme (Ctgf) and basal epithelium (Cdh3, Trp63 and Krt5). The branch also included FCNA, which is highly expressed in the MRUCD200/CD200R1 cells and binds proteins of the complement system (Endo et al., 2012), as well as the regulatory protein WNT10A. The latter is induced in early embryogenesis and enhances the self-renewing phenotype in cancer cells (Long et al., 2015). Also present were the WNT receptor FZD7 that mediates WNT/β-catenin signaling for pluripotency (Fernandez et al., 2014), and HUNK, which is involved in the formation and differentiation of adult intestinal stem cells (Luu et al., 2013).

Figure 4

Analysis of the Stem Cell Subset for Interactions between the Encoded Proteins Demonstrates Two Complexes and a Linking Branch

The genes presented were analyzed by STRING, which determines known and predicted protein–protein interactions. Bold red lines mark interactions between proteins encoded by genes that appeared in the previously published list of stem cell markers (Lim et al., 2010). Bold dashed line marks interaction with markers of stem cells in other tissues.

See also Table S3.

To test the relevance of the selected list of 114 genes to mammary stemness, it was compared with an external dataset defining genes that are upregulated in an enriched mammary stem cell gene subset in both mice and humans (Lim et al., 2010, and Table 2, row 1). An overlap of 39 genes allowed a relatively highly significant p value. Seven additional genes, Ltbp2, Hal, Wnt10A, Bst1, Fzd7, Lgr4, and Trnp1, which are upregulated in enriched stem cell subsets of other tissues or in embryonic stem cells, were also identified (Table S3). Most of the common genes were assembled into the aforedescribed complexes and connecting branch (Figure 4). Cdh3, which links the two complexes via the branch, was one of the highly expressed genes in MRUCD200/CD200R1 cells compared with the rest of the MRU-defined cells (1.77-fold). Interestingly, two of the non-mammary-associated genes, Lgr4 and Lgr6 that mark stemness in the skin and intestine, respectively (Li et al., 2015, Lough et al., 2013), were not connected to the complexes. Fzd7 and Hall, which were also not associated with the complexes, were the most highly expressed genes in MRUCD200/CD200R1 versus MRUnot CD200/CD200R1 cells (1.97-fold and 4.89-fold, respectively).

Table 2

Comparative Analysis of Genes that Mark Cell Status in the Hierarchy in This and Other Studies

No.	Aim	Source of Gene Set 1	Number of Genes in Gene Set 1	Source of Gene Set 2	Number of Genes in Gene Set 2	Number of Overlapping Genes	Significance (in Set of 20,000 Genes)
1	validating proposed cell hierarchy	current analysis: stem cell subset	114a	enriched mammary stem cell subset (Lim et al., 2010)	489	39b	p = 2.75 × 10−34
2		current analysis: stem cell subset	114	transcriptional regulators of mammary stem cells (Chakrabarti et al., 2014)	26	6	p = 6.3 × 10−9
3	validating involvement of luminal progenitors	current analysis: progenitor subset	262	enriched mammary luminal progenitor subset (Lim et al., 2010)	58	0	not significant
4	validating involvement of differentiated luminal cells	current analysis: differentiated cell subset	97	enriched mammary luminal differentiated cell subset (Lim et al., 2010)	116	0	not significant
5	comparing with gene expression in the breast	current analysis: stem cell subset	114	upregulated genes in basal mammary epithelial cells compared with luminal ones (Huper and Marks, 2007)	54	4	p = 0.0004
6		current analysis: progenitor subset	262	upregulated genes in basal mammary epithelial cells compared with luminal ones (Huper and Marks, 2007)	54	0	not significant
7		current analysis: differentiated cell subset	97	upregulated genes in basal mammary epithelial cells compared with luminal ones (Huper and Marks, 2007)	54	1	p = 0.2
8	comparing with gene expression of breast cancer subtypes	current analysis: stem cell subset	114	upregulated genes in basal breast cancer samples (Smid et al., 2008)	777	14	p = 0.0001
9		current analysis: progenitor subset	262	upregulated genes in basal breast cancer samples (Smid et al., 2008)	777	10	p = 0.12
10		current analysis: differentiated cell subset	97	upregulated genes in basal breast cancer samples (Smid et al., 2008)	777	4	p = 0.19

See also Tables S3 and S4.

105 annotated.

Additional stem cell-related genes were identified.

An additional comparison with a smaller dataset of 26 mammary stem cell transcriptional regulators (Table 2, row 2) identified six genes in common with the putative stem cell list of 114 genes, thus supporting their association with stemness.

The significant number of common genes shared by the external stem cell datasets and the list of genes with decreased expression from MRUCD200/CD200R1 via MRUnot CD200/CD200R1 toward CD200- and CD200R1-expressing cells supports the concept of a single-cell hierarchy that assembles the MRUCD200/CD200R1 and the MRUnot CD200/CD200R1 populations as stem cells and their progenitors. A related question was further raised: does the downstream cell hierarchy toward the CD200+ and CD200R1+ cells encompass the complete epithelial cell population?

To answer this, an additional set of 262 genes was identified with highest expression in MRUnot CD200/CD200R1, aimed to mark mammary progenitors (Figure 3D). Indeed, this set of genes had 28 orders of magnitude lower compatibility to the stem cell subset (Lim et al., 2010) compared with the putative stem cell list defined here (Figure 3C and Table S4). Interestingly, this progenitor set did not share any common genes with a signature of 58 mouse/human luminal progenitors (Lim et al., 2010 and Table 2, row 3). The lack of luminal orientation was confirmed by the absence of common genes between an additional subset of genes highly expressed in the CD200+ or CD200R1+ differentiated cells (Figure 3E) and a list of genes upregulated in differentiated luminal cells (Lim et al., 2010 and Table 2, row 4). As expected, this subset had 30 orders of magnitude lower compatibility with the reported stem cell subset (Lim et al., 2010) as compared with the stem cell list described here (Table S4).

The latter results indicated that the proposed cell hierarchy stemming from MRUCD200/CD200R1 does not include a defined luminal lineage. To test for basal orientation, a comparison was performed between the stem, progenitor, or differentiated subsets defined in the current study and gene sets that favor basal versus luminal expression in the mammary gland (Table 2, rows 5–7) or breast cancer (Table 2, rows 8–10). This comparison only yielded a significant number of common genes for the stem cell subset defined here (Table 2, rows 5 and 8). This implies that the cell hierarchy originating from the MRUCD200/CD200R1 cells is indeed topped by a highly enriched stem cell population with basal characteristics. However, the downstream hierarchy toward the CD200- and CD200R1-enriched populations involves a relatively narrow lineage that does not include specific basal versus luminal characteristics. This lineage may, however, include specific immune-associated characteristics of the CD200 and CD200R1 cells. Indeed, within 12 relevant immune-related databases that might shed light on the downstream characteristics of cell hierarchy (Subramanian et al., 2005), the one that provided the most significant number of common genes (p = 0.0058) related to upregulated genes in normal hematopoietic progenitors by RUNX1–RUNX1T1 translocation (Tonks et al., 2007). This translocation causes one of the most common molecular abnormalities in acute myelogenous leukemia that involves induction of CD200 expression. An additional significant value (p = 0.014) was obtained by comparing the progenitor list with the signature of genes that are upregulated in double-positive lymphocytes, expressing both CD4 and CD8 coreceptors versus their double-negative counterparts (Dik et al., 2005).

In Vitro and In Vivo Studies of the Relationship between the MRUCD200/CD200R1 and MRUnot CD200/CD200R1 Populations

CD200 and CD200R1 are expressed in epithelial cells of the basal and luminal layers of the mammary gland (Figure S3). Here we focused on the relevance of their expression levels in MRUs to unveiling the top of the mammary cell hierarchy. Sorted cells of the two MRU subpopulations were challenged for in vitro mammosphere generation under non-adherent conditions that indicate stem cell activity. The MRUCD200/CD200R1 subpopulation generated a significantly higher number of mammospheres than its MRUnot CD200/CD200R1 counterpart (Figure 5A). The proportion of basal, CK14-stained colonies that developed after seeding diluted cells of these populations was highest for the MRUCD200/CD200R1 cells (Figure 5B). A significant decrease of 16% in basal colony proportion was noted for the MRUnot CD200/CD200R1 cells, and a further 11% decline was detected for the unsorted cell population. The combined number of CK14- and CK18-stained colonies was higher by 18%–25% for the MRUnot CD200/CD200R1 subpopulation compared with MRUCD200/CD200R1.

Figure 5

MRUCD200/CD200R1 Generates More Mammospheres and Higher Relative Number of Basal Colonies than MRUnot CD200/CD200R1 in Culture

However, both populations initiate outgrowths that differentiate into lobuloalveolar-like structures with β-casein-synthesizing capability during pregnancy.

(A) Analysis of non-adherent mammospheres generated by the two MRU subpopulations. Chi-square values were calculated and significant differences (p < 0.002), indicated by an asterisk, were determined using Pearson test. Inset: typical mammosphere. Bar, 50 μM. Two additional experiments showed 2.7- and 2.0-fold higher numbers of mammospheres, respectively, in cultures of MRUCD200/CD200R1 compared with MRUnot CD200/CD200R1.

(B) Relative number of basal and luminal colonies generated by the two MRU subpopulations. Number of colonies/well in these experiments ranged from 24 to 166. Seven wells/population were analyzed by t test of three experiments. Inset: the smallest colony that was included in the analysis contained three cells. Bar, 50 μM.

(C–E) H&E-stained mammary paraffin sections from 17-day intact pregnant gland (control) and from outgrowths developed from the MRU subpopulations. The latter were transplanted into the cleared mammary fat pad of a virgin female that was mated and became 17 days pregnant. Bar, 50 μM.

(F–H) Alpha smooth muscle actin (αSMA) immunostaining of the structures described in (C–E). Bar, 50 μM.

(I–K) β-casein (β-Cas). immunostaining of the structures described in (C–E). Bar, 50 μM.

To further elucidate the relationship between the MRU subpopulations, we showed that transplantation of MRUCD200/CD200R1 populations gives rise to outgrowths containing MRUnot CD200/CD200R1 cells (Table S5). We then asked whether outgrowths developed from the downstream MRUnot CD200/CD200R1 population encompass MRUCD200/CD200R1 cells. Table S5 shows high CD200/CD200R1 expressors in the MRU fraction of outgrowths developed from MRUnot CD200/CD200R1. Next we examined putative differences in differentiation potential between outgrowths developed from MRUCD200/CD200R1 and MRUnot CD200/CD200R1. Both populations developed lobuloalveolar-like structures in late pregnancy (Figures 5C–5E). They were composed of αSMS-stained basal cells and luminal cells with β-casein-synthesizing capability (Figures 5F–5K). A much narrower lumen characterized the lobuloalveolar structures in the outgrowths relative to the intact gland. This difference may indicate a stem cell memory of their hormone-derived (virgin) state (Asselin-Labat et al., 2010), or result from lack of draining capacity via the nipple.

Of note, at this late stage of pregnancy, there was substantial infiltration of immune cells in the developing outgrowths in 42% of the glands, independent of their origin, which was not observed in the normal gland (Figures 5C–5E). It has been shown that immune cells are attracted by the chemokine receptor CCR6 (Boyle et al., 2015) and promote epithelial growth throughout the gland. Estrogen-stimulated macrophages are more potent than non-stimulated ones in promoting human outgrowth development in the de-epithelized mouse mammary gland (Fleming et al., 2012).

Discussion

Characterizing the top of the epithelial cell hierarchy in the mammary gland is highly important from several perspectives: understanding developmental and regenerative processes of the gland, identifying potential tumor-initiating cells, and manipulating the capabilities of the gland toward higher milk production. In this study, we identified a heterogeneous CD200highCD200R1high epithelial cell population. Transplantation assays demonstrated that ∼50% of its constituents enriched the CD24medCD49fhigh (MRU) fraction for cells with relatively high mammary repopulating activity. The higher repopulating potential of MRUCD200/CD200R1 compared with MRUnot CD200/CD200R1, which encompasses the rest of the MRUs, was reflected by the larger area of the reconstituted mammary epithelium. However, it did not include changes in the number of reconstituted glands. This minor, but significant difference became highly pronounced when gene expression was profiled. RNA-seq analysis of gene expression in 1,000 pooled mammary epithelial cells clearly distinguished between the two MRU subsets, as well as additional CD200+CD200R1− and CD200−CD200R1+ populations. With this tool, two major aspects were addressed: (1) a possible hierarchical relationship between the MRU subpopulations versus the concept of two different types of stem cells; and (2) a putative linkage between immunorelated gene expression and stemness.

The much higher number of activated genes and metabolic pathways detected in the MRUnot CD200/CD200R1 indicated the coexistence of two metabolically divergent populations at the top of the mammary cell hierarchy. Among the numerous pathways supporting the higher metabolic activity of MRUnot CD200/CD200R1, the most highly expressed ones were involved in differentiation and function of progenitors, suggesting that the MRUnot CD200/CD200R1 are not stem cells (Figure S4). For example, eNOS and nitric oxide signaling have been associated with osteogenic differentiation through activation of the downstream canonical WNT/β-catenin signaling pathway (Bandara et al., 2016). Higher levels of reactive oxygen species have been recently identified as a marker of progenitors when compared with stem cells in the mammary gland (Diehn et al., 2009). The strongly activated Tec kinase signaling is also highly important for T cell differentiation (Andreotti et al., 2010).

Deeper insight into the core of the mediating proteins provided further support for pro-differentiating activity in the MRUnot CD200/CD200R1. For example, MAPK, which was highly activated in this population, has been associated with mammary branching morphogenesis and the expansion of cells positive for K6, a marker of hyperproliferative progenitor cells (Fata et al., 2007). Elevated MAPK levels also mark common progenitors that are downstream of the stem cell population in wild-type mice (Godde et al., 2014). PI3 kinase is also highly activated in the MRUnot CD200/CD200R1. Its expression in the mammary gland is regulated by the signal transducer and activator of transcription STAT5A (Schmidt et al., 2014), which marks mammary progenitor activity (Yamaji et al., 2009), and negatively regulates the expression of pTEN, which is more highly expressed in MRUCD200/CD200R1. PKC is also known to promote epithelial progenitor cell differentiation (Rieger et al., 2016) and RAS promotes osteoprogenitor cell proliferation and bone formation (Papaioannou et al., 2016). These data converge into the IPA-defined higher metabolic activity of MRUnot CD200/CD200R1, which is needed for cell growth, proliferation, development, and movement. Validating the concept of a single hierarchy in which MRUCD200/CD200R1 cells give rise to common MRUnot CD200/CD200R1 progenitors that differentiate into CD200+CD200R1− and CD200R1+CD200− cells yielded a list of 144 genes exhibiting decreasing expression patterns. The high compatibility with the previously characterized subset of stem cell markers in mouse and human glands (Lim et al., 2010), as well as with stem cell regulators (Chakrabarti et al., 2014), further supports this hypothesis. Importantly, most of the genes that were common to these two lists converged into specific collagen-associated and muscle-associated gene complexes that are probably involved in determining stem cell characteristics. Cdh3, a central linking gene for these complexes, was one of the most divergently expressed genes, at 1.8-fold higher levels in MRUCD200/CD200R1 compared with MRUnot CD200/CD200R1. Its potency as a stem cell marker is therefore worth analyzing. Taken together, the 39 common genes, together with 7 genes that serve as stem cell markers in other tissues, may provide a core list for mammary stem cell identification.

While the use of CD200 and its receptor appeared suitable for identifying stem cells and close progenitors in the mammary epithelial cells, the downstream hierarchy defined by their expression did not include all epithelial cell populations or even defined luminal or basal ones. Rather, it was restricted to an as-yet undefined immune-related cell population.

The major differences between MRUCD200/CD200R1 and MRUnot CD200/CD200R1 gene-expression profiles and resulting pathway activities led to an accumulated phenotypic outcome (Figure S4). Transplantation of MRUCD200/CD200R1 resulted in better filling of the cleared fat pad. Successful serial transplantation characterized only this population, indicating stem cell activity (Asselin-Labat et al., 2010). MRUCD200/CD200R1 also maintained a 3-fold higher potential to generate mammospheres under non-adherent conditions than its MRUCD200/CD200R1 counterpart. A more basal origin, depicted by a higher proportion of CK14-stained colonies that decreased toward the MRUnot CD200/CD200R1 and further on toward the unsorted cells, may be also associated with the MRUCD200/CD200R1. Nevertheless, de-differentiation capability of the MRUnot CD200/CD200R1 (Clevers and Watt, 2018) or a steady state between the two subpopulations (Lanner and Rossant, 2010) at an early developmental stage cannot be excluded, due to the ability of MRUnot CD200/CD200R1 to re-form outgrowths containing MRUCD200/CD200R1. Both MRU subpopulations formed outgrowths containing basal and functional luminal cells, depicting their close hierarchical state and multipotent capabilities. As such, the data presented here assemble the MRUCD200/CD200R1 and MRUnot CD200/CD200R1 as putative stem cells and common progenitors, respectively (Figure S4), that together contribute to mammary homeostasis (Rios et al., 2014, Sun et al., 2014).

The second aspect of this study involves the possible link between innate immune-related gene expression and stemness. The poor metabolic activity of the MRUCD200/CD200R1 population was accompanied by significantly higher activity of the complement system as compared with cells representing the rest of the MRUs. The complement system has long been perceived as an effector arm of innate immunity that mediates important immunoregulatory and inflammatory functions (Mastellos and Lambris, 2002). In this respect, higher activity of both CD200 and complement may be difficult to reconcile, because the former is considered to be an inhibitor of the complement system (Elward and Gasque, 2003) and the CD200–CD200R1 interaction on myeloid cells has been reported to attenuate immune activity (Minas and Liversidge, 2006). Other studies, however, suggest that the complement system may assume a non-inflammatory role in certain biological settings (Mastellos and Lambris, 2002). Indeed, complement has been implicated as a mediator of lens and limb regeneration in lower vertebrates (Kimura et al., 2003) and of mammalian liver regeneration via C3A and C5A activation (Strey et al., 2003). Others have shown a critical role for complement proteins (mainly C5a) in hematopoietic stem cell mobilization from their bone marrow niche to the blood, which was associated with the complement system's circadian activity (Ratajczak, 2015). Here we demonstrated that at least six complement proteins of the classical pathway that composes the membrane attack complex—C1 to C6—are highly expressed in the stem cell core of the mammary gland. Whether these proteins are involved in maintaining mammary epithelial stemness or have a role in the initial steps of stem cell differentiation remains to be elucidated.

Experimental Procedures

Mice

FVB/N or C57BL mice were housed under a 12-hr light/dark cycle and given food and water ad libitum. For all surgical procedures, mice were anesthetized with isoflurane (Abbott Laboratories, Maidenhead, England) mixed with O2 using a veterinary anesthesia machine. All animals used in this study were treated humanely. Study protocols were in compliance with the regulations of the Israeli Ministry of Health and local institution policies.

Dissociation of Mammary Tissue into Single-Cell Suspension and Flow Cytometry

Mouse mammary cells were dissociated from the fourth inguinal mammary glands of 6- to 8-week-old FVB/N virgin females. The dissociation procedure was as previously described (Rauner and Barash, 2012). Further details can be found in Supplemental Experimental Procedures.

Cell Transplantation

The endogenous mammary epithelium was surgically removed bilaterally from #4 mammary glands of 21-day-old female mice weighing 10–12 g (i.e., “clearing”). Sorted mammary epithelial cell populations were transplanted as detailed in Supplemental Experimental Procedures.

Serial Transplantation

A similar transplantation procedure was performed for serial transplantations of CD200high/CD200R1high MRUs. Host mice were C57BL and donor mice were C57/GFP (ubiquitin-EGFP mice [Reichenstein et al., 2016]). Further details can be found in Supplemental Experimental Procedures.

Outgrowth Analysis

For whole-mount examination, transplanted mammary fat pads were excised from sacrificed mice, fixed, and stained with Carmine Red as previously described (Rauner et al., 2013). Stained whole mounts were visualized and photographed using the binocular equipped with CellSens standard 1.4 software.

RNA Extraction and RNA-Seq Analysis

Dispersed mouse mammary cells were sorted into a 384-well plate as described in Supplemental Experimental Procedures. An aliquot of 1,000 cells in 1 μL from each population was collected into 9.5 μL RNasin lysis buffer (SMART-Seq V4 Ultra Low Input RNA Kit; Clontech, Mountain View, CA), incubated for 5 min at room temperature and kept at −80°C. For RNA-seq analysis, samples were thawed. Reverse transcription and cDNA amplification (12 cycles) of the full transcriptome were performed with the SMART-Seq V4 Ultra Low Input RNA Kit according to the manufacturer's protocol. Clean-up reactions were performed with Ampure XP beads (Beckman Coulter, Brea, CA). The amplified cDNA products were sheared by Covaris E220X (Woburn, MA) and 4.5–12.0 ng of sheared amplified cDNA from each sample was processed as previously described (Blecher-Gonen et al., 2013). An individual barcode was ligated to each sample to allow multiplexing of 20 libraries on two sequencing lanes. Between 10 and 12 million single-end 60-base pair reads were sequenced per sample on an Illumina HiSeq 2500 v4 instrument.

Histological Analysis and Immunostaining

Immunostaining of paraformaldehyde-fixed cells in culture and of paraffin-embedded or frozen tissue sections were performed using relevant antibodies (Table S6) as previously described (Rauner and Barash, 2012) with the modifications detailed in Supplemental Experimental Procedures.

Clonal and Mammosphere Assays

These assays were essentially performed as previously described (Rauner and Barash, 2012). Further details can be found in Supplemental Experimental Procedures.

Bioinformatics Analysis

Reads were trimmed using cutadapt (Martin, 2011) and mapped to the GRCm38 genome using STAR (Dobin et al., 2013) v2.4.2a (default parameters). Counting proceeded over genes annotated in Ensembl release 82 using htseq-count (Anders et al., 2015) (mode: intersection-strict). Differential expression analysis was performed using DESeq2 (Love et al., 2014) with the betaPrior. Cook's distance cutoff and independent filtering parameters were set to False.

Statistics

Unless otherwise indicated, Student's t test was performed for statistical analyses.

Author Contributions

G.R., T.K., and I.B. designed the concept of the experiments. G.R. and T.K. performed the flow cytometry. S.G. performed the RNA-seq. G.H. performed the core of the bioinformatics analysis. T.K. performed data analysis. I.B. performed additional bioinformatics analysis and wrote the manuscript.