Authors: Fanny A. Pelissier Vatter, Denis Schapiro, Hang Chang, Alexander D. Borowsky, Jonathan K. Lee, Bahram Parvin, Martha R. Stampfer, Mark A. LaBarge, Bernd Bodenmiller, James B. Lorens
Categories: Article, aging, breast cancer, heterogeneity, human mammary epithelia, mass cytometry, single-cell analysis
Source: Cell Reports
Aging is associated with tissue-level changes in cellular composition that are correlated with increased susceptibility to disease. Aging human mammary tissue shows skewed progenitor cell potency, resulting in diminished tumor-suppressive cell types and the accumulation of defective epithelial progenitors. Quantitative characterization of these age-emergent human cell subpopulations is lacking, impeding our understanding of the relationship between age and cancer susceptibility. We conducted single-cell resolution proteomic phenotyping of healthy breast epithelia from 57 women, aged 16–91 years, using mass cytometry. Remarkable heterogeneity was quantified within the two mammary epithelial lineages. Population partitioning identified a subset of aberrant basal-like luminal cells that accumulate with age and originate from age-altered progenitors. Quantification of age-emergent phenotypes enabled robust classification of breast tissues by age in healthy women. This high-resolution mapping highlighted specific epithelial subpopulations that change with age in a manner consistent with increased susceptibility to breast cancer.
Keywords: human mammary epithelia, aging, mass cytometry, single-cell analysis, heterogeneity, breast cancer
Adult tissue stem and progenitor epithelial cells generate differentiated daughter cells for tissue remodeling and homeostasis (Biteau et al., 2008, Mansilla et al., 2011). Evidence suggests skewed stem cell function contributes to diseases of aging (Sharpless and DePinho, 2007). Human breast epithelium, comprised of apical luminal epithelium (LEP) and basal myoepithelium (MEP) cell layers surrounded by a basement membrane, undergoes remarkable growth and remodeling between puberty and menopause and during lactation, supported by stem and progenitor cells. The greatest risk factor for breast cancer is age and exclusively genetic explanations are inadequate (LaBarge et al., 2016, Stephens et al., 2012). Differentiation-defective progenitor cells accumulate and tumor-suppressive MEPs decline with age, while older LEPs display basal properties, such as nuclear-localized YAP and MEP gene expression (Chen et al., 2014, Garbe et al., 2012, Pelissier et al., 2014, Skibinski et al., 2014). We hypothesize that these age-associated changes elevate cancer risk. Congruently, LEPs from women with high cancer risk (e.g., BRCA1/2 carriers) show basal characteristics, and luminal progenitors with a basal phenotype are suggested cells of origin for murine mammary adenocarcinoma (Lim et al., 2009, Molyneux et al., 2010, Proia et al., 2011). We therefore sought to gain insight into molecular changes in the mammary epithelium during aging and comprehensively catalog age-emergent phenotypic diversity using mass cytometry (Bandura et al., 2009) in samples from women aged 16–91 years old.
To measure age-emergent phenotypic diversity in the human breast, we used mass cytometry to obtain single-cell proteomic profiles of cryopreserved normal primary human mammary epithelial cell (HMEC) strains at passage four, from 44 women of ages 16 to 91 years old (Figure 1; Table S1). A 29-antibody panel recognizing human mammary epithelial lineage markers and intracellular signaling proteins was used to establish high-dimensional phenotypes of single HMECs (Figure S1; Table S2) (dos Santos et al., 2013, LaBarge et al., 2007, Lim et al., 2009, Regan et al., 2012, Taylor-Papadimitriou et al., 1989, Villadsen et al., 2007). Simultaneous analysis of ≥20,000 HMECs from each of the women from three age groups (<30 years, n = 16; >30 years < 50 years, n = 13; and >50 years, n = 15) measured 29 protein epitope dimensions (Ds) (Figure 1A). Non-linear dimensionality reduction, t-distributed stochastic neighbor embedding (tSNE) (Amir et al., 2013), created a 2D map of the entire dataset at single-cell resolution with similar phenotypes proximal to each other (Figure 1B). Distinct LEP (K19+/K7+/K8/18+/CD133+) and MEP (K14+/K5/6+) cell populations were distinguishable on the tSNE map (Figure 2A). Several signaling markers showed lineage dependence, especially in MEPs; CD44, YAP, phospho-epidermal growth factor receptor (pEGFR), pStat1, pS6, and p-Phospholipase C Gamma 2 (pPLCγ2), previously implicated in myoepithelial function and contractility (Pasic et al., 2011), were prevalent in the MEP cell population (Figures 2A, S1B, and S1C). The partial superposition of K14+ and K5/6+ cell population revealed the high degree of heterogeneity. A small subpopulation of cells with low marker expression also was noted, which may have corresponded to the K14−/K19− epithelia previously described (Villadsen et al., 2007), and it was not further examined.
Figure 1 Mass Cytometry Analysis of Human Mammary Epithelial Cells(A) Summary of experimental design. 57 samples of HMECs from women aged 16 to 91 years old (comprising 13 uncultured breast epithelia samples and 44 primary cultured HMEC strains at passage 4 [p4]) were barcoded and stained using a panel of 29 antibodies labeled with isotope tags and analyzed using mass cytometry.(B) Strategy to analyze high-dimensional single-cell data and identify lineage and age-related phenotypic divergence.See also Tables S1 and S2 and Figure S1.
Figure 2 Collective tSNE Analysis Distinguishes Major Luminal and Myoepithelial Lineages(A) The raw data have been transformed with arcsinh with the cofactor of 5. tSNE maps from HMECs at p4 from women <30 years old (merged and subsampled at 50,000 cells, n = 16), >30 < 50 years old (n = 13), and >50 years old (n = 15).(B) Log
2fold change of marker expression of LEP over MEP manually gated from tSNE projection map in HMECs from women <30 years old, >30 < 50 years old, and >50 years old. Data are log2of ratio of median ± SD.(C) K19 and K14 expression in LEP as a function of age. 250MK, 90P, 245AT, 173T, and an outlier 42P were excluded from the analysis.See also Figures S1–S3.
To address how aging affects the human mammary epithelium, we first demarcated LEP and MEP populations using a manual K19 gate (Figure 2A). K19, K7, K8/18, CD133, and cKit expression was high in LEPs, while K14 and K5/6 expression was higher in MEPs (Figure 2B). LEPs exhibited higher phospho-nuclear factor κB (pNF-κB), which is implicated in mouse mammary epithelial proliferation and branching (Brantley et al., 2001). MEPs had increased expression of basal markers (Axl, pS6, pPLGγ2, pEGFR, CD44, pGsk3, and pStat1) involved in myoepithelial homeostasis (Pasic et al., 2011). The phenotypic space projections on the tSNE maps were similar among the three age groups (Figure 2A), however, the expression levels of a number of markers changed significantly with age (Figures S2 and S3). The most prominent age-related difference was observed in the LEP population, where K14 and YAP expression increased and K19 and K7 decreased with age (Figures 2C and S3B). Overall these results revealed remarkable phenotypic heterogeneity within the mammary epithelia.
The tSNE map displayed regions of cell density (Figure 2A) and K14, K5/6, pRb, and CyclinB1 expression within the MEP lineage, indicative of distinct cellular subpopulations (Figures 2A and S1C). Intra-lineage subpopulations were identified as distinct clusters of cells with shared phenotypes using PhenoGraph (Levine et al., 2015) (Figures 3A and S4A). There were four LEP (LEP1–LEP4) and nine MEP (MEP1–MEP9) clusters (Figure 3B). One subpopulation, denoted double positive (DP, between 120 and 503 non-doublet cells), co-expressed K14 and K19 in a separate phenotypic space between the LEP and MEP populations. This DP population likely comprised epithelial progenitors (Villadsen et al., 2007). A small subpopulation of cells with low marker expression was not further examined (<0.52 ion counts per cell). The LEP3 cluster showed high levels of pRb and CyclinB1, indicative of higher proliferation compared to the other LEP subpopulations (reviewed in Giacinti and Giordano, 2006). Clusters MEP4 and MEP7 expressed higher levels of CyclinB1 that correlated with higher DNA content (iridium intercalator counts; Figure S4B) compared to the other MEP subpopulations.
Figure 3 Age-Related Phenotypic Divergence in the Landscape of HMECs(A) Heatmaps of marker expression in PhenoGraph clusters of HMECs from women <30 years old (Z score scale, merged, n = 16) (excluding 250MK, 90P and 245AT, 173T).(B) tSNE projection of the PhenoGraph clusters identified with PhenoGraph identified in (A), colored by cluster.(C and D) Heatmaps of marker expression in each PhenoGraph cluster in HMECs from (C) women >30 and <50 years old and (D) women >50 years old, normalized to values from <30-year-old women.(E) Plots of cell percentage in each PhenoGraph cluster (excluding 250MK, 90P and 245AT, 173T). Data are mean ± SEM.(F) Intra-sample heterogeneity for each woman is represented graphically by a horizontal bar in which segment lengths represent the proportion of the sample assigned to each cluster, colored accordingly (excluding 250MK).(G) The first two components of correspondence analysis (CA), accounting for 70% of the co-association structure between PhenoGraph subpopulations and different strains. Proximity among women and among clusters indicates similarity, however, only a small angle connecting a woman and a cluster to the origin indicates an association. The angle between women <50 years old and LEP was statistically smaller than the angle between women <30 years old and women >30 and <50 years old and LEP (t test, p < 0.001). PhenoGraph subsets are displayed as triangles and HMEC samples as circles.(H) Contributions of the PhenoGraph subpopulations to CA-1 and CA-2.See also Figure S4.
Age-related changes in marker expression were observed mainly within the LEP subpopulations. Heatmaps of marker expression in each PhenoGraph cluster, in HMECs from women >30 and <50 years old (Figure 3C) and women >50 years old (Figure 3D), were normalized to values from <30-year-old women to highlight age-related changes. Increased K14 and decreased K19 expression was observed with age in LEP2, LEP3, and LEP4 clusters from women >30 and <50 years old and in all LEP subpopulations from women >50 years old. In addition to phenotypic changes with age, the abundance of the LEP clusters significantly increased, whereas abundance of MEP2, MEP5, and MEP8 clusters significantly decreased with age (Figure 3E). This trend was observed at the individual level, with high inter-sample heterogeneity (Figure 3F). We previously reported age-related changes in LEP and MEP cells in vivo based on K14/K19 staining, and 4 lineage markers (Garbe et al., 2012) did not discern the degree of heterogeneity apparent in this new analysis. Prominent changes in marker expression and abundance occurred in three of four LEP types as early as middle age, and all four types change beyond 50 years. Indeed, the abundance of LEP1 increased more than 3-fold. Decreased abundance of MEP also was type specific.
Correspondence analysis (CA) provided a global understanding of the relationships between all PhenoGraph clusters and the age factor (Härdle and Simar, 2007). CA reduces high-dimensional observations to a smaller set of explanatory components, allowing visualization of data on each woman and PhenoGraph subsets in the same space (Figure 3G). Women >50 years old were associated with LEP1–4 subsets and women <30 years old were associated with MEP1–9 subsets, probably reflecting the relative abundance of those lineages with age. The DP subset, which represents progenitor cells, was associated mainly with older women. The first component, contributing 43.2% and comprising mainly LEP1, captured the tendency of older women to have more LEP (Figures 3G and 3H). The second component (27.5%) provided a different ordering. Altogether, there was a significant association between an age-dependent luminal subset and the chronological age of the primary epithelial cells.
Unsupervised agglomerative hierarchical clustering (Citrus) was used to examine age-dependent changes in an orthogonal manner. Multidimensional single-cell data were distilled to a hierarchy of marker expression-related clusters, and cluster-specific cell frequency changes were determined (Bruggner et al., 2014). Seven clusters were identified (Figures 4A–4C) that were significantly more abundant with age (prediction error of 26% as estimated by cross-validation and a p value < 0.05 using a Student’s t test) (Figure 4A; Figure S4D), all of which represented the LEP compartment. Figure S4C illustrates the agglomerative clustering. The LEP subpopulations that showed age-dependent changes had specific marker expression signatures consistent with acquired MEP/basal-like characteristics (Figures 4A and 4B; Figure S4G). The age-emergent LEP clusters were all higher in K14 compared with the <30-year LEP. Cluster A, residing at the apex of the hierarchy, was K19^low^ and K14^high^ (Figures S4E and S4G). Clusters B, C, and D showed higher YAP, HER2, cKit, Axl, pS6, pPLCγ2, pEGFR, CD44, pGSK3, pNF-κB, pAkt, pERK1/2, pMEK1/2, pStat1, pStat3, and pStat5 expression than <30-year LEP. Most of these markers are associated with proliferation and migration and are mainly expressed in young MEP. Each cluster had defining characteristics, e.g., cluster B had the highest pRb and CyclinB1 expression that correlated with higher DNA content (Figure 4B; Figure S4F). Only cluster H decreased in abundance with age (Figures 4A, 4C, and S4D), and it mapped to the MEP compartment of the tSNE landscape. Those cells expressed low levels of K14, pS6, CyclinB1, and pRb (Figure S4F), possibly indicative of a quiescent, terminally differentiated MEP.
Figure 4 Evidence of Age-Dependent Phenotypic Divergence in the Luminal PopulationThe Citrus algorithm was applied to identify cell populations by hierarchical clustering of phenotypically similar cells from an aggregate dataset from all samples (excluding 250MK, 90P and 245AT, 173T). A defining characteristic of each cluster is denoted as cluster A, K14^high^/K19^low^; cluster B, proliferative LEP; cluster C, basal LEP; cluster D, pS6 luminal; cluster E, Keratin^low^; cluster F, Keratin-low S; cluster G, Axl^high^; and cluster H, sub-MEP.(A) Boxplots of cell abundance in each age-related cluster and its representative tSNE phenotypic projection. Each data point on these graphs represents the proportion of the cluster cell number compared to the total cell number in a single sample. The log10 scale represents an abundance of cells from 0 to 1.(B) Heatmaps of marker expression of each cluster normalized to LEP from <30-year-old women for clusters A to G and MEP from <30-year-old women for cluster H.(C) Hierarchical tree of agglomerative clusters obtained with the Citrus analysis. Node sizes are scaled on the basis of frequency of cells in each cluster.See also Figure S4.
Collectively, these results indicated that a subset of LEP acquires a basal phenotype and accumulates while a subset of MEP decreases in abundance with age.
To confirm our findings, we conducted mass cytometry profiling using the 29-antibody panel on epithelial cells (≥10,000 cells per sample) derived from uncultured breast epithelia samples obtained from 13 women of different ages (<30 years, n = 7 and >50 years, n = 6; Table S1). As predicted by the HMEC analysis, the tSNE map revealed extensive phenotypic heterogeneity in the breast epithelia (Figure 5A; Figure S5A). Unsupervised clustering identified four distinct phenotypes of LEP (LEP1–4), seven types of MEP (MEP1–7), a DP subpopulation (between 13 and 719 non-doublet cells), and a low-expressing cell phenotype (<18.32 ion counts per cell) (Figure 5B). Protein expression patterns were consistent with these phenotypic designations and with the HMEC analysis (Figure 5C). K19, K7, K8/18, CD133, and cKit expression was high in LEP, while MEP showed higher expression of basal markers (K14, K5/6, Axl, pS6, CD44, pEGFR, and pStat1). The abundance of the LEP1 subpopulation significantly increased, whereas the abundance of MEP2 significantly decreased with age (Figure 5D), a trend also observed at the individual level (Figure 5E). Citrus analysis identified three clusters that were significantly more abundant with age (Figure 5F; Figure S5B), all residing within the LEP compartment of the tSNE phenotypic landscape. All three LEP clusters (A, B, and C) showed age-dependent changes in specific marker expression signatures, consistent with acquired MEP/basal-like characteristics (Figure 5G; Figure S5C). To quantify the extent of acquired MEP-like/basal phenotype, we calculated the geometric distance between the breast epithelia or HMEC LEP Citrus clusters and their respective MEP populations. This demonstrated that the LEP-MEP phenotypic distance was reduced by 26.7% and 32% in the breast epithelia and HMECs with age, respectively (Figure 5H). Collectively, this finding supports the notion that age-emergent epithelial cells derived from uncultured breast epithelia samples showed phenotypes that matched those identified in the primary cultured HMECs.
Figure 5 Age-Related Phenotypic Divergence in Uncultured Breast Epithelia(A) tSNE maps from dissociated uncultured breast epithelia from women <30 years old (merged and subsampled at 50,000 cells, n = 7) and >50 years old (merged and subsampled at 50,000 cells, n = 6). The pGsk3 channel was removed from the analysis due to a technical issue.(B) tSNE projection of the PhenoGraph clusters. The tSNE projection (right panel) of women <30 years old (blue) and women >50 years old (green) is shown.(C) Heatmaps of Z score of marker expression in PhenoGraph clusters of uncultured breast epithelia from women <30 years old (merged, n = 7).(D) Plots of cell percentage in each PhenoGraph cluster. Data are mean ± SEM.(E) Intra-sample heterogeneity for each woman is represented graphically by a horizontal bar in which segment lengths represent the proportion of the sample assigned to each cluster, colored accordingly.(F) Boxplots of cell abundance in each age-related Citrus cluster and its representative tSNE phenotypic projection.(G) Heatmaps of marker expression of each cluster normalized to LEP from <30-year-old women for clusters A to C and MEP from <30-year-old women for clusters D to G.(H) The geometric distance was calculated using the square root of the sum of the squared differences between the median of each marker for each subpopulation.(I) Representative human breast sections immunostained for K14 (red), K19 (green), and DAPI (blue) from a 17-year-old, 36-year-old, and 58-year-old woman (left to right, respectively). Scale bar represents 100 μm.(J) Plots show classification performance of 171 breast sections from 50 women (<30 years n = 52, >30 < 50 years n = 86, and >50 years n = 33), analyzed using morphometric context with increasing training set size.(K) Plot shows Citrus classification performance using a training set of 10 women (<30 years n = 5, >50 years n = 5). The black and white circles indicate whether an incorrectly assigned sample was from a peripheral non-tumor mastectomy (P), milk (MK), a tumor (T), or a tissue with no history, respectively. Data are mean ± SEM.See also Figures S5 and S6.
As the K14^high^K19^low^ clusters from both breast tissue and HMECs formed the apex of the age-dependent cluster hierarchy, we hypothesized that the expression pattern of these cytokeratins could be used to predict the approximate age of normal breast tissue. Human breast sections were stained with anti-K14 and anti-K19 (<30 years, n = 52 [10 women]; >30 < 50 years, n = 86 [25 women]; and >50 years, n = 33 [15 women]) (Figure 5I), and a classification model was built using morphometric context (Chang et al., 2013). At least 1,000 cells per section were analyzed. This computational approach relied on automated cell segmentation, with manual curation, to define different epithelial cells prior to quantification of single-cell K14 and K19 levels and morphometric features. The machine learning-based classification model correctly assigned more than 50% of the samples into their correct age group, as compared with a random guess of 33.3% (Figure 5J), based on the higher level of K14 and lower level of K19 in LEP with age, as observed on the tissue sections (Figure 5I). These data validated predictions from the mass cytometry data, and the in situ analysis demonstrated quantifiable age-related changes in LEP in breast tissue.
Next, we used age-dependent phenotypic divergence to build a second classification model to test the hypothesis that age-related changes in marker expression from our statistical analysis would generalize to an independent dataset. This second model was based on the totality of the mass cytometry data, and it was not restricted to K14 and K19. In general, classification models use cross-validation to avoid testing hypotheses suggested by the data (type III errors). Using a training set of 5 HMECs each from women in the <30-year and >50-year age groups, we successfully assigned 13/16 women <30 years old and 12/15 women >50 years old (Figure 5K). The classification performance was increased with the number of training samples (Figure S6A). Strikingly, the HMEC strains that were incorrectly classified as “old” were derived from tissue peripheral to a breast tumor, harbored a known high cancer risk mutation, or had LEP proportions above the mean, while a >50-year HMEC strain with decreased LEP proportion was incorrectly assigned as “young.” Hence, the classification model validated the hypothesis that subsets of luminal cells change with age and can be used as age predictors, and it suggested that this information could be a relevant indicator for cancer risk.
Next we investigated the functional consequences associated with the age-related basal phenotypic changes in the mammary luminal compartment. MEPs form cell-cell contacts with both LEPs and other MEPs and adhere to the basement membrane (Bergstraesser et al., 1995, Pitelka et al., 1973). These adhesive properties are characteristic of MEP, whereas LEP-extracellular matrix interactions are relatively minimal (Cerchiari et al., 2015). We hypothesized that the increased basalness of older LEPs would affect their cell adhesion and migratory capacity. Cell migration kinetics was measured by real-time impedance in LEP and MEP cells, sorted via fluorescence-activated cell sorting (FACS), from 6 different primary HMEC strains. LEPs isolated from women <30 years old migrated faster than isogenic MEPs (Figures S6B and S6C). In contrast, LEPs from women >50 years old migrated much slower and at a rate comparable to the cognate MEPs (Figures S6B and S6C), consistent with increased basal adhesion properties.
In addition, EGFR-mitogen-activated protein kinase (MAPK) pathway activation (pEGFR, pMEK, pErk, and pAkt) was higher in MEPs than in LEPs (Figure 2A; Figure S1), and it was higher in the LEP clusters that changed in abundance with age (Figure 4B). Therefore, we evaluated age-dependent differences in HMEC responses to EGF. HMECs from 3 women <30 years old and 3 women >50 years old were treated with EGF combined with vanadate, and signal transduction was measured by mass cytometry (Figures S6D–S6F). Both LEPs and MEPs exhibited EGFR pathway activation within the time course (Figure S1B); however, pStat, pEGFR, pErk, pMEK, and pPLCγ2 levels were increased in LEPs from women >50 years old compared to younger women. tSNE analysis revealed a subpopulation of HMECs with activated EGFR (Figure S6E). LEPs derived from older women were more prevalent in this population, consistent with increased EGF-signaling capacity (Figure S6F). Thus, older LEPs acquired myoepithelial-like adhesion and migration characteristics, as well as an increased EGF signaling.
We next asked whether the prominent age-related changes in LEPs were the result of changes in luminal progenitors (LPs). Luminal-biased progenitors expressing cKit (Lim et al., 2009) were FACS enriched from HMEC strains derived from 3 women <30 years old and 3 women >50 years old at passage 4, and they were analyzed by mass cytometry using 12 mammary epithelial progenitor and lineage markers (Figure 6A). Unsupervised agglomerative clustering identified two LP clusters that were more abundant (A and B) and one cluster that was less abundant (C) with age (Figure 6B; Figure S7A). These clusters localized within the luminal compartment of the tSNE phenotypic space, and they displayed specific marker signatures (Figure 6C; Figure S7D): cluster A was K7^high^ while cluster B was K19^high^. Both of these clusters displayed high expression of CD133, cKit, HER2, and YAP, establishing an age-emergent LP marker signature; thus, cKit+/CD133− and cKit+/CD133+ LP cells were FACS enriched from HMECs at passage 4. A higher proportion of cKit+/CD133+ LPs were present in older women (Figure 6D). The older cKit+/CD133+ LPs generated more acini in a Matrigel/collagen 3D embedding assay (Figure 6E), which incorporated more 5-ethynyl-2′-deoxyuridine (EdU), a proxy for cell proliferation (Figure 6F), compared to the corresponding younger LPs. In addition, cKit+CD133+ LPs were more luminally biased than cKit+CD133− LPs in younger women (Figure 6G; Figure S7B), as demonstrated by a higher proportion of K19+/K14− cells in the organoids. In contrast, older cKit+CD133+ LPs showed a higher proportion of basal K14+/K19− cells (Figure 6H; Figure S7C). Overall, the older cKit+/CD133+ LPs had higher clonogenic activity in vitro and gave rise to cells with more MEP/basal-like characteristics, which is consistent with the interpretation that altered LPs give rise to the LEPs that bear the phenotypic hallmarks of aging mammary epithelial cells.
Figure 6 Evidence of Age-Dependent Phenotypic Divergence in HMEC Progenitors(A) tSNE maps from FACS-enriched HMEC cKit+ progenitors from 3 women <30 years old and 3 women >50 years old (merged and subsampled at 50,000 cells). The lower right tSNE map shows the spatial projection of women <30 years old (blue) and women >50 years old (green).(B) Boxplots of cell abundance of age-dependent clusters identified with Citrus and their representative tSNE spatial projection.(C) Heatmaps of marker expression of each cluster compared to the background.(D) Proportions of cKit+CD133− and cKit+CD133+ as a function of age (n = 3; t test, young cKit+CD133− versus cKit+CD133+ p = 0.0046, and young cKit+CD133+ versus old cKit+CD133+ p = 0.037).(E) Acini formation potential of cKit+CD133− and cKit+CD133+ in Matrigel/collagen I gels as a function of age (n = 3; t test, p = 0.0123).(F) Proportions of acini that were incorporating EdU as a function of age. Data are means ± SEM. An acinus was quantified as EdU positive if at least one cell was incorporating EdU (n = 3; t test, p = 0.0452). Data are means ± SEM.(G and H) Histograms represent log
2-transformed ratios of K14 to K19 protein expression in single cells of acini (G) from a representative woman <30 years old (240L, 19 years) and (H) from a representative woman >50 years old (029, 68 years). Histograms are heat mapped to indicate the phenotypes of K14−/K19+ LEP (green), K14+/K19+ progenitors (yellow), and K14+/K19− MEP (red). Insets show representative HMEC organoids immunostained for K14 (red), K19 (green), and DAPI (blue). Scale bar represents 50 μm.See also Figure S7.
Normal epithelial cells must bypass tumor-suppressive barriers to give rise to malignancies. Pre-stasis HMECs can be efficiently immortalized in a two-step process that bypasses Rb function (by CCND1 expression or CDKN2A knockdown) and reactivates telomerase activity (indirectly by MYC expression), while incurring no gross genomic changes (Garbe et al., 2014). Mass cytometry analysis was conducted on 6 immortalized HMEC cell lines and visualized by tSNE. Immortalization via CCND1 overexpression to bypass stasis was associated with a luminal subtype, whereas knockdown of CDKN2A was associated with a basal subtype, and age >60 years, independent of CDKN2A, favored more luminal subtypes to emerge, consistent with our previous report (Lee et al., 2015) (Figures 7A–7C). Strikingly, five of the six immortal cell lines exhibited high expression of the basal markers YAP, Axl, pS6, pPLCγ2, pEGFR, CD44, and pGSK3 (Figures 7D and S7E), which is the same protein cluster observed in the subset of LEPs that accumulated with age (Figures 4B and 5G). Thus, the specific marker expression signatures found in age-emergent LEP subpopulations resemble the immortalized derivatives of older HMECs. This is consistent with the hypothesis that accumulation of altered LPs and LEPs with basal traits during mammary gland aging reflects a breast cancer susceptibility phenotype.
Figure 7 Aged Epithelia Resemble Immortalized Epithelial Cells(A) tSNE map of immortalized HMECs (left, merged, 6,000 cells per sample, n = 6). Right: each color represents a strain.(B) Five selected markers are shown (K19, K4, K7, Axl, and YAP), with knockdown of CDKN2A: p16sMY and overexpression of CCND1: D1MY.(C) Plots show percentage of LEP and MEP in each strain according to the gating strategy.(D) Heatmap of Z score of median of marker expression of each strain.See also Figure S7.
At single-cell resolution, these data show the dynamic phenotypic heterogeneity of human mammary epithelia spanning eight decades of life. Phenotypic diversity is present within all cell populations, and en bloc-averaged behavior may not represent that of individual cells (Altschuler and Wu, 2010). Mapping the normal diversity of cellular phenotypes within an adult tissue is key to understanding organ-level function and cell-level functionality. Using unbiased computational analyses of 29-parameter mass cytometry, we interrogated epithelial cell lineage diversity in HMEC and uncultured human breast epithelia samples from 57 individuals. We show that the mammary epithelium comprises a complex population of cells residing in phenotypically and functionally diverse states that change with age, and it is more dynamic and heterogeneous than previously perceived (Santagata and Ince, 2014, Taylor-Papadimitriou et al., 1989, Villadsen et al., 2007). This unique data resource provides a repository of single-cell proteomic data combined with cell functional and in situ tissue validation to better understand the aging process in mammary epithelia. An important outcome of the high-dimensional comparison, between breast epithelia and primary cultures of HMEC strains, was the excellent correspondence between lineage representation and phenotypes of aging. In general, a challenge of aging research in human tissues is that age-dependent changes can be cataloged (e.g., transcriptional, proteomic, and biochemical), but it is difficult to study the functional consequences of those changes. As shown here, age-related functional changes in migration, EGF signaling, and 3D colony formation, which were predicted by measurements of the 29 epitopes in breast epithelia, could be explored using primary cultures. Further development of these types of approaches, which combine careful primary HMEC culture with advanced analytical instrumentation and computational analysis, represents a means of accessing human tissue biology in a new and meaningful way.
Older women accumulate mammary multipotent progenitors, the putative cancer cell of origin, with altered differentiation potential, which we and others hypothesize is one of the mechanisms that underlies increased susceptibility to breast cancer with age (Choudhury et al., 2013, Garbe et al., 2012, LaBarge et al., 2016, Proia et al., 2011). Their LEP daughters are not fully lineage committed and acquire a basal phenotype. By using a test set of appropriate size and machine-learning tools, we identified a unique subset of LEPs that accumulate with age, which have skewed differentiation and a basal-like phenotype. These age-emergent LEPs are more adherent to the extracellular matrix (ECM), and they exhibit increased EGF signaling compared to younger LEPs. An unsupervised classification model based on these age-related markers validates the hypothesis that specific age-related changes in LEP subpopulations are good age group predictors. Concomitant age-related changes in LP cells support the concept that the altered stem/progenitor cell populations that accumulate during aging give rise to the age-dependent LEP populations.
The unique constellation of protein levels and modification states that enable classification of mammary epithelia according to age constitutes a signature of aging in the mammary gland. A primary component that distinguishes LEP clusters by age is increased expression of the basal cytokeratin K14 with decreased luminal K19 expression. In situ validation of the mass cytometry-derived classification model, by K14 and K19 quantification in tissue sections of normal breast, robustly assigned most breast biopsy samples into their correct age group, further supporting the observation that increased basalness of the LEP compartment is a hallmark characteristic of aging breast tissue. The reduction of LEP migration rate with age, approaching that of MEP, was consistent with increased ECM engagement. Change in the LEP-ECM binding energy is predicted to impair the ability to maintain normal epithelial bilayers (Cerchiari et al., 2015). Indeed, LEP-ECM interactions in older LEP could inhibit the Hippo pathway (Cordenonsi et al., 2011) and activate YAP/TAZ (Naylor et al., 2005), which may help explain our previous observation of increased nuclear YAP in post-menopausal LEP in vivo (Pelissier et al., 2014).
Accumulation of stem/progenitor cells with skewed differentiation and function is a hallmark of aging in a number of tissues (Encinas et al., 2011, Garbe et al., 2012, Lugert et al., 2010), and it may confer increased susceptibility to oncogenic events. BRCA1 basal-like breast cancers may originate from cKit-expressing progenitors (Lim et al., 2009, Molyneux et al., 2010). The key changes in the subpopulations of LPs that accumulated with age involved a specific marker signature where cKit, CD133, YAP, and HER2 expression was increased and CD44 expression was decreased. These findings are congruent with the observation that cKit overexpression prevents normal differentiation in murine mammary epithelial progenitors (Regan et al., 2012), and they may explain the accumulation of cKit+ LPs with age. CD133 (Prominin 1), another putative LP marker (Hilton et al., 2014, Raouf et al., 2008), changed the most with age (increased 18.25- and 6.36-fold and decreased 5-fold in protein abundance in each respective cluster compared to <30-year cKit+). In some breast cancers CD133 positivity correlates with a restricted subgroup of tumor stem cells in BRCA1-deficient mammary tumors (Wright et al., 2008). YAP expression in LPs might result in incompletely differentiated LEPs with basal traits (Pelissier et al., 2014). Amplification of HER2 has been shown to play an important role in the development and progression of certain aggressive types of breast cancer (Ménard et al., 2000). Hence, these age-related phenotypes correspond with the cancer-relevant hypothetical effects of cKit, CD133, HER2, and YAP marker expression signature. Older cKit+CD133+ LPs formed basal proliferative colonies in 3D Matrigel/collagen gels, resembling the reported activity of LPs from BRCA1 carriers (Lim et al., 2009). Moreover, hormone treatment in luminal breast cancer reduced estrogen receptor α levels and promoted a cancer stem cell phenotype (CD133^high^ and CD44^low^) (Sansone et al., 2016). This observation correlates well with the fact that augmented CD133^high^ and CD44^low^ cell proportion is associated with an increase in malignancy (Pece et al., 2010) and probably breast cancer risk (Garbe et al., 2012). The cKit+CD133+ LPs are considered potential cells of origin for breast cancer, and we propose that their accumulation with age represents an important facet of age-related susceptibility to breast cancer. Moreover, based on the involvement of these signature proteins (i.e., cKit, CD133, YAP, HER2, and CD44) in a number of breast cancer-relevant contexts, it is tempting to speculate that their dysregulation in older LPs also is related to increased susceptibility to breast cancer with age.
We identified shared expression signatures between age-emergent luminal cell subpopulations and immortalized derivatives of older HMECs, which have overcome at least two major barriers to tumor progression. Five of the six immortal cell lines exhibited high expression of the core markers (YAP, Axl, pS6, pPLCγ2, pEGFR, CD44, and pGSK3 [inactivated when phosphorylated]), which defines a putative signature of transformed HMECs that was also part of the core changes in normally aging LEPs. All of those core markers were reported to have a role in breast cancer progression (Gjerdrum et al., 2010, Kassis et al., 1999, Kim et al., 2015, Ko et al., 2016, Louderbough and Schroeder, 2011, Masuda et al., 2012, Yanai et al., 2015). Machine-learning algorithms that we used to validate the aging signature incorrectly classified nine young and middle-aged primary strains as old. Six of those strains were from women with mutations in BRCA1 or ATM, or they were derived from normal-appearing tissue that was peripheral to a tumor, or they were the normal cells cultured out directly from a tumor. It is tempting to speculate that epithelia in women who are innately at higher risk will exhibit an effective age that is older than the chronological age. A similar idea was proposed by epidemiologists that took into account breast tissue age, using only measures of hormones and childbirths (Pike et al., 1983). Here we provide evidence of a cell- and molecular-level manifestation of this concept. A biological explanation of advanced effective age in epithelia could be that molecular states that are associated with risk exert field effects (Deng et al., 1996), which alter the epithelia through a dynamic and reciprocal communication between stroma and epithelia. Ultimately, further investigation may lead to the development of novel approaches for prevention, patient stratification, and therapeutic interventions to combat age-associated breast cancers.
All cell culture was in M87A medium with 0.1 nM oxytocin (X) and cholera toxin (CT) at 0.5 ng/mL (Garbe et al., 2009). Primary HMEC strains were generated and maintained as described (Labarge et al., 2013). All tissues were obtained with proper oversight from the Lawrence Berkeley National Laboratory institutional review board. Breast tissue from reduction mammoplasty was manually dissected to enrich for gland-containing material. Stromal tissue was separated from epithelial fragments using a brief treatment with collagenase. The uncultured breast epithelia samples were dissociated as single cells with trypsin. All the pre-stasis HMEC strains were used at fourth passage (Table S1). Fibroblasts were removed and collected separately by differential trypsinization during the first passage. During the functional assay, HMECs were treated with EGF (Sigma E-9644, 0.1 μg/mL) and sodium orthovanadate (Sigma 13721-39-6, 12.5 mM) for 1 hr. Samples were harvested with TrypLE, fixed with 1.6% paraformaldehyde (PFA) for 10 min at room temperature (RT), and frozen as a pellet at −80°C for further analysis.
Antibodies were obtained in carrier protein-free PBS and then prepared using the MaxPAR antibody conjugation kit (Fluidigm), according to the manufacturer’s protocol. After determining the percentage yield by measurement of absorbance at 280 nm, the metal-labeled antibodies were diluted in Candor PBS Antibody Stabilization solution (Candor Bioscience) for long-term storage at 4°C (Table S2). Antibodies were titrated and validated beforehand using both positive and negative cell controls (Table S2; Figure S1). Extensive antibody validation has been performed and published previously (Chevrier et al., 2017, Giesen et al., 2014).
HMEC strains were incubated with cisplatin (WR International, Cat# 89150-634, 25 μM) for 1 min to assess cell viability (Fienberg et al., 2012), fixed in 1.6% PFA for 10 min at RT, and washed once with Cell Staining Media (CSM, PBS with 0.5% BSA and 0.02% NaN3 with 0.03% saponin). The cells were then resuspended in PBS, and DMSO stocks of the barcoding reagent were added as described (Bodenmiller et al., 2012, Zivanovic et al., 2014). The cells were incubated at RT for 30 min, washed three times with CSM, and then pooled into a single FACS tube for staining with metal-labeled antibodies for 1 hr at RT. A staining volume of 800 μL was used (∼30 × 10^6^ cells/mL). After antibody staining, the cells were washed twice with CSM and once with PBS, and then incubated for 20 min at RT or overnight at 4°C with an iridium-containing intercalator (DVS Sciences) in PBS with 1.6% PFA. The cells were then washed three times with CSM and once with PBS, diluted with water to ∼10^6^ cells/mL, and filtered through a 40-μm membrane just before analysis by mass cytometry.
The scale used before analysis is the arcsinh with the cofactor of 5 (x_transf = asinh(x/5)). After gating out viable and iridium-labeled events, the data were analyzed by applying tSNE. This non-linear dimensionally reduction technique is implemented via Barnes-Hut approximations in the MATLAB toolbox cyt (Amir et al., 2013). We used the default parameters (initial dimensions, 110; perplexity, 30; and theta, 0.5). Each sample contained 20,000 cells, when merged, 320,000 cells from HMEC <30 years, 260,000 cells from HMEC >30 < 50 years, and 300,000 cells from HMEC >50 years. In tSNE, each cell is represented as a point in high-dimensional space. Each dimension is one parameter (the expression level of each protein in our case).
The unsupervised PhenoGraph algorithm in cyt has been used to group cells that are phenotypically similar and cluster these subpopulations using modularity optimization (Levine et al., 2015). tSNE and PhenoGraph were performed only on surface markers. A number of neighbors of 800 was selected. This parameter was chosen based on prior knowledge of the underlying cell types. Lower values for nearest neighbors result in an overclustering and higher values an underclustering.
The Citrus toolbox in R was used to identify clusters that changed in abundance with age (Bruggner et al., 2014) in an unsupervised manner. Therefore, clusters were identified using a hierarchical clustering and linked to clinical data for characterization. The minimal selected cluster size was 0.1% of the total analyzed data. Stratifying clusters were learned by using regularized unsupervised learning methods. Heatmaps were obtained with MATLAB and Cytobank. The results were reproduced with strains obtained from reduction mammoplasties (RMs) only.
Citrus was implemented in the PAM package for R and used nearest shrunken centroids as a predictive model to identify properties that are predictive of sample class. The prediction model was based on the initial training data model. Therefore, new samples were mapped and later assigned to the initial clusters for prediction. Using a training set of 5 HMECs from women <30 years old and 5 women >50 years old (n = 10), Citrus efficiently assigned most of the test set to young or old. The training set was changed to n = 8 and to n = 12 with similar efficiency. After randomization of the training set, the classification failed. The R code is found in the Supplemental Experimental Procedures.
Each image was represented as its Cellular Morphometric Context (Chang et al., 2013), which was constructed as the histogram of cellular morphometric subtypes derived from the cellular morphometric features (K14/K19 signals) through K-Means (dictionary size = 1,024). Homogeneous kernel map (Vedaldi and Zisserman, 2012) was then applied on the Cellular Morphometric Context representation, so that linear support vector machine (SVM) (Yang et al., 2009) could be adopted for efficient and effective differentiation among age groups.
GraphPad Prism, R, and MATLAB were used for all statistical analyses. Standard linear regression and t tests were used. Grouped analyses were performed with Bonferroni-Holm correction for multiple comparisons. Significance was established when ^∗^p < 0.05, ^∗∗^p < 0.01, and ^∗∗∗^p < 0.001.
We thank Dr. Thibault Vatter, Dr. Robert Bruggner, and Dr. Adeeb Rahman for expert advice. The imaging and FACS were performed at the Molecular Imaging Center (MIC) at the University of Bergen. F.A.P.V. was supported by a University of Bergen pre-doctoral fellowship and travel grants. D.S. was supported by the Forschungskredit of the University of Zurich (grant FK-74419-01-01) and the BioEntrepreneur-Fellowship of the University of Zurich (BIOEF-17-001). J.B.L. is supported by grants from the Research Council of Norway (240130), Centres of Excellence funding scheme (223250), and Helse Vest Health Authority (911794). M.A.L. is supported by NIH R00AG033176 and R01AG040081 and a Congressionally Directed Medical Research Programs Breast Cancer Research Program Era of Hope Scholar Award BC141351. B.B. is supported by the Swiss National Science Foundation (SNSF) R’Equip grant 316030-139220, a SNSF Assistant Professorship grant PP00P3-144874, the PhosphonetPPM SystemsX grant (UC4 DK108132), and funding from the European Research Council (ERC) under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement 336921. D.S. was supported by a Forschungskredit Fellowship of the University of Zurich.
F.A.P.V., J.B.L., B.B., and M.A.L. designed the research. F.A.P.V. performed experiments. F.A.P.V., D.S., H.C., A.D.B., J.K.L., and B.P. analyzed the data and provided conceptual input. M.R.S., A.D.B., and M.A.L. provided cell strains, tissue sections, and other key reagents. B.B. provided the CyTOF platform. F.A.P.V., D.S., B.B., J.B.L., and M.A.L. wrote the manuscript.
The authors declare no competing interests.
Published: April 24, 2018
Mark A. LaBarge, Email: mlabarge@coh.org.
Bernd Bodenmiller, Email: bernd.bodenmiller@imls.uzh.ch.
James B. Lorens, Email: jim.lorens@uib.no.
The accession number for the CyTOF data reported in this paper is Mendeley Data: https://doi.org/10.17632/j7mrbgt3hh.1.