Intersection of FOXO- and RUNX1-mediated gene expression programs in single breast epithelial cells during morphogenesis and tumor progression

Edited by Yoshiaki Ito, Institute of Molecular and Cell Biology, Singapore, and accepted by the Editorial Board July 28, 2011 (received for review March 3, 2011)
August 22, 2011
108 (40) E803-E812

Abstract

Gene expression networks are complicated by the assortment of regulatory factors that bind DNA and modulate transcription combinatorially. Single-cell measurements can reveal biological mechanisms hidden by population averages, but their value has not been fully explored in the context of mRNA regulation. Here, we adapted a single-cell expression profiling technique to examine the gene expression program downstream of Forkhead box O (FOXO) transcription factors during 3D breast epithelial acinar morphogenesis. By analyzing patterns of mRNA fluctuations among individual matrix-attached epithelial cells, we found that a subset of FOXO target genes was jointly regulated by the transcription factor Runt-related transcription factor 1 (RUNX1). Knockdown of RUNX1 causes hyperproliferation and abnormal morphogenesis, both of which require normal FOXO function. Down-regulating RUNX1 and FOXOs simultaneously causes widespread oxidative stress, which arrests proliferation and restores normal acinar morphology. In hormone-negative breast cancers lacking human epidermal growth factor receptor 2 (HER2) amplification, we find that RUNX1 down-regulation is strongly associated with up-regulation of FOXO1, which may be required to support growth of RUNX1-negative tumors. The coordinate function of these two tumor suppressors may provide a failsafe mechanism that inhibits cancer progression.

Author Summary

Fig. P1.
Stochastic profiling of FOXO target genes reveals cross-talk with the transcription factor RUNX1. (A) Stochastic sampling strategy used to examine single-cell heterogeneities in FOXO target genes (red) for matrix-attached cells in 3D breast epithelial cultures. Housekeeping genes that exhibit more modest biological variations (blue) will not fluctuate as substantially between random 10-cell samplings. (B) Functional relationships between FOXO and RUNX1 expression in breast epithelia. RUNX1 loss or down-regulation causes a hyperproliferative phenotype in the presence of functional FOXO signaling, but it leads to oxidative stress-induced proliferative suppression when FOXO signaling is lost or inhibited.
Our work in breast epithelial cells reveals a functionally important cross-talk between two recognized tumor suppressors that have largely been studied in other tissue contexts. FOXO–RUNX1 cross-talk was uncovered by using stochastic profiling to observe correlated cell to cell fluctuations in mRNA expression. The study illustrates how in-depth, quantitative observations of this kind can help generate testable hypotheses. It remains unclear what upstream signals give rise to the observed heterogeneities in FOXO expression and RUNX1 phosphorylation in 3D culture. Our hope is that continued profiling experiments combined with mechanistic studies will enumerate all of the important transcriptional states that matrix-attached cells occupy during 3D morphogenesis. With this information, we can begin to examine how the different subpopulations relate to one another and explore whether similar subpopulations exist in breast tumors.
Finally, we asked whether the FOXO–RUNX1 dependencies seen in 3D culture were recapitulated in clinical breast cancer specimens. Both RUNX1 and FOXOs are recognized tumor suppressors, but our results predicted that reduction or loss of RUNX1 expression would require FOXO activity to enable cells to withstand the resulting oxidative stress. In a retrospective study of hormone-negative tumors that lacked amplification of the human epidermal growth factor receptor 2 oncogene, we found that tumors with the lowest RUNX1 expression had among the highest levels of FOXO1 expression. This negative correlation raises the possibility that the FOXO and RUNX1 transcriptional programs may intersect in vivo like we observed in 3D culture.
To test whether RUNX1- and FOXO-mediated gene expression were functionally coupled, we looked for phenotypes that emerged in 3D culture when RUNX1, FOXO, or both were perturbed. We found that RUNX1 knockdown caused cells to proliferate for longer in 3D culture and that disruption of FOXO function blocked the sustained proliferation. This genetic interaction was surprising, because FOXO target genes are often antiproliferative. To solve this puzzle, we examined the levels of reactive oxygen species during morphogenesis, because FOXOs have been implicated in protecting cells from oxidative stress (2). We found that FOXO disruption primed cells for substantial increases in reactive oxygen species when RUNX1 was knocked down. The resulting oxidative catastrophe caused a secondary proliferation arrest that restored normal 3D morphology, thereby explaining the observed genetic interaction between RUNX1 and FOXOs (Fig. P1B).
How could the regulation of the two FOXO groups differ? A simple explanation was that one of the groups was under control of a second class of transcription factors. Using bioinformatics, we identified a shared promoter sequence that distinguished one group from the other. This shared sequence contained a DNA binding site for the RUNX family of transcription factors. The breast epithelial cells used for 3D culture express RUNX1 and RUNX2, but RUNX1 mRNA was expressed at >15-fold higher levels than RUNX2, suggesting that RUNX1 was the dominant isoform. Furthermore, we found that RUNX1 promoter binding was strong in one FOXO group but not the other, which is in agreement with the bioinformatics-based prediction.
In our first application of stochastic profiling (1), we examined transcriptional heterogeneities among matrix-attached cells in a 3D culture model of breast epithelial cells. One candidate that we identified and provisionally validated was FOXO1. Because FOXO1 is itself a transcription factor that could lead to secondary heterogeneities in mRNA expression, we began our study here with a focused panel of 15 reported FOXO target genes that are expressed during 3D cell culture (Fig. P1A). After measuring the 10-cell sampling fluctuations of the panel, we discovered that most FOXO genes fell into one of two groups. Both groups were found to be heterogeneous, but their pattern of expression among single cells was markedly different.
To examine cell to cell differences in mRNA expression, we used a recently developed technique called stochastic profiling (1). This method examines the statistical fluctuations of gene expression measurements collected from random samplings of 10 cells (Fig. P1A). Among repeated 10-cell measurements, genes with levels that fluctuate more than expected are predicted to be heterogeneously expressed. After extracting candidate heterogeneities, we can cluster genes based on their pattern of sampling fluctuations to identify heterogeneous transcriptional programs.
Gene expression in cells is controlled by the recruitment of specific proteins to so-called promoter regions of the genome that coordinate the synthesis of mRNA from DNA, a process called transcription. Understanding transcriptional regulation is challenging, because protein occupancy on a gene's promoter is dynamic, complex, and variable between cells. Here, we exploited cell to cell heterogeneities in transcription to uncover a specific pattern of regulation between two transcription factor proteins, called Forkhead box O (FOXO) and Runt-related transcription factor 1 (RUNX1), which act as gene switches. In breast epithelial cells, FOXO and RUNX1 regulate transcription such that one or the other protein, but not both together, can act as a tumor suppressor. This study provides insight into the tumor-suppressive functions of proteins that modulate gene expression.
This article is a PNAS Direct Submission.
See full research article on page E803 of www.pnas.org.
Cite this Author Summary as: PNAS 10.1073/pnas.1103423108.

References

1
KA Janes, CC Wang, KJ Holmberg, K Cabral, JS Brugge, Identifying single-cell molecular programs by stochastic profiling. Nat Methods 7, 311–317 (2010).
2
GJ Kops, et al., Forkhead transcription factor FOXO3a protects quiescent cells from oxidative stress. Nature 419, 316–321 (2002).
Genetically identical mammalian cells often display patterns of gene protein expression that are profoundly different (13). Cell to cell heterogeneity has been recognized in cancer for nearly half a century (4, 5). However, only lately have heterogeneous cell populations been exploited as a means for uncovering new mechanisms of biological regulation (68). With recently developed techniques that can interrogate single-cell states, we are now poised to embrace heterogeneity rather than average it out (911).
Nonuniformities emerge at the earliest steps of gene expression (1214). Therefore, measurements of mRNA expression heterogeneity may help to uncover new biology if combined with systems approaches for analysis (11). This possibility offers a particularly rich opportunity for unraveling transcriptional networks, where complex combinations of factors work together to mediate programs of gene expression (15). Early pioneering studies initially focused on describing the fundamental kinetics and statistics of transcription in single cells (1, 1214, 1619). Now, an open question is how methods that provide single-cell information should be applied to help understand the coordination of transcriptional programs and their role in cell phenotype (20, 21).
For this purpose, we recently developed a technique called stochastic profiling (2). Stochastic profiling does not attempt to quantify mRNA levels in single cells but instead, gleans single-cell information by analyzing the statistical fluctuations in 15–20 repeated measurements of 10 cells. The increased amount of starting material allows high-quality, transcriptome-wide data to be obtained with conventional oligonucleotide microarrays after PCR-based amplification of cDNA from 10 microdissected cells (2, 22, 23). By surveying hundreds of cells overall, stochastic profiling captures the most reproducible cell to cell expression heterogeneities in a cell population. Individual transcripts that are nonuniformly expressed are then organized by the pattern of their expression fluctuations to reveal coordinated single-cell programs.
In this work, we report a detailed analysis centering around one class of transcription factors—the FOXOs—whose expression was found to be strongly nonuniform in our initial study (2). FOXO proteins are a subgroup of the Forkhead family of transcriptional regulators that play important roles in cell cycle arrest, stress responses, and cell death (24). All FOXO isoforms recognize a common (A/G)TAAA(T/C)A DNA consensus, which is frequently observed in the extended promoters of many genes (SI Appendix, Table S1) (25, 26). However, FOXOs are often not redundant, and isoform-specific functions have been widely documented (2729). FOXO transcriptional activity is regulated both positively and negatively by phosphorylation and ubiquitination (3034). Although FOXOs predominantly function as transcriptional activators, they may also act as repressors in certain contexts (35). Last, FOXOs can interact with at least a dozen other transcription factors, resulting in altered binding specificity and transcriptional activity (SI Appendix, Table S2) (36). Thus, in many ways, FOXO proteins are an archetype for the complexity that lies between signal transduction and gene expression.
Here, we asked whether a focused examination of endogenous FOXOs by stochastic profiling would yield insights into their function at the network level. Instead of using constitutively active FOXO alleles to homogenize signaling (SI Appendix, Table S3), we quantified fluctuations in FOXO expression and activity that occur naturally during 3D organotypic culture of breast epithelial cells (2, 37). By mapping these fluctuations onto a panel of FOXO target genes, we discovered that the measured FOXO expression signature divides into two groups with distinct single-cell patterns. One of these groups receives a key second input from another transcription factor, RUNX1, which is variably active. RUNX1 is required for the proper timing of proliferative suppression in 3D spheroids, and its stable knockdown results in abnormal hyperproliferative acinar structures. Interestingly, this phenotype requires native single-cell regulation of FOXO function to withstand the increased oxidative stress caused by RUNX1 knockdown. Combined inhibition of RUNX1 and FOXO leads to an acute state of oxidative stress that induces proliferation arrest and restores normal 3D morphology. In a large study of estrogen receptor (ER)-, progesterone receptor (PR)-, and HER2-negative breast cancers, we find that reduced RUNX1 expression correlates with FOXO1 up-regulation, which presumably enables tumor progression. Our results illustrate how careful single-cell analysis of gene expression can reveal functional interactions within transcriptional networks that are important for cancer-relevant cell phenotypes.

Results

Dissecting Single-Cell FOXO Coregulation by Stochastic Sampling.

We previously cataloged the transcriptional heterogeneities that emerge among individual matrix-attached breast epithelial cells during acinar morphogenesis in a 3D culture model (2, 37). One gene that was predicted with high confidence to be nonuniformly expressed was the transcription factor FOXO1. The pattern of single-cell FOXO1 expression correlated with many transcripts that had been previously linked to oxidative stress and proliferative suppression, suggesting a coordinated cellular response. Because FOXO family transcriptional activity is itself linked to oxidative stress and proliferation arrest (38, 39), we chose to investigate the interplay between the expression of FOXOs and FOXO-dependent target genes.
The MCF10A clone 5E (MCF10A-5E) used for 3D culture expresses not only FOXO1 but also FOXO3 (SI Appendix, Fig. S1), and we found that both FOXO proteins were heterogeneously expressed and localized in individual acini during morphogenesis (Fig. 1A). Together with the two FOXOs, we monitored time-dependent changes in the levels of 18 transcripts during morphogenesis: eight genes were validated or reported FOXO targets (BTG1, CAV1, CDKN1A, FBXO32, SEMA3C, SESN1, SOD2, and SOX4) (SI Appendix, Table S3), five genes were constitutively expressed (GAPDH, HINT1, PRDX6, S100A6, and UBC), two genes had nonuniform expression that correlated with FOXO1 at the single-cell level (BCL2L13 and CCNI) (2), and three genes were positively (CDKN1C and KRT10) or negatively (CCNB1) correlated with the time-dependent induction of many FOXO genes but have not been reported to be regulated by FOXOs (Fig. 1B). We found that most FOXO target genes were up-regulated together with FOXO1 and FOXO3, peaking at around day 10, whereas others showed delayed kinetics or were modestly down-regulated. This finding suggested that the gene panel could be induced by FOXOs and perhaps, also modulated by other transcriptional inputs during 3D morphogenesis.
Fig. 1.
Focused stochastic sampling of a FOXO expression dichotomy in 3D breast epithelial cultures. (A) Heterogeneous expression and localization of FOXO1 and FOXO3 proteins at day 10 of morphogenesis. Frozen sections were stained with the indicated antibodies and imaged by wide-field immunofluorescence. (Scale bar: 25 μm.) (B) Time-dependent induction of a FOXO gene panel during morphogenesis. MCF10A-5E cells were placed under morphogenesis conditions for the indicated times, and population-level mRNA measurements were performed by qPCR. Validated and reported FOXO target genes are shown in red and orange, respectively. Loading controls for qPCR normalization are shown in gray, and genes not previously implicated as FOXO targets are shown in black. Data are shown as the mean ± SEM of triplicate biological samples. (C) Stochastic sampling and qPCR profiling of the FOXO gene panel in matrix-attached cells at day 10 of morphogenesis. qPCR cycle thresholds for the indicated genes were mean-centered, scaled by their amplification efficiency, and clustered by Euclidean distance with average linkage. Samplings are numbered according to their organization in E. (D) Identification of candidate heterogeneities by stochastic sampling. Local gene neighborhoods were compared against a log-normal distribution with 30% coefficient of variation at a false discovery rate equal to 0.05 (yellow) as described (2). Note the clear distinction between genes predicted to be heterogeneously expressed (purple) and those genes whose fluctuations are consistent with background biological variation (gray). (E) Candidate coexpression groups organized by stochastic sampling. Sampling data from the candidate heterogeneities identified in D were scaled to unit log-normal variance and clustered by Euclidean distance with Ward's linkage.
We next sought to isolate the genes in the panel that showed strong expression heterogeneities in matrix-attached cells at day 10 near the peak of presumed FOXO activity (Fig. 1B). We adopted a stochastic profiling approach that involves quantitative gene expression measurements in random samplings of 10 microdissected cells (2). The global profiling method that originally identified FOXO1 is useful to broadly survey expression heterogeneities in cell populations (2). However, it is less powerful for a focused study, in which a modest collection of genes must be analyzed with high sensitivity and quantitative accuracy. Therefore, the stochastic-sampling principle was adapted to gene expression measurements obtained by real-time quantitative PCR (qPCR). The sensitivity of qPCR allowed us to omit the secondary reamplification step required for oligonucleotide microarrays, and the resulting measurements were quantitatively accurate across dozens of genes with differing abundances (SI Appendix, Fig. S2) (2). Combined with the superior dynamic range of qPCR (40), we expected that these modifications would reveal more subtle nonuniformities and coregulatory patterns that were overlooked by the global method (2).
We exponentially amplified cDNA from 18 groups of 10 matrix-attached cells and quantified expression fluctuations for the 20-gene panel (Fig. 1C). SEMA3C (a reported FOXO target gene) (35) and all of the constitutively expressed genes showed small sampling to sampling fluctuations that were not distinguishable from estimated biological variability (30% log-normal coefficient of variation) (41). However, the remaining 14 genes showed much larger fluctuations than could be explained by background variability, suggesting that these transcripts were heterogeneously expressed (Fig. 1D).
Next, we standardized the qPCR measurements of the candidate heterogeneities to have uniform variance and clustered these genes based on the pattern of their fluctuations (Fig. 1E). SOD2 (a validated FOXO gene) (42, 43) and two of the genes with no reported connection to FOXO fell outside the major cluster that contained both FOXO transcripts along with nine other genes. SOD2 expression is also controlled by NF-κB (44, 45), and our previous global analysis had shown SOD2 to cluster at the single-cell level with several other NF-κB inducers and target genes (2), suggesting an explanation for its distinct fluctuation pattern. Together, these findings reduced our 20-gene panel to 11 transcripts with strong single-cell heterogeneities that were associated with FOXO expression and activity.

FOXO Genes Show Two Distinct Patterns of Expression in Single Cells.

To validate the proposed organization of the gene panel (Fig. 1E), we developed a three-color RNA FISH procedure for quantifying the extent of coexpression between genes in single cells. We improved on two-color RNA FISH (2) by adding a third fluorescence reading from an equal mixture of probes against GAPDH, HINT1, and PRDX6 (Methods). These three genes were uniformly expressed (Fig. 1 B and C) and collectively provided an internal control, which could account for cell to cell differences in probe penetration and total RNA levels. Normalizing the two-color RNA FISH data to the loading control allowed us to distinguish pairs of genes with strong, weak, or no coexpression among matrix-attached cells (Fig. 2 AC).
Fig. 2.
Three-color RNA FISH reveals two groups of FOXO genes with decoupled single-cell expression patterns. (A–C) Multicolor RNA FISH provides a semiquantitative means for characterizing (A) strong, (B) weak, and (C) no coregulation of gene pairs within individual cells. Haptenylated riboprobes complementary to the indicated genes were hybridized to day 10 frozen sections of MCF10-5E acini and imaged by wide-field immunofluorescence. Fluorescence intensities for the individual genes were pseudocolored (Upper), and a mixture of three housekeeping genes was used as a loading control (Lower Left). (Scale bar: 25 μm.) The loading-normalized fluorescence intensities of each gene pair were correlated for single cells of three independent acini to quantify the extent of coregulation (Rsingle cell; Lower Right). (D) FOXO genes are coregulated at the single-cell level within groups 1 and 2 but not between groups. Pairwise single-cell correlations (Rsc) measured by RNA FISH were analyzed as in A–C and compared with the dendrogram derived from stochastic sampling (Fig. 1E). The complete RNA FISH dataset is shown in SI Appendix, Fig. S3. (E and F) Single-cell distributions of RNA FISH intensities from FOXO groups 1 and 2 are heavy tailed, indicating two expression dichotomies. The data from SI Appendix, Fig. S3 are shown as a quantile–quantile plot displaying the observed RNA FISH quantiles (after log transformation) vs. the quantiles of a standard normal distribution. If the data are consistent with a log-normal distribution, the quantiles will plot as a straight line, which was observed for the RNA FISH intensities for the homogeneously expressed gene UBC (SI Appendix, Fig. S3T). Both FOXO groups deviated significantly from a log-normal distribution because of more frequent than expected observations of very high and very low RNA FISH fluorescence (P < 0.001 by Jarque–Berra test).
Using RNA FISH, we systematically evaluated the extent of coexpression between gene pairs based on the organization proposed by stochastic sampling (SI Appendix, Fig. S3 A–S). As expected, many pairs within the large FOXO cluster were tightly correlated, and their coexpression with genes outside the cluster was weak or nonexistent (Fig. 2 A and D). More surprising was the clear partition of the FOXO cluster into two groups (Fig. 2D). Genes within each FOXO group were coexpressed, but between group correlations were measurably weaker or undetectable (Fig. 2 B–D). Although the split between the two FOXO groups did not exactly coincide with a branch point on the clustering dendrogram, the boundary was consistent with the overall dendrogram organization (46). We also retroactively inspected the original dataset and identified multiple 10-cell samplings, which indicated that the two FOXO groups were not expressed concordantly (Fig. 1E, samplings 1, 9, and 13). By compiling the measured RNA FISH intensities from thousands of single cells, we found that neither group conformed to a log-normal distribution expected for background biological variation (Fig. 2 E and F and SI Appendix, Fig. S3T) (47, 48). Genes from the two FOXO groups instead showed increased frequencies of high and low expression, which were presumably the cell states detected by stochastic sampling. The results from stochastic sampling and RNA FISH together suggested that the FOXO genes in the panel were subject to one of two modes of acute single-cell regulation.

Control of FOXO Group 1 Genes by the Transcription Factor RUNX1.

The simplest explanation for the distinction between FOXO groups 1 and 2 was that one group was receiving an inducible input from a transcription factor other than FOXO. The N termini of FOXOs can bind various other transcription factors (24), raising the possibility that a specific FOXO heterodimer might be important to one FOXO group. However, for no reported FOXO dimerization partner were we able to find a pattern of consensus sites within 5 kb of the transcription start that could discriminate FOXO groups 1 and 2 (SI Appendix, Tables S1 and S2) (49). This finding suggested that the second input might be entirely independent of FOXO localization and activity, thereby explaining the uncoupled regulation of the two FOXO groups.
To identify candidate inputs, we used the bioinformatics algorithm multiple expectation maximization for motif elicitation (MEME) to search for recurring sequence motifs in upstream promoter regions that discriminated the two FOXO groups (50). A guanine-thymine-rich motif was uncovered specifically among group 1 genes, which aligned with the TG(T/C)GGT consensus binding site for the RUNX family of transcription factors (Fig. 3A). Other RUNX family members have been implicated in breast cancer as oncogenes (RUNX2) or tumor suppressors (RUNX3) (5154). By contrast, RUNX1 has been mostly investigated in hematopoietic cells (55), where it also acts as a tumor suppressor whose function is commonly disrupted in leukemias (56, 57). Very recently, however, the RUNX1 genomic locus was shown to be lost in one malignant variant of the parental MCF10A line used for 3D culture in our study (58). This report also found that high-grade breast tumors have decreased RUNX1 mRNA levels compared with lower-grade tumors, although the functional consequences of RUNX1 down-regulation were not examined. Finally, multiple RUNX isoforms have been shown to interact with a constitutively active FOXO3 when overexpressed (59, 60). Using qPCR with a shared genomic DNA standard, we found that RUNX1 mRNA was expressed at >15-fold higher levels compared with RUNX2, and RUNX3 was absent (Fig. 3B). Taken together with our bioinformatics analysis, these studies suggested that RUNX1-mediated gene expression might have a role in acinar morphogenesis and could possibly intersect with FOXOs.
Fig. 3.
Heterogeneous RUNX1 signaling decouples FOXO groups 1 and 2. (A) MEME analysis identifies a discriminatory motif for group 1 genes that contains a consensus binding site for RUNX family transcription factors; 5 kb upstream and 1 kb downstream of the transcription start site for the genes in the two FOXO groups were analyzed for discriminative motif discovery by MEME (50). The TRANSFAC database was then searched for aligning consensus sites by TOMTOM, with the motif P value shown for RUNX (89). A similarly informative motif was not identified for FOXO group 2. (B) Relative mRNA copy numbers of RUNX family members in MCF10A-5E cells. Levels of endogenous RUNX1, RUNX2, and RUNX3 were compared with a universal genomic DNA standard and normalized to measured RUNX2 expression. HeLa cells, which express all RUNX family members, were used as a positive control. Data are shown as the mean ± SEM of triplicate biological samples. n.d., not detected. (C and D) RUNX1 protein is homogeneously expressed but heterogeneously phosphorylated at day 10 of morphogenesis. Frozen sections were stained with the indicated antibodies and imaged by wide-field immunofluorescence. (Scale bar: 25 μm.) (E) ChIP-based survey of FOXO and RUNX1 promoter occupancy for the genes in the two FOXO groups. FOXO and RUNX1 consensus sites within 5 kb on either side of the transcription start site are shown in green and yellow, respectively; 150- to 200-bp qPCR amplicons used in the study are shown flanked by 800-bp whiskers indicating the extent of chromatin shearing and thus, the region conceivably detected by each amplicon. (F) Summary of FOXO and RUNX binding sites validated by ChIP. Binding sites were considered validated if their median enrichment over a control IgG ChIP was twofold or more across four independent ChIPs. The number of binding sites between groups was compared by Wilcoxon rank sum test. Similar results were obtained when a more redundant RUNX consensus site, YGYGGTY (91), was used. (G) Knockdown of endogenous RUNX1 without compensatory changes in RUNX2. MCF10A-5E cells were infected with the indicated viral hairpins, and RUNX levels were determined by immunoblotting. β-tubulin was used as a loading control. (H and I) RUNX1 knockdown causes the two FOXO groups to become coregulated. Representative genes from FOXO groups 1 (BCL2L13) and 2 (CDKN1C) were analyzed by RNA FISH as described in Fig. 2 A–C. Single-cell correlations (Rsingle cell) across three independent acini are shown with 90% confidence intervals in parenthesis. (Scale bar: 20 μm.)
In our earlier global survey of expression heterogeneities during 3D morphogenesis, we detected RUNX1 mRNA but did not predict that its expression would be nonuniform among matrix-attached cells (2). At day 10 when stochastic sampling was performed, we observed homogeneous expression of RUNX1 protein by immunofluorescence (Fig. 3C). This finding confirmed the prediction of our earlier survey and suggested that group 1 nonuniformities could not be caused by heterogeneities in RUNX1 expression. RUNX1 activity is also modulated posttranslationally by proline-directed kinases such as ERK1/2 and cyclin-dependent kinases (CDKs) (6163). We consequently stained for RUNX1 phosphorylated at S276 (p-RUNX1) and found highly variable immunoreactivity among matrix-attached cells (Fig. 3D). p-RUNX1 levels were strongly reduced when acini were treated with the CDK inhibitor roscovitine, but the staining pattern was unaffected by treatment with U0126 to inhibit ERK1/2 activation (SI Appendix, Fig. S5). Thus, cell to cell differences in p-RUNX1 downstream of CDKs (63) could conceivably provide the basis for a second heterogeneity that could decouple the group 1 FOXO genes in the panel.
Although a RUNX1-containing motif discriminated FOXO group 1 from group 2 (Fig. 3A), RUNX1 consensus sites were located throughout the promoters of all genes in the panel (Fig. 3E). To distinguish actual RUNX1 and FOXO binding events, we used ChIP with antibodies specific for endogenous RUNX1, FOXO1, or FOXO3 (Methods). An important technical hurdle was to reconcile the typical sample requirements for ChIP (∼107 cells) with the typical 3D morphogenesis sample at day 10 (∼250,000 cells). As a substitute for 3D culture, we seeded MCF10A-5E cells on tissue culture plastic overlaid with a diluted concentration of the same reconstituted basement membrane protein gel that is used in 3D morphogenesis assays (37). Remarkably, we found that long-term 2D culture under these conditions mimicked the observed 3D gene expression kinetics for the FOXO gene panel (SI Appendix, Fig. S4). The overall FOXO expression program, therefore, did not require standard 3D culture explicitly and could be examined at the ChIP scale by 2D culture under these conditions.
For ChIP, we designed qPCR primers scanning nearly the entire 5-kb promoter region of each gene after considering the size of the input chromatin fragments (∼800 bp) (Fig. 3E). Binding sites were scored as ChIP-validated if their median enrichment from four independent ChIPs was at least twofold above a paired IgG control. The summary of the resulting 1,128 qPCR measurements indicated a significant enrichment for RUNX1 binding events in the promoters of group 1 genes compared with group 2 genes (P < 0.01 by Wilcoxon rank sum test) (Fig. 3F). Conversely, FOXO binding was observed for all genes in the panel, and there was no discrimination between groups based on the total number of binding sites. The heterogeneous activation of RUNX1 (Fig. 3D) together with its selective promoter occupancy on group 1 genes (Fig. 3F) strongly supported RUNX1 as the second input to the FOXO gene panel.
If RUNX1 activity is indeed critical for decoupling FOXO groups 1 and 2, then removing RUNX1 should cause the two groups to become coexpressed in single cells. We tested this hypothesis by stably knocking down RUNX1 and quantifying whether coexpression between the two FOXO groups was affected. We validated a lentiviral shRNA that uniformly reduced RUNX1 protein without compensatory changes in RUNX2 levels (Fig. 3G and SI Appendix, Fig. S6) and generated 3D acinar structures with cells stably expressing shRUNX1. For two intergroup gene pairings that showed weak to no correlation in control acini (BCL2L13CDKN1C and CAV1BTG1), we found that single-cell coregulation increased significantly with RUNX1 knockdown (P < 0.05 by Fisher z-transformed confidence intervals) (Fig. 3 H and I and SI Appendix, Fig. S7). These data strongly support the conclusion that RUNX1 is the dominant second input that uncouples the single-cell expression of FOXO genes as first detected by stochastic sampling.

FOXO1 and RUNX1 Jointly Coordinate Proliferative Arrest In Breast Epithelial Acini.

Aside from the perturbation of gene expression in FOXO groups 1 and 2, we did not observe any changes in proliferation, apoptosis, or polarity of MCF10A-5E acini with RUNX1 knockdown at day 10 of morphogenesis. By day 14, however, there was a clear distinction between shRUNX1 cells and controls. Whereas most control cultures had reduced proliferation by day 14 based on phospho-Rb (pRb) staining, acini derived from shRUNX1 cells were still actively proliferating (Fig. 4 A and B). Proliferative suppression was restored with roscovitine, implying that CDK activity was sustained in shRUNX1 acini (Fig. 4C). This phenotype was specific to RUNX1, because add back of RNAi-resistant murine Runx1 restored proliferative suppression (Fig. 4D). Besides sustained proliferation, acini derived from shRUNX1 cells were also often larger and abnormally shaped by day 14 (Fig. 4B). This finding suggested an inappropriate coordination of proliferation and arrest in developing acini at the single-cell level.
Fig. 4.
RUNX1 knockdown causes proliferative and morphogenetic defects that are blocked by homogenization of FOXO activity. (A) Add back of RUNX1 levels in shRUNX1 cells with murine RNAi-resistant Runx1. MCF10A-5E cells stably expressing murine Runx1 (lane 3) or vector control (lanes 1 and 2) were infected with shGFP (lane 1) or shRUNX1 (lanes 2 and 3), and RUNX1/Runx1 levels were determined by immunoblotting. β-tubulin was used as a loading control. (B) RUNX1 knockdown delays growth arrest and causes abnormally shaped acinar structures. MCF10A-5E acini were fixed at day 14, stained for pRb and HA-tagged murine Runx1 (Runx1), and analyzed by confocal immunofluorescence. (C) Delayed growth arrest caused by RUNX1 knockdown is blocked by the CDK inhibitor roscovitine. MCF10A-5E acini were treated with 5 μM roscovitine on day 10, fixed at day 14, and stained for pRb. The percentage of acini with five or greater pRb-positive cells was quantified in >200 acini per replicate. (D) pRb quantification of the RUNX1 knockdown and add-back experiment described in B. (E) Overexpression of a DN-FOXO1 in the context of RUNX1 knockdown. MCF10A-5E cells stably expressing an HA-tagged DN-FOXO1 (lanes 3 and 4) or vector control (lanes 1 and 2) were infected with shRUNX1 (lanes 2 and 4) or shGFP control (lanes 1 and 3), and RUNX1 levels were determined by immunoblotting. Expression of the DN-FOXO1 was verified by HA expression, and β-tubulin was used as a loading control. (F and G) DN-FOXO1 restores normal acinar shape and growth arrest in acini derived from shRUNX1 cells. Acini from the MCF10A-5E lines described in E were fixed at day 14, stained for pRb and HA-tagged DN-FOXO1, and analyzed by confocal immunofluorescence as described in B and C. Note that acini lacking DN-FOXO1 (orange) remain hyperproliferative and are misshapen. For C, D, and G, data are shown as the mean ± SEM of quadruplicate 3D cultures at day 14. Significance of the interaction between shRUNX1 and DN-FOXO1 in F was determined by two-way ANOVA. (Scale bar: B and F, 25 μm.)
Because RUNX1 down-regulation induced coordinate expression of the two FOXO groups (Fig. 3 H and I) as well as escape from growth arrest (Fig. 4 B and D), we asked whether FOXO was necessary for the shRUNX1-induced hyperproliferative phenotype. To address this question, we coexpressed shRUNX1 and a truncated FOXO1 that acts as a dominant negative (Fig. 4E and SI Appendix, Fig. S8) (64). In the context of RUNX1 knockdown, dominant-negative FOXO1 (DN-FOXO1) restored the normal proliferative suppression and shape of maturing acini (Fig. 4 F and G). Surprisingly, overexpression of DN-FOXO1 alone had no effect on growth arrest or acinar morphology, suggesting that the genes induced by heterogeneous FOXO signaling become critical only when RUNX1 function is impaired.

FOXO1 Protects Cells from Oxidative Stress Caused by RUNX1 Loss and Permits Disruption of Epithelial Function.

The observed FOXO requirement for shRUNX1-mediated escape from proliferative suppression seems paradoxical, because FOXO itself can induce many growth arrest genes (38). However, other FOXO genes may be more important for proliferation in the context of RUNX1 knockdown. For example, FOXO also regulates a program of genes that protect against oxidative stress (39), and if RUNX1 knockdown were to cause oxidative stress, then proliferation of knockdown cells could be contingent on proper FOXO function.
To examine the redox state of MCF10A-5E cells during morphogenesis, we stained live 3D cultures with the oxidation-sensitive dye, chloromethyl-2′,7′-dichlorofluorescein diacetate (DCFDA). In early control cultures (day 6), we consistently observed a fraction of acini whose matrix-attached cells stained heterogeneously for DCFDA, indicating nonuniform levels of reactive oxygen species (ROS) (Fig. 5 A and B). This observation differs from an earlier report showing no DCFDA staining in matrix-attached cells at this point during MCF10A morphogenesis (65). We attribute the discrepancy to our use of the 5E clone (2), which behaves somewhat differently in this context compared with the parental MCF10A line (SI Appendix, Fig. S9). Notably, however, both the 5E clone and parental line showed heterogeneous DCFDA staining in matrix-attached cells at day 10 when they were profiled for heterogeneities (SI Appendix, Fig. S9). This finding suggests that the timing of nonuniform ROS accumulation may be accelerated in the 5E clone, but the heterogeneity itself is a property of 3D MCF10A culture.
Fig. 5.
Dual inhibition of RUNX1–FOXO function causes oxidative stress that leads to proliferative suppression. (A and B) RUNX1 knockdown plus DN-FOXO1 synergistically increases global ROS levels in MCF10A-5E acini at day 6. Cells expressing shRUNX1 or an RFP-tagged DN-FOXO1 (Inset, red) were labeled with DCFDA (green) as described in Methods and counterstained with DRAQ-5 (blue) to label nuclei. The percentage of acini staining heterogeneously or homogeneously for DCFDA on a cell by cell basis is shown in B. The remaining acini were negative for DCFDA staining. (C–E) Treatment with the antioxidant Trolox restores hyperproliferation in shRUNX1 + DN-FOXO1 cells. Acini from the MCF10A-5E lines described in Fig. 4E were fixed at day 14, stained for pRb and HA-tagged DN-FOXO1, and analyzed by confocal immunofluorescence as described in Fig. 4 B and C. For B and C, data are shown as the mean ± SEM of quadruplicate 3D cultures at (B) day 6 or (C) day 14. Significance of the interaction between shRUNX1 and DN-FOXO1 was determined by two-way ANOVA. (Scale bar: A, 20 μm; D and E, 25 μm.)
In MCF10A-5E acini at day 6, we found that RUNX1 knockdown or DN-FOXO1 individually did not affect the frequency of heterogeneous DCFDA-positive acini. However, shRUNX1 cells showed a noticeable increase in baseline ROS levels, and the extent of cell to cell heterogeneity DN-FOXO1 cells was clearly exaggerated. In striking contrast, combined expression of shRUNX1 and DN-FOXO1 significantly increased the number of acini that stained strongly and homogeneously for DCFDA (Fig. 5 A and B). This finding indicates that dual inhibition of RUNX1 and FOXO1 synergistically causes widespread oxidative stress during 3D morphogenesis in vitro.
Sublethal oxidative stress causes proliferative suppression in many cell types, suggesting a mechanism for the observed 3D phenotype (Fig. 4 F and G) (6668). To determine whether ROS levels were important for growth inhibition in shRUNX1+DN-FOXO1 cells, we used the synthetic antioxidant Trolox. Prolonged 3D culture in Trolox had no effect on proliferation in control cultures at day 14 but eliminated suppression observed in shRUNX1+DN-FOXO1 acini (Fig. 5 C–E). Thus, normal FOXO function is required to dampen elevated levels of ROS caused by RUNX1 knockdown.

Loss of RUNX1 Expression Correlates with FOXO1 Up-Regulation in Human Breast Cancers.

Both RUNX1 and FOXOs are known to act as tumor suppressors in different contexts (56, 69). However, our in vitro results suggested that loss of RUNX1 or FOXO1 in breast cancer could be mutually exclusive, in that FOXO activity would be required to stabilize transformed cells with low RUNX1 expression. To test this prediction, we retrospectively analyzed a high-quality microarray study that included ER-, PR-, and HER2-negative breast cancers (Gene Expression Omnibus accession no. GSE6861). This triple-negative subset enriches for basal-like cancers, which have a molecular expression profile that is similar to the MCF10A lineage (70, 71). We used FOXO1 mRNA to read out both FOXO1 levels and overall FOXO activity as a target gene (72, 73) and examined the relative levels of expression averaged across two distinct probe sets. The chosen study was uniquely reliable in its estimates of RUNX1 and FOXO1 (R > 0.4 between probes for the same gene), possibly because the microarrays used were explicitly optimized for clinical specimens.
Across 89 triple-negative breast cancers, we found a significant anticorrelation between FOXO1 expression and RUNX1 expression (Fig. 6A) (P < 0.0005, Student t test). The relationship between FOXO1 and RUNX1 was notable, because negative correlations are, in general, less likely to be spurious and can often suggest functional connections (74). Among the 10% of tumors with the highest and lowest expression of RUNX1, we observed significant differences in FOXO1 levels: low RUNX1 expression was associated with high FOXO1 expression and vice versa (Fig. 6B) (P < 0.005, Wilcoxon rank sum test). The strong anticorrelation and reciprocal differences between RUNX1 and FOXO1 were not observed in the subset of 71 breast cancers that were not triple negative (SI Appendix, Fig. S10). Taking these clinical data together with our in vitro studies, we predict that FOXO1 expression and activity is critical for triple-negative breast cancers with low RUNX1 expression.
Fig. 6.
Low RUNX1 expression is associated with high FOXO1 expression in triple-negative breast cancer specimens. (A) RUNX1 and FOXO1 mRNA expression is anticorrelated in triple-negative breast cancers. Correlated probe sets from the Gene Expression Omnibus accession no. GSE6861 dataset were standardized as z scores and averaged to give the estimated relative levels of RUNX1 and FOXO1. The 10th and 90th percentiles of RUNX1 expression (dashed boxes) were further analyzed in B. (B) Tumors with the lowest RUNX1 mRNA expression have the highest FOXO1 mRNA expression. Box and whisker plots from the samples in the dashed boxes in A are shown with notches to indicate 90% nonparametric confidence intervals. Significance for the increase in FOXO1 expression was determined by the Wilcoxon rank sum test.

Discussion

The results of this study point to an unanticipated coupling between gene expression programs mediated by FOXOs and RUNX1. These transcription factors act as standalone tumor suppressors in endothelial and myeloid cells (56, 69) but converge during morphogenesis and tumor progression in breast epithelia. Our results suggest an isoform-specific role for RUNX1 in tumor suppression, because others have shown that RUNX2 can collaborate with the MYC oncogene to promote tumorigenesis (75). The connection discovered between FOXOs and RUNX1 was made possible by an in depth fluctuation analysis of single-cell gene expression in a tissue-like context. This work illustrates how biological mechanisms can be revealed simply by observing complex molecular patterns in situations where spurious correlations are unlikely (76).

High-Sensitivity Analysis of Single-Cell Expression Heterogeneity by Focused Stochastic Sampling.

Combining the principle of stochastic sampling with qPCR rather than microarrays was advantageous for the FOXO gene panel, because it allowed higher-sensitivity detection of transcriptional heterogeneity. Although indispensable for global analyses, microarrays tend to dampen expression differences compared with when the same samples are analyzed by qPCR (40). Indeed, one-half of the genes scored as heterogeneous in this study were missed by our earlier microarray-based survey (2). All of these genes had been clearly detected on the microarray platform and showed observable 10-cell fluctuations. However, the magnitude of sampling to sampling variation was deemed insignificant when tested amid the thousands of other genes on the array. The improved sensitivity and reduced false discovery of qPCR-based sampling make it an attractive alternative to microarrays for a focused study.

Molecular and Functional Interactions Between FOXOs and RUNX1.

During segmentation in Drosophila, the RUNX ortholog Runt interacts genetically with various transcriptional modulators to activate or repress expression of pair-rule genes (77). The genetic interactions that we observed between FOXOs and RUNX1 raise the question of whether these transcription factors might interact at the molecular level. FOXO3 has been reported to bind to RUNX1 and a related family member at the BCL2L11 promoter (59, 60). However, these earlier studies overexpressed RUNX proteins together with a constitutively active FOXO3 mutant. Despite repeated attempts, we were unable to coimmunoprecipitate endogenous RUNX1 with endogenous FOXOs, even when using modified lysis conditions with chemical cross-linking and gentle detergents. Furthermore, among the qPCR amplicons with measurable enrichment in our ChIP experiments, we found that only 33% were simultaneously enriched in anti-FOXO ChIPs and anti-RUNX1 ChIPs. Of these amplicons, 75% had DNA binding sites for both FOXO and RUNX1, consistent with two separate binding events. Together, these observations suggest that a direct FOXO–RUNX1 interaction might occur only when RUNX1 expression or FOXO activity is very high.
Instead, we favor a mechanism in which FOXOs and RUNX1 converge on a common set of genes to promote transcription independently. Thousands of genes have tandem FOXO–RUNX1 binding sites that are conserved across mammals (78, 79). Those sites with the closest spacing (<35 bp) are enriched for various gene ontologies associated with development or signaling (SI Appendix, Table S4). Secondary functional interactions are also likely to exist. For example, CDK-mediated phosphorylation of RUNX1 (SI Appendix, Fig. S5) would be subject to indirect regulation through CDK inhibitors, such as CDKN1A and CDKN1C. Transcription of these inhibitors may, in turn, depend on FOXOs and RUNX1 (Fig. 3F). Dual FOXO–RUNX1 coordination may be important in developmental contexts where rapid but precise expansion of cells is required, such as during 3D epithelial acinar morphogenesis.

Implications for Breast Cancer Progression and Therapy.

The inverse relationship that we observed between FOXO1 and RUNX1 levels in triple-negative breast tumors raises general questions about oncogenesis and cancer therapy. First, why is loss of RUNX1 function so commonly associated with acute myeloid leukemia (AML) (80)? This finding would seem to contradict the codependency between RUNX1 and FOXOs proposed here, because FOXOs are absent or excluded from the nucleus in nearly all AML patients (81). Because ROS levels are high in AML cells (82), one could speculate based on our studies that AML cells have adopted other mechanisms for tolerating oxidative stress. In support of this possibility, AML cells show substantial up-regulation of superoxide dismutases (SODs), which neutralize ROS (83). Interestingly, our earlier data in 3D cultures suggest that the manganese SOD gene (SOD2) is predominantly controlled by NF-κB rather than FOXO (2), and AML cells often show constitutive NF-κB activation (84). The tissue specificity of RUNX1 as a tumor suppressor may derive from the other aberrations in AML cells that allow them to withstand oxidative stress in the absence of FOXO activity.
How then might this influence therapy for breast cancers in which RUNX1 loss depends on FOXO activity? We predict that supplementary antioxidants would be particularly detrimental, because they would relieve the FOXO requirement and facilitate tumor progression that would not otherwise occur (65). One unconventional strategy might be to try evoking an oxidative catastrophe by acutely disrupting FOXOs and possibly, NF-κB. Blocking endogenous antioxidant responses may be most effective at forestalling aggressive tumors that are proliferating rapidly, which is the case in breast cancers with triple-negative status (85). Encouragingly, targeted methods for increasing ROS levels in solid tumors are already under development (86).

Methods

Plasmids.

pBabe DN-FOXO1-HA neo, pBabe RFP-DN-FOXO1-HA neo, and pBabe HA-Runx1-neo were constructed by PCR or subcloning with standard approaches. Cloning details are available in SI Appendix, SI Methods.

Cell Lines.

The MCF10A-5E clone was previously reported (2) and was cultured in 3D as described for MCF10A cells (37).

Lentiviral RNAi.

pLKO.1 shRUNX1 puro was obtained through the RNAi Consortium and contained the hairpin CCGGCCTCGAAGACATCGGCAGAAACTCGAGTTTCTGCCGATGTCTTCGAGGTTTTT, where the 21-mer sequence targeting human RUNX1 is underlined. pLKO.1 puro and pLKO.1 shGFP puro (Addgene) have been described previously (87). Lentiviruses were prepared and used to infect MCF10A-5E cells by the standard approaches described in SI Appendix, SI Methods.

Retroviral Overexpression.

Retroviruses were prepared in 293T cells (ATCC) by double transfection of the pBabe construct together with pCL ampho (Addgene) according to standard procedures. Viral particles were collected at 48 h and passed through a 0.45-μm filter before use. MCF10A-5E cells were infected as described in SI Appendix, SI Methods.

Immunoblotting.

MCF10A-5E cells expressing the indicated constructs were lysed in radioimmunoprecipitation assay (RIPA) buffer (50 mM Tris, pH 8.0, 150 mM NaCl, 5 mM EDTA, 1% Nonidet P-40, 0.1% SDS, 0.5% sodium deoxycholate); 20–30 μg clarified extract were separated on an 8% or 10% SDS/PAGE gel and transferred to PVDF (Millipore). Membranes were blocked with 5% nonfat skim milk in Tris-buffered saline-Tween (TBS-T) (20 mM Tris, pH 7.6, 150 mM NaCl, 0.1% Tween-20), washed one time in TBS-T, and incubated overnight at 4 °C in 5% BSA (Sigma) containing one of the following primary antibodies: anti-AML1/RUNX1 (1:1,000; Cell Signaling), anti-RUNX1 (1:1,000; Abcam), anti-RUNX2 (1:1,000; MBL), or anti-HA (clone 3F10, 1:500; Roche). Membranes were washed three times for 5 min each in TBS-T and incubated for 1 h at room temperature in 5% nonfat skim milk in TBS-T containing HRP-conjugated goat anti-rabbit or goat anti-rat secondary antibody (1:5,000, Santa Cruz). Membranes were washed three times for 5 min each in TBS-T and imaged by chemiluminescence on an LAS-3000 detection system (FujiFilm).

Frozen Sectioning, Laser Capture Microdissection, and Small-Sample mRNA Quantification.

Embedding, sectioning, microdissection, and small-sample mRNA amplification were performed at day 10 of 3D morphogenesis exactly as described elsewhere (2).

qPCR.

qPCR was performed on amplified small-sample material or with purified RNA as described elsewhere (2, 88). Primer sequences are available on request.

Immunofluorescence.

Immunofluorescence in frozen sections was performed as described elsewhere (2) using the following primary antibodies: FOXO1 (1:100 dilution; Cell Signaling), FOXO3 (1:200 dilution; Upstate), E-cadherin (1:500 dilution; BD Biosciences), AML1/RUNX1 (1:100 dilution; Cell Signaling), phospho-AML1/RUNX1 [S249 (S276 in the longer splice variant), 1:500 dilution; Cell Signaling], phospho-Rb (S807/811, 1:200 dilution; Cell Signaling and T821, 1:1,000 dilution; Epitomics), anti-HA (clone 3F10, 1:500 dilution; Roche), and phospho-ERK1/2 (T202/Y204, 1:200 dilution; Cell Signaling). Immunofluorescence on coverslips was performed as described in SI Appendix, SI Methods.
Whole-mount immunofluorescence of day 14 acini was performed essentially as described previously (37), except that 1× Western Blocking Reagent (Roche) was substituted for 10% goat serum in the primary block solution and the secondary block was omitted.

Multicolor RNA FISH.

Multicolor RNA FISH was performed as described elsewhere (2) with the addition of a mixture of Alexa 647-labeled probes (GAPDH, HINT1, and PRDX6 used in combination as a loading control) supplemented to the hybridization mixture at a concentration of 200 ng/mL. To prepare this probe mixture, aminoallyl-labeled riboprobes were first synthesized with 80% aminoallyl-UTP (Ambion) and 20% unlabeled UTP as described previously (2). Fluorescent labeling of aminoallyl-labeled riboprobes was performed with Alexa 647 reactive dye (Invitrogen) as recommended by the manufacturer. After purification on a Purelink PCR column (Invitrogen), labeled riboprobes were ethanol-precipitated, resuspended in RNase-free water to 0.2 μg/mL, and stored at −80 °C. By spectrophotometry, riboprobes were determined to contain 2–3 Alexa 647 dye molecules per probe.

ChIP.

MCF10A-5E cells were plated at 12,500 cells/cm2 in 15-cm dishes and cultured in morphogenesis medium (assay medium + 2% matrigel + 5 ng/mL EGF) for 10 d according to 3D culture conditions (37). Cells were fixed for 5 min at room temperature by adding formaldehyde to the culture medium to a final concentration of 1% (wt/vol), which was later quenched for 5 min at room temperature with 1/20 volume of 2.5 M glycine. Cells were washed two times in cold PBS, lysed in 300 μL lysis buffer [1% (wt/vol) SDS, 5 mM EDTA, 50 mM Tris⋅Cl, pH 8.0, 1× protease inhibitors], and transferred to a microcentrifuge tube. Lysates were incubated on ice for 10 min and then sonicated at 4 °C on a Bioruptor (Diagenode) for six 8-min rounds of 25-s pulses followed by 35-s intervals. After centrifugation at 12,000 × g for 20 min, the supernatant was collected, and 20 μL soluble chromatin were kept as the input fraction. Soluble chromatin was diluted 10-fold in dilution buffer (1% Triton X-100, 2 mM EDTA, 150 mM NaCl, 20 mM Tris⋅Cl, pH 8.0, 1× protease inhibitor) and then precleared by incubating 1 mL diluted chromatin with protein A-Sepharose (50 μL 50% slurry in 10 mM Tris⋅Cl, pH 8.1, 1 mM EDTA) for 2 h at 4 °C. ChIP-grade anti-FOXO1 (Abcam), anti-FOXO3 (Abcam), anti-AML1/RUNX1 (Abcam), or an equivalent amount of normal rabbit IgG was added to the diluted chromatin and incubated overnight at 4 °C with agitation; 50 μL protein A-Sepharose were added to the immune complexes and incubated for 2–4 h at 4 °C. Sepharose beads were collected and washed sequentially for 10 min each in RIPA (1×), RIPA plus 500 mM NaCl (3×), LiCl buffer (0.25 M LiCl, 1% Nonidet P-40, 1% deoxycholate, 10 mM Tris⋅Cl, 1 mM EDTA; 2×), and Tris-EDTA buffer (2×). DNA from the beads and the input fraction was eluted by reversing methylene cross-links with 500 μL elution buffer [0.5% (wt/vol) SDS, 200 mM NaCl, 10 mM Tris⋅Cl, 1 mM EDTA] at 65 °C overnight. Samples were treated with 100 μg/mL RNase for 30 min at 37 °C and 200 μg/mL proteinase K for 90 min at 50 °C, followed by extraction with phenol-chloroform. The aqueous fraction was ethanol-precipitated, washed one time in 70% ethanol, air dried, and dissolved in nuclease-free water. qPCR with genome-specific primers was used to amplify the ChIP-enriched DNA. Primer sequences are available on request.

ROS Imaging.

At day 6 of 3D morphogenesis, acini were washed one time for 5 min in PBS and incubated with 5 μM DCFDA (freshly prepared as a 1,000× stock in DMSO, Invitrogen) and 2.5 μM DRAQ5 (Cell Signaling) for 1 h at 37 °C. The acini were then washed one time for 5 min in PBS and fixed with 2% PFA for 20 min at room temperature. After three 5-min washes in PBS, slides were mounted with Prolong antifade medium (Invitrogen) and sealed with nail polish after curing.

Image Segmentation and Quantification.

Single cells from RNA FISH images were segmented by hand based on wheat-germ agglutinin staining (DAPI channel). Traced image segments were then applied to the digoxigenin- and dinitrophenyl-labeled riboprobe stainings (FITC and Cy3 channels) as well as the loading control staining (Cy5 channel). Mean fluorescence intensities per cell for each riboprobe were normalized to the mean fluorescence intensity of the loading control for that cell. This normalized fluorescence intensity was used as a relative measure of transcript abundance for quantifying single-cell coregulatory patterns of gene expression.
Acini in whole-mount specimens were scored positive for pRb staining if the structure contained at least five pRb-positive cells. For ROS imaging, acini were scored based on whether they stained negative, heterogeneously positive, or homogeneously positive for DCFDA fluorescence; 100–200 acini were scored per replicate for four independent morphogenesis cultures.

Bioinformatics Analysis.

To identify discriminatory motifs between the two FOXO groups, 5 kb upstream and 1 kb downstream of the transcription start site of each gene were analyzed by using MEME (http://meme.nbcr.net/) (50) in discriminative mode. Enriched DNA motifs were then fed into TOMTOM (89) to search the TRANSFAC database for similar transcription factor binding motifs.
For clinical bioinformatics, probe sets for FOXO1 (g9257221_3p_a_at and Hs.170133.0.A1_3p_s_at) and RUNX1 (g1932819_3p_a_at and g1932819_3p_x_at) were extracted from the triple-negative breast cancer samples in Gene Expression Omnibus accession no. GSE6861. Each probe set was standardized as a z score, and the mean z score across two probe sets was used to estimate the relative expression level per sample.

Statistical Analyses.

Hypothesis testing for single-cell heterogeneities was performed by the χ2 goodness of fit test on the binned sampling data at a false discovery rate of 5% as previously described (2). The RNA FISH data were tested for log normality by log-transforming the normalized, median-scaled measurements and using the Jarque–Berra test, which is a χ2 test based on a weighted sum of statistical moments estimated from the distribution (90). Pairwise testing of ChIP binding sites and differences in FOXO1 expression in triple-negative breast cancers were performed using the Wilcoxon rank sum test. Confidence intervals for Rsingle cell were calculated by the Fisher z transformation. Phenotypic interactions for the shRUNX1 + DN-FOXO1 experiments were examined by two-way ANOVA. The significance of the anticorrelation between RUNX1 and FOXO1 expression in triple-negative breast cancers was evaluated by the t test after the following transformation:
where ρ is the initial Pearson correlation of n samples and ρt is the transformed value. All tests were performed at α = 0.05.

Notes

The authors declare no conflict of interest.
See Author Summary on page 16499.

Acknowledgments

We thank Deirdre O'Toole for help with the cloning of the RNA FISH riboprobes, Silvia LaRue for cryosectioning, Laura Selfors for assistance with the tumor sample bioinformatics, and Dr. Richard Iggo for comments on the microarray data from the European Organisation for Research and Treatment of Cancer 10994 clinical trial (GSE6861). This work was supported by National Institutes of Health Grant 5-R01-CA105134 (to J.S.B.), National Institutes of Health Director's New Innovator Award Program 1-DP2-OD006464 (to K.A.J.), the Mary Kay Ash Charitable Foundation (K.A.J.), the Pew Scholars Program in the Biomedical Sciences (K.A.J.), and the David and Lucile Packard Foundation (K.A.J.).

Supporting Information

Appendix (PDF)
Supporting Information

References

1
JM Levsky, SM Shenoy, RC Pezo, RH Singer, Single-cell gene expression profiling. Science 297, 836–840 (2002).
2
KA Janes, CC Wang, KJ Holmberg, K Cabral, JS Brugge, Identifying single-cell molecular programs by stochastic profiling. Nat Methods 7, 311–317 (2010).
3
DK Singh, et al., Patterns of basal signaling heterogeneity can distinguish cellular populations with different drug sensitivities. Mol Syst Biol 6, 369 (2010).
4
IJ Fidler, ML Kripke, Metastasis results from preexisting variant cells within a malignant tumor. Science 197, 893–895 (1977).
5
CM Southam, A Brunschwig, Quantitative studies of autotransplantation of human cancer. Cancer 14, 971–978 (1961).
6
LH Loo, et al., Heterogeneity in the physiological states and pharmacological responses of differentiating 3T3-L1 preadipocytes. J Cell Biol 187, 375–384 (2009).
7
SV Sharma, et al., A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations. Cell 141, 69–80 (2010).
8
SL Spencer, S Gaudet, JG Albeck, JM Burke, PK Sorger, Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis. Nature 459, 428–432 (2009).
9
A Brock, H Chang, S Huang, Non-genetic heterogeneity—a mutation-independent driving force for the somatic evolution of tumours. Nat Rev Genet 10, 336–342 (2009).
10
M Niepel, SL Spencer, PK Sorger, Non-genetic cell-to-cell variability and the consequences for pharmacology. Curr Opin Chem Biol 13, 556–561 (2009).
11
SJ Altschuler, LF Wu, Cellular heterogeneity: Do differences make a difference? Cell 141, 559–563 (2010).
12
I Golding, J Paulsson, SM Zawilski, EC Cox, Real-time kinetics of gene activity in individual bacteria. Cell 123, 1025–1036 (2005).
13
A Raj, CS Peskin, D Tranchina, DY Vargas, S Tyagi, Stochastic mRNA synthesis in mammalian cells. PLoS Biol 4, e309 (2006).
14
J Yu, J Xiao, X Ren, K Lao, XS Xie, Probing gene expression in live cells, one protein molecule at a time. Science 311, 1600–1603 (2006).
15
S Komili, PA Silver, Coupling and coordination in gene expression processes: A systems biology view. Nat Rev Genet 9, 38–48 (2008).
16
D Zenklusen, DR Larson, RH Singer, Single-RNA counting reveals alternative modes of gene expression in yeast. Nat Struct Mol Biol 15, 1263–1271 (2008).
17
JR Chubb, T Trcek, SM Shenoy, RH Singer, Transcriptional pulsing of a developmental gene. Curr Biol 16, 1018–1025 (2006).
18
J Elf, GW Li, XS Xie, Probing transcription factor dynamics at the single-molecule level in a living cell. Science 316, 1191–1194 (2007).
19
SM Janicki, et al., From silencing to gene expression: Real-time analysis in single cells. Cell 116, 683–698 (2004).
20
SJ Gandhi, D Zenklusen, T Lionnet, RH Singer, Transcription of functionally related constitutive genes is not coordinated. Nat Struct Mol Biol 18, 27–34 (2011).
21
A Raj, SA Rifkin, E Andersen, A van Oudenaarden, Variability in gene expression underlies incomplete penetrance. Nature 463, 913–918 (2010).
22
NN Iscove, et al., Representation is faithfully preserved in global cDNA amplified exponentially from sub-picogram quantities of mRNA. Nat Biotechnol 20, 940–943 (2002).
23
MR Emmert-Buck, et al., Laser capture microdissection. Science 274, 998–1001 (1996).
24
EL Greer, A Brunet, FOXO transcription factors at the interface between longevity and tumor suppression. Oncogene 24, 7410–7425 (2005).
25
S Pierrou, M Hellqvist, L Samuelsson, S Enerbäck, P Carlsson, Cloning and characterization of seven human forkhead proteins: Binding site specificity and DNA bending. EMBO J 13, 5002–5012 (1994).
26
Z Xuan, MQ Zhang, From worm to human: Bioinformatics approaches to identify FOXO target genes. Mech Ageing Dev 126, 209–215 (2005).
27
ES Kwon, SD Narasimhan, K Yen, HA Tissenbaum, A new DAF-16 isoform regulates longevity. Nature 466, 498–502 (2010).
28
Y Zhao, et al., Cytosolic FoxO1 is essential for the induction of autophagy and tumour suppressor activity. Nat Cell Biol 12, 665–675 (2010).
29
T Hosaka, et al., Disruption of forkhead transcription factor (FOXO) family members in mice reveals their functional diversification. Proc Natl Acad Sci USA 101, 2975–2980 (2004).
30
A van der Horst, et al., FOXO4 transcriptional activity is regulated by monoubiquitination and USP7/HAUSP. Nat Cell Biol 8, 1064–1073 (2006).
31
H Huang, et al., Skp2 inhibits FOXO1 in tumor suppression through ubiquitin-mediated degradation. Proc Natl Acad Sci USA 102, 1649–1654 (2005).
32
DR Plas, CB Thompson, Akt activation promotes degradation of tuberin and FOXO3a via the proteasome. J Biol Chem 278, 12361–12366 (2003).
33
MA Essers, et al., FOXO transcription factor activation by oxidative stress mediated by the small GTPase Ral and JNK. EMBO J 23, 4802–4812 (2004).
34
A Brunet, et al., Akt promotes cell survival by phosphorylating and inhibiting a Forkhead transcription factor. Cell 96, 857–868 (1999).
35
S Ramaswamy, N Nakamura, I Sansal, L Bergeron, WR Sellers, A novel mechanism of gene regulation and tumor suppression by the transcription factor FKHR. Cancer Cell 2, 81–91 (2002).
36
KE van der Vos, PJ Coffer, FOXO-binding partners: It takes two to tango. Oncogene 27, 2289–2299 (2008).
37
J Debnath, SK Muthuswamy, JS Brugge, Morphogenesis and oncogenesis of MCF-10A mammary epithelial acini grown in three-dimensional basement membrane cultures. Methods 30, 256–268 (2003).
38
KK Ho, SS Myatt, EW Lam, Many forks in the path: Cycling with FoxO. Oncogene 27, 2300–2311 (2008).
39
DN Gross, AP van den Heuvel, MJ Birnbaum, The role of FoxO in the regulation of metabolism. Oncogene 27, 2320–2336 (2008).
40
L Shi, et al., The MicroArray Quality Control (MAQC) project shows inter- and intraplatform reproducibility of gene expression measurements. Nat Biotechnol 24, 1151–1161 (2006).
41
JR Newman, et al., Single-cell proteomic analysis of S. cerevisiae reveals the architecture of biological noise. Nature 441, 840–846 (2006).
42
GJ Kops, et al., Forkhead transcription factor FOXO3a protects quiescent cells from oxidative stress. Nature 419, 316–321 (2002).
43
Z Tothova, et al., FoxOs are critical mediators of hematopoietic stem cell resistance to physiologic oxidative stress. Cell 128, 325–339 (2007).
44
KC Das, Y Lewis-Molock, CW White, Thiol modulation of TNF alpha and IL-1 induced MnSOD gene expression and activation of NF-kappa B. Mol Cell Biochem 148, 45–57 (1995).
45
M Karin, Y Ben-Neriah, Phosphorylation meets ubiquitination: The control of NF-[kappa]B activity. Annu Rev Immunol 18, 621–663 (2000).
46
KA Janes, MB Yaffe, Data-driven modelling of signal-transduction networks. Nat Rev Mol Cell Biol 7, 820–828 (2006).
47
M Bengtsson, A Ståhlberg, P Rorsman, M Kubista, Gene expression profiling in single cells from the pancreatic islets of Langerhans reveals lognormal distribution of mRNA levels. Genome Res 15, 1388–1392 (2005).
48
L Warren, D Bryder, IL Weissman, SR Quake, Transcription factor profiling in individual hematopoietic progenitors by digital RT-PCR. Proc Natl Acad Sci USA 103, 17807–17812 (2006).
49
TS Keshava Prasad, et al., Human Protein Reference Database—2009 update. Nucleic Acids Res 37, D767–D772 (2009).
50
TL Bailey, C Elkan, Fitting a mixture model by expectation maximization to discover motifs in biopolymers. Proc Int Conf Intell Syst Mol Biol 2, 28–36 (1994).
51
QC Lau, et al., RUNX3 is frequently inactivated by dual mechanisms of protein mislocalization and promoter hypermethylation in breast cancer. Cancer Res 66, 6512–6520 (2006).
52
J Pratap, et al., Regulatory roles of Runx2 in metastatic tumor and cancer cell interactions with bone. Cancer Metastasis Rev 25, 589–600 (2006).
53
J Pratap, et al., Ectopic runx2 expression in mammary epithelial cells disrupts formation of normal acini structure: Implications for breast cancer progression. Cancer Res 69, 6807–6814 (2009).
54
P Shore, A role for Runx2 in normal mammary gland and breast cancer bone metastasis. J Cell Biochem 96, 484–489 (2005).
55
T Okuda, J van Deursen, SW Hiebert, G Grosveld, JR Downing, AML1, the target of multiple chromosomal translocations in human leukemia, is essential for normal fetal liver hematopoiesis. Cell 84, 321–330 (1996).
56
FP Silva, et al., Identification of RUNX1/AML1 as a classical tumor suppressor gene. Oncogene 22, 538–547 (2003).
57
H Miyoshi, et al., t(8;21) breakpoints on chromosome 21 in acute myeloid leukemia are clustered within a limited region of a single gene, AML1. Proc Natl Acad Sci USA 88, 10431–10434 (1991).
58
M Kadota, et al., Delineating genetic alterations for tumor progression in the MCF10A series of breast cancer cell lines. PLoS One 5, e9201 (2010).
59
Y Yamamura, WL Lee, K Inoue, H Ida, Y Ito, RUNX3 cooperates with FoxO3a to induce apoptosis in gastric cancer cells. J Biol Chem 281, 5267–5276 (2006).
60
GM Wildey, PH Howe, Runx1 is a co-activator with FOXO3 to mediate transforming growth factor beta (TGFbeta)-induced Bim transcription in hepatic cells. J Biol Chem 284, 20227–20239 (2009).
61
Y Zhang, JR Biggs, AS Kraft, Phorbol ester treatment of K562 cells regulates the transcriptional activity of AML1c through phosphorylation. J Biol Chem 279, 53116–53125 (2004).
62
T Tanaka, et al., The extracellular signal-regulated kinase pathway phosphorylates AML1, an acute myeloid leukemia gene product, and potentially regulates its transactivation ability. Mol Cell Biol 16, 3967–3979 (1996).
63
JR Biggs, LF Peterson, Y Zhang, AS Kraft, DE Zhang, AML1/RUNX1 phosphorylation by cyclin-dependent kinases regulates the degradation of AML1/RUNX1 by the anaphase-promoting complex. Mol Cell Biol 26, 7420–7429 (2006).
64
J Nakae, T Kitamura, DL Silver, D Accili, The forkhead transcription factor Foxo1 (Fkhr) confers insulin sensitivity onto glucose-6-phosphatase expression. J Clin Invest 108, 1359–1367 (2001).
65
ZT Schafer, et al., Antioxidant and oncogene rescue of metabolic defects caused by loss of matrix attachment. Nature 461, 109–113 (2009).
66
T Russo, et al., A p53-independent pathway for activation of WAF1/CIP1 expression following oxidative stress. J Biol Chem 270, 29386–29391 (1995).
67
DA Clopton, P Saltman, Low-level oxidative stress causes cell-cycle specific arrest in cultured cells. Biochem Biophys Res Commun 210, 189–196 (1995).
68
AG Wiese, RE Pacifici, KJ Davies, Transient adaptation of oxidative stress in mammalian cells. Arch Biochem Biophys 318, 231–240 (1995).
69
JH Paik, et al., FoxOs are lineage-restricted redundant tumor suppressors and regulate endothelial cell homeostasis. Cell 128, 309–323 (2007).
70
RM Neve, et al., A collection of breast cancer cell lines for the study of functionally distinct cancer subtypes. Cancer Cell 10, 515–527 (2006).
71
EA Rakha, JS Reis-Filho, IO Ellis, Basal-like breast cancer: A critical review. J Clin Oncol 26, 2568–2581 (2008).
72
A Essaghir, N Dif, CY Marbehant, PJ Coffer, JB Demoulin, The transcription of FOXO genes is stimulated by FOXO3 and repressed by growth factors. J Biol Chem 284, 10334–10342 (2009).
73
B Al-Mubarak, FX Soriano, GE Hardingham, Synaptic NMDAR activity suppresses FOXO1 expression via a cis-acting FOXO binding site: FOXO1 is a FOXO target gene. Channels (Austin) 3, 233–238 (2009).
74
MD Fox, et al., The human brain is intrinsically organized into dynamic, anticorrelated functional networks. Proc Natl Acad Sci USA 102, 9673–9678 (2005).
75
K Blyth, et al., Runx2 and MYC collaborate in lymphoma development by suppressing apoptotic and growth arrest pathways in vivo. Cancer Res 66, 2195–2201 (2006).
76
K Hu, L Ji, KT Applegate, G Danuser, CM Waterman-Storer, Differential transmission of actin motion within focal adhesions. Science 315, 111–115 (2007).
77
JC Wheeler, K Shigesada, JP Gergen, Y Ito, Mechanisms of transcriptional regulation by Runt domain proteins. Semin Cell Dev Biol 11, 369–375 (2000).
78
I Ovcharenko, MA Nobrega, Identifying synonymous regulatory elements in vertebrate genomes. Nucleic Acids Res 33, W403–W407 (2005).
79
YC Lin, et al., A global network of transcription factors, involving E2A, EBF1 and Foxo1, that orchestrates B cell fate. Nat Immunol 11, 635–643 (2010).
80
M Osato, Point mutations in the RUNX1/AML1 gene: Another actor in RUNX leukemia. Oncogene 23, 4284–4296 (2004).
81
N Chapuis, et al., IκB kinase overcomes PI3K/Akt and ERK/MAPK to control FOXO3a activity in acute myeloid leukemia. Blood 116, 4240–4250 (2010).
82
FL Zhou, et al., Involvement of oxidative stress in the relapse of acute myeloid leukemia. J Biol Chem 285, 15010–15015 (2010).
83
M Kato, et al., Superoxide radical generation and Mn- and Cu-Zn superoxide dismutases activities in human leukemic cells. Hematol Oncol 21, 11–16 (2003).
84
ML Guzman, et al., Nuclear factor-kappaB is constitutively activated in primitive human acute myelogenous leukemia cells. Blood 98, 2301–2307 (2001).
85
BP Schneider, et al., Triple-negative breast cancer: Risk factors to potential targets. Clin Cancer Res 14, 8010–8018 (2008).
86
J Fang, T Seki, H Maeda, Therapeutic strategies by modulating oxygen stress in cancer and inflammation. Adv Drug Deliv Rev 61, 290–302 (2009).
87
J Moffat, et al., A lentiviral RNAi library for human and mouse genes applied to an arrayed viral high-content screen. Cell 124, 1283–1298 (2006).
88
K Miller-Jensen, KA Janes, JS Brugge, DA Lauffenburger, Common effector processing mediates cell-specific responses to stimuli. Nature 448, 604–608 (2007).
89
S Gupta, JA Stamatoyannopoulos, TL Bailey, WS Noble, Quantifying similarity between motifs. Genome Biol 8, R24 (2007).
90
DJ Sheskin Handbook of Parametric and Nonparametric Statistical Procedures (Chapman & Hall, 4th Ed, New York, 2007).
91
IN Melnikova, BE Crute, S Wang, NA Speck, Sequence specificity of the core-binding factor. J Virol 67, 2408–2411 (1993).

Information & Authors

Information

Published in

Go to Proceedings of the National Academy of Sciences
Proceedings of the National Academy of Sciences
Vol. 108 | No. 40
October 4, 2011
PubMed: 21873240

Classifications

Submission history

Published online: August 22, 2011
Published in issue: October 4, 2011

Keywords

  1. heterogeneity
  2. triple-negative
  3. stochastic
  4. systems biology

Acknowledgments

We thank Deirdre O'Toole for help with the cloning of the RNA FISH riboprobes, Silvia LaRue for cryosectioning, Laura Selfors for assistance with the tumor sample bioinformatics, and Dr. Richard Iggo for comments on the microarray data from the European Organisation for Research and Treatment of Cancer 10994 clinical trial (GSE6861). This work was supported by National Institutes of Health Grant 5-R01-CA105134 (to J.S.B.), National Institutes of Health Director's New Innovator Award Program 1-DP2-OD006464 (to K.A.J.), the Mary Kay Ash Charitable Foundation (K.A.J.), the Pew Scholars Program in the Biomedical Sciences (K.A.J.), and the David and Lucile Packard Foundation (K.A.J.).

Notes

This article is a PNAS Direct Submission.
See full research article on page E803 of www.pnas.org.
This article is a PNAS Direct Submission. Y.I. is a guest editor invited by the Editorial Board.

Authors

Affiliations

Lixin Wang1
Department of Biomedical Engineering, University of Virginia, Charlottesville, VA 22908; and
Joan S. Brugge
Department of Cell Biology, Harvard Medical School, Boston, MA 02115
Kevin A. Janes2,1 [email protected]
Department of Biomedical Engineering, University of Virginia, Charlottesville, VA 22908; and

Notes

2
To whom correspondence should be addressed. E-mail: [email protected].
Author contributions: L.W., J.S.B., and K.A.J. designed research; L.W. and K.A.J. performed research; L.W., J.S.B., and K.A.J. analyzed data; and J.S.B. and K.A.J. wrote the paper.
1
L.W. and K.A.J. contributed equally to this work.

Competing Interests

The authors declare no conflict of interest.

Metrics & Citations

Metrics

Note: The article usage is presented with a three- to four-day delay and will update daily once available. Due to ths delay, usage data will not appear immediately following publication. Citation information is sourced from Crossref Cited-by service.


Citation statements




Altmetrics

Citations

If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

Cited by

    Loading...

    View Options

    View options

    PDF format

    Download this article as a PDF file

    DOWNLOAD PDF

    Get Access

    Login options

    Check if you have access through your login credentials or your institution to get full access on this article.

    Personal login Institutional Login

    Recommend to a librarian

    Recommend PNAS to a Librarian

    Purchase options

    Purchase this article to access the full text.

    Single Article Purchase

    Intersection of FOXO- and RUNX1-mediated gene expression programs in single breast epithelial cells during morphogenesis and tumor progression
    Proceedings of the National Academy of Sciences
    • Vol. 108
    • No. 40
    • pp. 16483-16857

    Media

    Figures

    Tables

    Other

    Share

    Share

    Share article link

    Share on social media