Effects of Different Tissue Microenvironments on Gene Expression in Breast Cancer Cells

In metastasis, circulating tumor cells penetrate the walls of blood vessels and enter the metastatic target tissue, thereby becoming exposed to novel and relatively unsupportive microenvironments. In the new microenvironments, the tumor cells often remain in a dormant state indefinitely and must adapt before they are able to successfully colonize the tissue. Very little is known about this adaptive process. We studied temporal changes in gene expression when breast cancer cells adapt to survive and grow on brain, bone marrow, and lung tissue maintained in an in vivo culture system, as models of the metastatic colonization of these tissues. We observed the transient activation of genes typically associated with homeostasis and stress during the initial stages of adaptation, followed by the activation of genes that mediate more advanced functions, such as elaboration of cell morphology and cell division, as the cells adapted to thrive in the host tissue microenvironment. We also observed the temporary induction of genes characteristic of the host tissue, which was particularly evident when tumor cells were grown on brain tissue. These early transient gene expression events suggest potential points of therapeutic intervention that are not evident in data from well-established tumors.


Introduction
Signals supplied by the local microenvironment control gene expression during development [1] and when cells are placed in novel tissue microenvironments for investigative purposes [2,3]. Similarly, cancer cells in primary tumors respond to signals present in their local microenvironment [4] and reciprocally alter gene expression in the local stroma [5]. Similar reciprocal effects can be expected between metastatic tumors and their target tissue microenvironments. The individual steps in the progression of primary cancer to metastatic cancer have been reviewed [6,7], and there are extensive gene and protein expression data from preinvasive and invasive primary tumors, metastatic tumors, and circulating tumor cells. The transitions between these physiological compartments, i.e. crossing the basement membrane in the epithelial-to-mesenchymal transition (EMT), entering the circulatory system, extravasating into and colonizing a target tissue, have all been studied using diverse strategies [7].
A common approach used to study extravasation and colonization of target tissue has been comparison of cells serially transplanted to select for those able to invade a specific tissue or organ with different efficiencies. For example, MDA-MB-231 cell lines have been selected to efficiently metastasize to the brain [8]. Gene expression analysis of cells selected in this manner was used to identify genes that mediate breast cancer metastasis to the brain [9]. Similarly, gene expression in 4T1 cells, which metastasize to bone, lung, and lymph nodes, was compared to gene expression in 4TO7 cells, which were derived from the same primary tumor, and which extravasate into but fail to colonize the lung [10], to identify genes that mediate colonization of lung [11]. In those experiments, cells introduced into the circulation were cultured from the target tissue, and this process was repeated, resulting in cell lines with progressively more aggressive ability to extravasate, survive, and colonize the target tissue. However, the efficiency of colonization depends on post-extravasation cell death and on the duration of post-extravasation dormancy [12]. Consequently, correlations between gene expression and metastatic potential for cell lines selected for differential metastatic potential by this approach do not resolve the individual contributions of extravasation, survival, and colonization. Each of these steps in metastasis engages distinct cell-physiological processes and represents conceptually distinct points of potential therapeutic intervention. In addition, once a cell has extravasated and learned to survive, the rate of colonization can increase progressively, and mechanisms underlying the acquisition of more aggressive growth have not been described. Events that ensure the survival of the cell immediately after extravasation are not understood.
We explored the time-dependence of differential gene expression when cancer cells were directly embedded in different tissue microenvironments, as a model for the immediate events that occur during post-extravasation colonization in different metastasis microenvironments. We observed early responses which include the up-regulation of genes normally involved in material homeostasis, whereas later changes include the up-regulation of genes normally involved in cell shape, division and motility, in accordance with gene expression data on established tumors. Of particular interest, we observed that breast cancer cells initially express genes that are characteristic of the surrounding tissue. This was most evident when breast cancer cells were co-cultivated with brain tissue, for which a large number of tissue specific genes can be identified, but may also occur when the co-cultivated tissue is bone marrow or lung. Many of the changes in gene expression induced by co-cultivated tissues and observed in early passages did not persist in later passages. Such transient gene expression must also occur in natural metastasis and cannot be discovered by examining gene expression in natural primary tumors or established natural metastatic tumors, inasmuch as they may precede tumor growth. Therefore, therapeutic targets to delay or prevent the further development of early metastatic lesions might be identified by exploring early gene expression events when tumor cells first encounter a foreign tissue environment.

Results
The processes by which disseminated tumor cells (DTC) adapt to different tissue microenvironments are poorly understood. Typically, disseminated tumor cells experience a dormant phase immediately after extravasation that may persist for many years. These disseminated tumor cells reside in foreign tissue microenvironments as single cells until they acquire the ability to divide. We designed an experiment to identify genes that are differentially regulated during this adaptive period in brain, bone marrow, and lung tissue. We circumvented the extravasation step, itself, by directly grafting tumor cell spheroids onto these tissues prepared using a dorsal skinfold chamber in vivo co-culture system. The dorsal skinfold chamber is a small stainless steel and glass device that is surgically attached to a fold in the skin on the back of a mouse, in such a way as reveal the underside (i.e. the hypodermis) of the skin on one side of the fold through hole in the other side, which is protected by a microscopy window. Minced tissues grafted to the hypodermis in dorsal skinfold chambers survive and revascularize, and retain their original character with respect to several molecular markers [13]. Tumor cells can be grown on these tissue substrates by grafting tumor cell spheroids, yielding a pseudo-orthotopic intravital microscopy model for tumor growth in which growth rate and other parameters can be monitored [14]. The tumor cells can be recovered from the dorsal skinfold chambers and analyzed for gene expression.
Tumor cell adaptation to different tissues Fig.1 presents an outline of the experimental design used in the present study. We prepared brain, bone marrow, and lung tissues in this manner and grafted the murine Her-2/neu-positive transgenic mouse breast cancer cell line, N202, transformed to express H2B-GFP, as described previously [15]. Three weeks after in vivo seeding, cells from the resulting N202 tumors were grown in vitro, in the presence of G418, which selected for the presence of the H2B-GFP cassette. Spheroids were prepared from these cells and were reintroduced into the in vivo culture system for about three weeks. This process was repeated for up to four cycles. Cells grown in this manner on brain tissue for four cycles initially grew more slowly in vitro than did the parental cell line, but eventually returned to or surpassed the initial growth rate. Growth in vitro and in vivo correlated well (Figs. S1A, D, G). In bone marrow, a steady increase in growth rate was observed both in vitro and in vivo (Figs. S1B, E, H), also with strong correlation. Although strong correlation between in vitro and in vivo growth was not anticipated, it supports the conclusion that some characteristics acquired during in vivo adaptation persist during in vitro growth. Adaptation to lung tissue was irregular and correlation between in vivo and in vitro growth was poor. The lack of correlation between growth in vivo vs. in vitro for lung tissue (Figs. S1C, F, I), appeared to be due to growth arrest caused by hemolysis in the IVM chamber at several time points (not shown). Increases in tumor growth rates with serial transplantation have been observed previously, and are generally understood to reflect adaptation to the microenvironment [16,17].
Cancer cells subjected to selection for growth in a foreign tissue microenvironment exhibit changes in gene expression some of which persist when the cells are removed from the foreign tissue microenvironment and grown in vitro. Gene expression profiles for parental cells (P0) and passages 1-4 (P1-P4) were generated from the in vitro cell cultures after cell sorting. Genes were assigned to groups by k-means, based on their expression trajectories, which varied as the cells adapted to the brain, bone marrow, and lung tissue microenvironments over serial passages. Fig. 2 shows graphs of k-means groups for cells grown on brain tissue at passages P0-4 assuming 25 k-means groups. These k-means groups Figure 1. Outline of the experiment. Tissue (brain, bone marrow, or lung) is taken from a donor mouse, minced with a razor, and placed in a dorsal skinfold chamber. After ,7 days, when the tissue has developed vasculature, tumor spheroids prepared initially from the parental (N202) cell line were grafted on top of the vasculaturized tissue bed. After 21 days, the tumor and tissue were placed in cell culture, allowing tumor cells and cells from the tissue to grow. These were selected using G418 and placed back in tissue culture to increase their numbers. These cells were converted into spheroids and reintroduced onto vascularized tissue in a new dorsal skinfold chamber or sorted on the basis fo GFP fluorescence and analyzed for differential gene expression. This was repeated for up to 4 cycles (P0-P4). Genes were grouped by k-means according to differential gene expression in P1-P4 relative to the parental cell line (P0). These k-means groups were then examined for enrichment of Gene Ontology Biological Process Terms. doi:10.1371/journal.pone.0101160.g001 reflect the coordinated regulation of the member genes over successive passages. The extent to which these changes reflect bona fide genetic selection versus epigenetic alterations is unknown, but either mechanism allows for selection of cells in which gene expression promotes or is at least neutral relative to survival and cell proliferation, which form the basis of selection. The corresponding graphs for bone marrow and lung can be found in Figs. S2B and C. The trajectories of individual genes and other sorting tools can be found at http://www.voxvill.org.
In this discussion, k-means groups for brain, bone marrow and lung are designated ''BN N X'', ''BM N X'', and ''LN N X'', where the subscript N refers to the number of k-means groups and X refers to a specific k-means group. There is no rigorous method to define the correct number of groups (see Materials and Methods). For most of the discussion that follows, N = 25 was chosen subjectively to resolve between clearly distinct temporal patterns, but to avoid separating similar temporal patterns on the basis of noise. Other choices for N can be explored at http://www.voxvill. org. GO Biological Process Terms (www.geneontology.org) associated with the genes in each k-means group at statistically significant frequencies were determined and are italicized in the text that follows. GO Biological Process Terms relevant to the discussion and the corresponding GO IDs are listed in Table 1, with p-values from Fisher's exact test [18] corrected to indicate false discovery rates (q-values) [19]. Several physiological and developmental functions commonly associated with the development of tumors, including homeostasis, cell cycle, motility and morphology, angiogenesis, and functions that reflect perturbations in protein dynamics and innate immunity, are apparent in these data and are listed in this table. Following a brief discussion of the induction of these various functions, we show that cancer cells grafted onto brain tissue initially respond to this microenvironment by temporarily elevating the expression of genes normally associated with neurobiological functions. Other significant GO terms, of which there are many, can be explored in detail at http://www.voxvill.org. The following discussion pertains to GO terms that were significant with qval #0.05.

Functional evolution during adaptation
Material Homeostasis The molecular determinants present in a metastatic tumor microenvironment are generally thought to be initially inhospitable to the extravasated disseminated tumor cells, possibly leading to dormancy. Extravasation necessarily implies the immediate proximity of microvasculature, therefore it cannot be taken for granted that extravasated disseminated tumor cells have limited access to nutrients due to a delay in the activation of neovascularization programs. In our experimental model, tumor spheroids consist of only a few tens of thousands of cells grafted onto a fully vascularized bed of tissue, suggesting that diffusion is unlikely to be limiting. Nevertheless, material balance, particularly ion balance, appears to be a key variable. Several k-means groups (Fig. 2) consist of genes that are up-regulated in P1 and P2. The genes in these groups manage early responses of the tumor cells to the novel tissue microenvironment and are often involved in ion transport. To explore this, we determined the GO Biological Process Terms that were disproportionately represented by genes in each k-means group. In these early response groups, the most significant GO terms (qval#0.05) in BN 25 6 ( Fig. 2 and Table 1) are sodium ion transport, L-amino acid transport, and ion transport. Ion transport, anion transport, and ATP hydrolysis coupled proton transport are highly significant in BN 25 15, iron ion homeostasis and arginine transport are significant in BN 25 18, and amine transport is significant in BN 25 19. All of these groups consist of genes that are sharply up-regulated in P1 or P2. The GO annotations of statistical significance associated with these groups can be found in Table  S1, along with a large collection of additional genes characterized by up-regulation in P1 and associated GO Biological Process Terms relating to transport that are statistically significant with qval#0.05. These results indicate that the cells first respond to this novel brain tissue microenvironment by regulating genes involved in transport of materials to maintain homeostasis. The immediate importance of ion homeostasis in adaptation to bone marrow can be seen most easily in the increase in expression of genes in BM 25 4 in which anion transport, sodium ion export, cellular cation homeostasis, elevation of cytosolic calcium ion concentration, positive regulation of potassium ion transport, sodium ion transmembrane transport, and others. In general, adaptation to this microenvironment appears to be less demanding with regard to ion homeostasis, displaying 36 transport-related GO Biological Process Terms compared to 68 such terms for the experiment using brain tissue.
In the lung tissue experiment, there are 23 such terms related to initial up-regulation of gene expression. Although the in vitro growth rate for cells grown on lung tissue did not correlate well with in vivo growth, several of the k-means groups indicate a general, if not smooth, adaptive response to the lung microenvironment, insofar as genes in groups LN 25 2, 7, 8, 20, 23, and 25 undergo net increases over adaptation time, while genes in groups LN 25 1, 10, 16, and LN 25 22 undergo net decreases. Early responses involving ion transport in lung include iron ion transport in LN 25 8, transmembrane transport and sodium ion transport in LN 25 12, and ion transport in LN 25 23 are up-regulated. These observations suggest that the brain tissue microenvironment presents a greater adaptive challenge to ion homeostasis than bone marrow or lung tissue, but that all three respond significantly by activating functions concerned with homeostasis.
Although adaptive response to ion imbalance is an obvious possible explanation for induction of these functions, another possibility is that their induction in the brain tissue experiment is due to a local microenvironment that vigorously controls the expression of genes involved in maintenance of the blood-brain barrier, which tightly regulates ion flow. Consistent with this interpretation is ,7-fold up-regulation of Jam2, ,4-fold upregulation of Cldn1, ,2-fold up-regulation of Cldn3, Cldn12, and Magi2, which are all tight junction components [20], signifying a response to the microenvironment that is not ostensibly homeostatic in nature. Inappropriate regulation of ion flow caused by defective or incomplete microenvironmental regulation of ion channels may actually increase rather than decrease stress, suggesting that only those cells that evolve to ignore these microenvironmental signals may be able to proceed.
Cell growth, morphology and locomotion. In natural disseminated tumor cell dormancy, the tumor cells fail to grow when they first encounter a foreign tissue microenvironment. We observed parallel behavior when N202 cells were transplanted onto brain tissue. The cell-physiological basis for this growth arrest is not well-understood, but a generic response to stress may be responsible. Alternatively, genes exist whose primary selected functions may include prevention of metastasis [21]. Gene expression in a third of the k-means groups (BN 25 2, 4, 5, 16, 20, and 25) is repressed during P1 and P2, during which time cell growth is slow and cell morphology is relatively simple, whereas expression in other groups (BN 25 6,13,15,18,19, and 24) is elevated during this initial phase and anti-correlated with growth rate and morphological complexity (Figs. S1 and S3). BN 25 4, which reaches its lowest point at P2 when cell cycle is at its slowest, contains significant GO terms most directly related to cell cycle ( Table 1). The genes annotated with the GO terms discussed in this and the next section are listed in Table S2. During adaptation, cell shapes became less complex in P1 and P2, and then became more complex in P3 and P4 (Fig. S3), consistent with the regulation of cytoskeletal processes by Rho GTPases [22]. These processes are enriched in BN 25 2 and BN 25 4, coincident with the timing of changes in proliferation and morphological complexity ( Table 1). .32-fold up-regulation places the Kit proto-oncogene from BN 25 19 among the most highly upregulated genes observed in these experiments, although here, the expression of Kit proto-oncogene is inversely correlated with rapid growth.
In the bone marrow experiment, cell-matrix adhesion and positive regulation of monocyte chemotaxis are significant in BM 25 4, in which increases in gene expression correlate well with cell growth rate, as does BM 25 16, which comprises genes associated with GO terms concerning proliferation and mobility including positive regulation of cell proliferation and positive regulation of cell motility. Although these correlations are quite evident and statistically significant in the experiment using bone marrow, both bone marrow and lung tissue appear to be able to immediately support growth of the parental N202 cells with minimal initial adaptation. In lung, growth does not trend steadily higher toward the end of the 4 adaptation cycles, and although many genes in LN 25 14 are associated with the term mitosis, cell cycle and cell division, the first of which achieves the highest significance level of any term in the experiment with lung tissue, this k-means group is distinguished in that is among the most weakly correlated with changes in growth, suggesting that genes associated with the mitotic machinery were left unaffected during this adaptation series. Regulation of cell migration is significant in LN 25 23, in which gene expression increases, overall, during adaptation.

Innate
immunity and endoplasmic reticulum stress. Improper regulation of Ca 2+ and other forms of homeostasis perturbations lead to endoplasmic reticulum stress [23], and there is some indication that the unfolded protein response and innate immunity share regulatory genes [24]. Genes in BN 25 15 and 19 reach their maxima at P2, and are enriched in terms related to endoplasmic reticulum stress and innate immune functions, respectively. In BN 25 15, genes that are induced in P1 and P2 and return to control levels by P3 are enriched in the terms response to endoplasmic reticulum stress, endoplasmic reticulum unfolded protein response, and positive regulation of proteasomal ubiquitin-dependent protein catabolic process. Concordantly, anti-apoptotic Bcl2, which is found in BN 25 6 and protects against endoplasmic reticulum stress [25], is at its highest level in P2, and has been associated with Bcell lymphomas and other cancers [26]. Several genes initially up-regulated in P1 and P2 in group BN 25 19 are associated with innate immune system functions such as opsonization, chemotaxis, positive regulation of macrophage activation, positive regulation of phagocytosis and others. In clinical tumor studies of gene expression, genes related to innate immunity are difficult to interpret correctly due to the infusion of macrophage and other circulating cell types, also observed in nude mouse animal models. In the present experiments, these gene expression measurements are performed on pure cultures of tumor cells sorted on the basis of GFP fluorescence, and therefore cannot reflect contamination by non-tumor cell types. Therefore, these innate immune responses are genuine tumor responses. Inflammatory response and immune response correlate with rapid growth in vitro in cells adapted to grow on lung tissue in LN 25 23. Angiogenesis. Angiogenesis promotes tumor survival and growth, and in this pseudo-orthotopic model, neovasculature arises from pre-existing vasculature in the local stromal tissue [14]. BN 25 20 contains genes that are sharply down-regulated during P1 and P2, and then rebound in P3 and P4. Several of the genes in this group, Plau, Col8a1, and Nrp1, are associated with angiogenesis and Ddah1 promotes endothelial cell proliferation [27]. In BN 25 25, Adamts1 and Sema4a, associated with negative regulation of angiogenesis, are sharply down-regulated. None of these genes is similarly regulated in the bone marrow or lung tissue microenvironments. In BN 25 17, Eya1 and Hey1, regulate sprouting angiogenesis [28], are associated with positive regulation of Notch signaling pathway, and are strongly down-regulated in all three tissues. Nevertheless, revascularization in all three tissues is robust (Figs. S4A,B), possibly indicating that revascularization may be driven by the three tissue microenvironments, rather than by the tumor cells. The tumor cells adapting to brain tissue do provide several pro-angiogenic functions temporarily, e.g. in BN 25 15, Plxnd1, Vangl2, Fzd2, Fzd5, and Vegfa promote patterning of blood vessels and vascular development. In LN 25 20, positive regulation of angiogenesis is a later event that correlates with more rapid growth.
The foregoing discussion in this and in previous sections indicates that tumor cells growing in the pseudo-orthotopic skinfold chamber model modulate gene expression in ways that make sense within the context of our understanding of cancer. The tumor cells must adjust to the immediate stresses of the foreign tissue microenvironment for survival. During this time, external signals received by the cells and cell-autonomous programs adjust to achieve compatibility vis-à-vis appropriate protein synthesis and trafficking. Inadequacies of the local microenvironment with respect to material exchange lead to promotion of angiogenesis. More complex cell physiology such as changes in cell growth and morphology occur only after these more basic problems are solved. These challenges to tumor cell growth appear to be more severe in brain tissue than in bone marrow or lung tissue.

Host tissue genes induced in tumor cells
We infer from these analyses that breast cancer-derived N202 cells transplanted on brain tissue also initially express genes characteristic of the brain tissue microenvironment. For this analysis, we defined GO Neurobiological Process Terms as GO Biological Process Terms containing one or more of the character strings 'axon', 'neuro', 'neural', 'brain', 'neocortex', 'nerv', 'glial', 'hippocampal', 'hippocampus', 'cerebellar', 'cerebral', 'dopa', 'synaptic', 'sensory', 'astrocyte', 'olfactory', and 'memory'. Table 2 lists significant (qval#0.05) GO Neurobiological Process Terms from k-means groups in which genes were initially upregulated in P1 by at least 1.5-fold, using N = 25. In which 34 genes associated with significant GNBP were up-regulated in P1, and 22 were down-regulated. Using various assumptions for N (i.e. N = 9, 16, 25, and 36), 71 genes associated with significant GO Neurobiological Process Terms were up-regulated in brain tissue by at least 1.5-fold in P1, compared with 14 and 20 in P1 in bone marrow and lung, respectively ( Table 3). Occasionally genes annotated with GO Neurobiological Process Terms were assigned to groups in which no enrichment of GO Neurobiological Process Terms was observed. Independently of k-means analysis, enrichment of GO Neurobiological Process Terms in genes up-regulated in P1 by $1.5-fold was observed in the brain tissue experiment relative to non-GO Neurobiological Process Terms (p = 2610 25 , Fisher's exact test) and in P2 (p = 2610 23 , Fisher's exact test), and in the lung tissue experiment in P3 (p = 0.02, Fisher's exact test).
These results indicate strong statistical support for the statement that brain tissue induces GO Neurobiological Process Terms in breast cancer cells in P1 and P2, that is, when the cells first encounter brain tissue. Some enrichment for GO Neurobiological Process Terms also occurs in lung, possibly related to lung innervation.
Comparison of these differential expression data with data from human breast cancer metastases to the brain revealed strikingly high concordance. Among genes included in both this study and in the study by Palmieri et al. [45], 18/21 (86%) of the increases or decreases that accompany the transition from primary to metastatic cancer in clinical samples changed in the same direction in the transition from P0 to P2 in the present data. Permutation analysis over randomly selected groups of these 21 genes showed that this concordance in the direction of change was highly significant (p = 0.0002) (Fig. S5). Therefore, grafting of breast cancer cells onto brain tissue growing in a murine dorsal skinfold chamber recapitulates significant events that occur in human breast cancer brain metastasis.
These results indicate that the N202 cells initially interpret signals from the brain tissue microenvironment by expressing genes that are characteristic of the surrounding host tissue. The cells eventually evolve away from this orthotopic gene expression phenotype and toward a more rapid growth phenotype. In a biological replicate experiment, genes associated with neurobiological functions and induced in P1 and P2 in the first experiment ( Table 1) were similarly induced in the replicate experiment by a ,2:1 margin (27:13 induced vs. repressed in repP1 and 25:15 induced vs. repressed in repP2). Among genes that were repressed in the first experiment, most were repressed in repP1 (5:19, induced vs. repressed), but many of these returned to baseline or above by repP2 (15:9, induced vs. repressed), similar to the later stage, P3 in the first experiment. Genes from BN 25 5 in the replicate experiment were induced immediately, rather than exhibiting the two-cycle delay of the previous experiment. Modeling based on least-squares indicated that repP1 and repP2 corresponded most closely with times P2.60 and P3.68 in the original time series, indicating more rapid adaptation in the second experiment (Fig. S6). Such variation in timing would be expected of adaptation driven by stochastic processes, and as in natural cancer, this comparison emphasizes the ''case study'' nature of each independent adaptive series.
This observation of orthotopic gene expression was facilitated by the cell sorting, which allowed us to identify tissue specific gene expression without interference from contaminating cells from the host tissue, and also by extensive tissue-specific gene expression characteristic of the brain. Functions associated with bone marrow are often characteristic of inflammation as well, and as such cannot be used as evidence of orthotopic gene expression in this model. Nevertheless, their absence would be surprising and, indeed, many genes expressed in bone marrow and associated with inflammation are expressed when N202 cells are placed in contact with bone marrow. For example, in BM 25 1, Tlr2, Ripk2, Myd88 are initially up-regulated and control inflammatory responses [46], such as positive regulation of interleukin-6 production, MyD88-dependent toll-like receptor signaling pathway, and positive regulation of tumor necrosis factor production (see http://www.voxvill.org). Immune response due to Csf3, Cxcl5, and Cxcl1, and inflammatory response due to Itgb6, Cxcl5, and Cxcl1, are associated with P1 of BM 36 25. Several genes that change in P1 have functions in bone unrelated to inflammation, per se: in BM 25 1, Adamts1 is important in inflammation, but also in normal bone physiology [47] and metastasis to bone [48]; Alp1 is a tissue-specific alkaline phosphatase expressed in bone, liver, and kidney [49]; Nr4a2   can be induced in bone-marrow-derived mesenchymal stem cells [50]; Bcl9 is up-regulated in osteoarthritis [51]; Sfrp1, and Tob1 and Srfp1 are significantly associated with GO terms bone trabecula formation and negative regulation of osteoblast differentiation, respectively; Stc2 controls bone growth [52]; and Loxl4 is up-regulated in bone marrow but not brain and lung, and is involved in crosslinking collagen. In BM 36 22, in which genes are up-regulated only in P1, hemopoietic progenitor cell differentiation is significant due to Fst and Sfrp1, as are negative regulation of bone remodeling, negative regulation of osteoclast differentiation, convergent extension involved in somitogenesis and bone trabecula formation. In lung, sterol biosynthetic process achieve very high significance (qval = 4610 25 ) in LN 25 12, an expression profile that tracks closely with in vitro growth. These data show that cancer cells exposed to foreign tissue microenvironments initially respond by expressing genes typical of the cells that comprise the microenvironment. Genes in this early response class represent potential therapeutic targets to delay or prevent brain metastasis using drugs already in use for other purposes, including Alpl, Arg2, Bcl2, C3, Chrm3, Kit, Maoa, and Odc1 (see http://www.voxvill.org), with the essential caveat that the contributions of these genes to the fitness of the metastatic tumor are unknown. A few genes that can be inhibited by known drugs or chemicals are up-regulated at early times and persist or are up-regulated at later times, including Alpl, Plat, and Pnp. At least three drugable genes, Hdac1, Jun, and Vdr, are initially down-regulated by the brain tissue microenvironment, but then return to control levels when more rapid growth is observed, suggesting a delay in the timing of the potential therapeutic window. Interestingly, vorinostat, which inhibits Hdac1, also inhibits brain metastatic colonization in a breast cancer model [53]. In addition to genes that are induced early during adaptation, others induced in brain tissue, bone marrow and lung reach their highest levels in later passages, and these genes may also be of interest vis-à-vis drug targeting, insofar as blocking their later elevated expression may affect tumor cell fitness (see Table  S3). Many other adaptive responses that occur in the serially transplanted tumor cells, including associations with stress and wound healing, innate immune responses, changes in metabolism and others, can be explored further at http://www.voxvill.org.

Discussion
We allowed a murine breast cancer-derived cell line to adapt to various co-cultured tissues in an in vivo culture system, and measured gene expression over time as the cells adapted to the different tissues. We then grouped genes according to similarities in their temporal expression behavior by k-means and correlated the resulting groups with GO Biological Process terminology. This allowed us to describe a sequence of steps taken by a cancerderived cell line during its adaptation to thrive in foreign tissue microenvironments. The sequential activation of functions relating, first, to homeostasis, then to management of endoplasmic reticulum stress, and finally to increased morphological complexity, migration, and cell division parallels the generally assumed model in which metastatic cells that have entered a novel tissue microenvironment by extravasation must first acquire the ability to survive before venturing to proliferate. Less attention has been paid to this early stage of metastatic colonization, in part because early metastatic tumors are microscopic and are therefore difficult to detect, and yield limited material for molecular analysis. In addition, contamination by the surrounding tissue can be difficult to exclude during direct molecular analysis, particularly in syngeneic models.
Breast cancer cells co-cultivated with brain tissue initially express a large number of genes typically expressed brain tissue, presumably due to microenvironmental signals emanating from the brain tissue, similar to stromal effects on gene expression during the epithelial-mesenchymal transition [54]. Induction of orthotopic gene expression was also observed when the cocultivated tissue was bone marrow or lung tissue, although the observation of orthotopic gene expression with brain tissue was more dramatic, possibly due to a richer ontology vocabulary for brain versus bone marrow or lung. Our finding of remarkably high correspondence between gene expression in murine breast cancer cells adapted to grow on minced brain tissue cultivated in vivo and gene expression in clinical breast cancer metastases to the brain suggests that temporal changes observed in our model system are likely to recapitulate temporal phenomena in clinical disease.
We also observed that gene expression evolves during tumor cell adaptation to these various tissues, and that some early changes are temporary. This suggests that analysis of gene expression in well-established clinical tumors is unlikely to accurately reflect earlier gene expression events, some of which may represent novel clinical targets. In our experiments, 89 genes reached their maximum expression of $2-fold in tumor cells immediately upon co-cultivation with brain tissue, and there are drugs and small molecule inhibitors for 26 of these. This suggests potential strategies to prevent metastatic dissemination of breast cancer to the brain, assuming that the higher expression of some of these genes favors survival or growth. Similarly, 205 genes reached a maximum expression $2-fold higher in the second passage on brain tissue, and there are small molecule inhibitors for 54 of these, suggesting a potential strategy to treat slightly more advanced brain-metastatic cancer, still prior to the rapid growth phase. Fiftytwo (52) genes reached their maximum expression of $2-fold in the third passage, and there are inhibitors for 13 of these. Similarly, 201 genes were up-regulated by $2-fold in at least one time point when the breast cancer cells were co-cultivated with bone marrow, 61 of which have known inhibitors. 204 genes were up-regulated by $2-fold in at least one time point when the breast cancer cells were co-cultivated with lung, for which there are also, coincidentally, 61 inhibitors. There were 101 genes up-regulated by $2-fold in both the bone marrow-adapted and lung-adapted cells, and there are inhibitors for 40 of these. Many of these inhibitors correspond to drugs already approved to treat cancer, while many others are currently used to treat maladies other than cancer, and might be repurposed to prevent or treat breast cancer metastasis.
The tissue-dependence of gene expression observed in cells cocultured briefly with various tissues, especially when the changes revert to parental levels in later passages, indicates a potential shortcoming in drug screening strategies for cytostatic and cytotoxic agents that use tumor-derived cell lines growing in vitro that have been maintained in culture for many generations. It would be surprising to discover that the changes induced in vivo were stable for many generations in vitro. In cases in which such induced genes imparted a growth advantage to the cells growing in vivo, one would generally not expect these advantages to persist during the in vitro stage of our protocol. Therefore, drugs that inhibit such genes may influence growth in the tissue environment, but would not be expected to influence growth in vitro, and would be missed in screens that do not have an in vivo component.
The large number of tumor cells grafted onto tissues in the dorsal skinfold chambers in these experiments suggests a potentially significant caveat: the microenvironment experienced by the grafted cells is partially determined by the tumor cells, themselves, whereas a genuine extravasated tumor cell exists as a single cell possibly through a very long period of dormancy. Therefore, by using a large number of cells, we increase the chances of one or more clones surviving and growing, but sacrifice fidelity with the true biological situation. Also, these studies are based on a small number of adaptation experiments, and therefore cannot be generalized without caution. However, the expression of stress-and tissue-specific genes in early stages of adaptation, followed by apparent adaptation to the novel environment, followed by more demanding cell-physiological processes such as morphological reticulation and growth, conform to the prevailing general notions of how cancer, particularly metastatic cancer, must progress. Also, while cells of the experimental tumors are serially transplanted in this experimental design, the organ tissue is not, thereby precluding any long-term adaptive response on the part of the microenvironment in these experiments.
A concerted effort to identify the time-resolved changes that occur during cancer progression in animal models may provide new insights into treatment strategies especially with regard to gene expression changes that occur early during metastasis and therefore represent likely drug targets for metastasis prevention or delay. Our experiments do not exclude the likely possibility that some of the changes in gene expression we observe may occur prior to bona fide metastasis when the tumor cells are in circulation.

Ethics statement
All procedures involving mice were approved by the Institutional Animal Care and Use Committee of Explora Biolabs (San Diego, CA), and performed under Protocol EB09-022. Experiments were carried out in accordance with recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. Every effort was made to minimize animal suffering, anxiety, and discomfort.

Intravital video microscopy
Platinum chambers were fitted in the dorsal skinfold of female nude mice by surgery as described extensively in [13,14]. Brain, bone marrow and lung tissue were obtained from donor mice (C57BL/6), minced with a razor blade and ,10 micrograms of minced tissue were introduced into each dorsal skinfold chamber. These tissues were allowed to revascularize for 7-10 days, to form ''pseudo-brain'', ''pseudo-bone marrow'', and ''pseudo-lung''. Spheroids were prepared as described previously [13,14], from the HER2-positive mouse breast cancer cell line clone N202.1A, which were derived from mammary carcinomas of FVB-NeuN #202 mice (H-2q), expressing a rat HER-2/neu proto-oncogene transgene, controlled by the MMTV promoter [55]. The N202 cells were further modified to express enhanced green fluorescent protein (eGFP) from a stably integrated H2B-GFP construct [15]. These N202 spheroids were grafted on top of the revascularized bone marrow, lung tissue, or brain tissue. Growth of the N202 cells expressing H2B-GFP was followed via the GFP using fluorescence microscopy through the window of the dorsal skin fold chamber. After 21 days, the developing tumors were removed from the chambers and placed in cell culture for 2-4 weeks, depending on the rate of cell growth. Cells expressing GFP were enriched by selection with G418, and again expanded in culture for 2-4 weeks. New spheroids were again prepared and a new cycle of grafting into freshly prepared dorsal skin fold chambers was initiated. This was repeated for a total of 4 cycles. To determine cell growth rates, each passage of N202 for the three tissues were seeded in vitro at ,30000 cells in triplicate and counted every three days during two weeks using a hemocytometer. The mean was calculated for each time point, using 6 or more values.

Gene expression analysis
At each cycle, cells were then sorted based on GFP, RNA was harvested from these highly purified cultures, and gene expression analysis was performed on Affymetrix GeneChip mouse gene 1.0 St array, following the manufacturer's protocol for target preparation, hybridization, washing and scanning. The cells were grown to 60-70% confluence, and total RNA was extracted using an RNeasy Mini Kit (Qiagen, Valencia, CA) following the manufacturer's instructions. RNA concentration was determined by spectroscopy, and integrity was assessed qualitatively by agarose gel electrophoresis. 100 ng of total RNA were used for target preparation.
The parental cell line served as a common control. Technical replicates starting with cells in culture for the common control (P0) and for the cells harvested at the four subsequent passages (P1-4) for cells adapted to brain tissue were performed. Cells harvested from bone marrow and lung tissue were analyzed once per time point. Gene expression data were processed using the RMA (Robust Multiarray Average) [56] algorithm in the Affymetrix Expression Console Software package (Fig. S7). Data were analyzed using routines written in R [57]. Genes were filtered for those having signal intensities exceeding 95% of the signals from negative hybridization controls, and this yielded 10,522 positive probe sets (Fig. S8). Gene expression data are deposited with the NCBI Gene Expression Omnibus (GSE54626).

Data analysis procedures
Ratios were obtained by dividing the signal log ratios of intensities obtained for each gene after each of the four serial passages described above, creating a 5 element vector for each gene of the form (0, x 1 , x 2 , x 3 , x 4 ) where x n are the ratios for each passage, and these vectors were sorted into k-means groups. When using k-means, a decision must be made regarding the number of groups, N, and there is no rigorously correct way to do this. One approach is to set a maximum acceptable fitting error (i.e. the 'sum of squares within'), another is to track the error assuming different numbers of groups and to choose a number of groups beyond which the error no longer decreases rapidly (i.e. identify an ''elbow''), and a third is to try different numbers of groups and test them post hoc for statistically significant biological correlations; we used the latter. We found statistically significant correlations with GO biological process terms using different assumptions for N. Pvalues were determined by Fisher's exact test and corrected to yield the false discovery rate (FDR) according to Benjamini and Hochberg [19]. This was done assuming N = 9, 16, 25, 36, 49, and 64, where N corresponds to the number of k-means groups. Duplicate microarrays were performed for the parental (unadapted) cell line (Fig. S9) and each time point of the brain experiment (Fig. S10), and duplicate signal measurements were averaged after RMA. For the bone marrow and lung experiments, a single microarray experiment was done for each time point.
Each gene was matched to GO Biological Process Terms. Annotations were the union of the original Affymetrix GeneChip mouse gene1.0 St array provided by the Affymetrix Expression Console software and annotations from Mouse Genome Informatics (http://www.informatics.jax.org). For each value of N, the frequency of each GO Biological Process Term within N k-means groups was calculated. Values for the expected frequencies were calculated based on the number of genes in each group. For each GO term and for each group, a 262 contingency table was generated containing the frequency of genes annotated and not annotated with the term within the k-means group, compared with the total number of genes annotated and not annotated with the term. Fisher's exact test was performed on each contingency table to give a list of terms and p-values associated with k-means group. This is a multiple testing situation, so false discovery rates (qval) were calculated by sorting p-values from largest to smallest and dividing by their rank, according to Benjamini and Hochberg [19]. These lists were then sorted by qval from smallest to largest and stored in a MySQL database. A web interface was built to access this database. Genes and GO terms are cross-referenced by links to NCBI/Gene and AmiGO (http://amigo.geneontology. org), respectively. We used Therapeutic Targets Database druggene database (http://bidd.nus.edu.sg/group/ttd/), The Drug Gene Interaction Database (http://dgidb.genome.wustl.edu/), and biodbnet (http://biodbnet.abcc.ncifcrf.gov/); the latter was used to recover official gene symbols for cross-species comparison. We include other drugs and inhibitors gleaned from the literature. These are usually referenced via a Pubmed link. Note that identical or similar drug effects across species cannot be assured.  Figure 2 in the text and is included here for convenience. As in Figure 2, the central black line in each frame represents no change, i.e. 1-fold changes, and the gray lines represent 4-fold increases or decreases, relative to the parental cell line. Horizontally, the five positions along the x-axis represent sequential in vivo passages, P0 (i.e. parental cells), P1, P2, P3, P4, and gene expression analysis was performed on cells expanded in vitro from each of these passages. (B) and (C) represent the corresponding clusters for bone marrow and lung. (TIF) Figure S3 Cell morphology during adaptation to brain tissue. Cell morphology during adaptation to brain tissue varies considerably. Three example fields of cells at each adaptive stage are shown. (P0) The parental cell line has short dendritic projections. (P1) Cells after the first passage become slightly less complex in shape. (P2) After the second passage, cells have very few dendritic projections. (P3) Cells of P3 have numerous, stumpy dendritic extensions, and this phenotype is more pronounced in (P4). (TIF) Figure S4 (A) Robust neovascularization in P1. Bright field images showing that neovascularization of the developing tumor (evident as a raised, bright nodule) becomes robust in the first passage (P1) in all three tissues, brain (bn), bone marrow (bm), and lung (ln). Pressure against the glass due to rapid growth reduces circulation in P1 ln in the picture on the right. (B) Robust neovascularization in all passages. Bright field images of N202 cells growing in three successive passages (P1-3) of independently adapting lineages on brain tissue, illustrating that vascularization is robust for all passages. (TIF) Figure S5 Concordance between clinical brain metastasis and this model system. These histograms show the correspondence between human metastatic breast cancer to the brain and the transition from P0 to P1, 2, 3 and 4 when N202 cells are trained to grow on brain tissue in the dorsal skinfold chamber. This analysis is based on gene expression data from Palmieri et al. [45], in which genes that are up-regulated or down-regulated in breast cancer metastasis to the brain are described. 21 out of 56 genes that were significantly different in metastasis vs. primary tumors identified in [45] were represented by orthologs in our data. We devised a statistic for concordance, c/n, in which c is the number of genes that change in the same direction in the two experiments and n is the number of genes being compared, i.e. n = 21. 18 out of 21 (i.e. c/n = 18/21 = 86%) of these genes that were significantly different in primary tumors vs. metastatic tumors were altered in the same direction in our data in the transition P0RP2, with p = 0.0002. A ratio near 0.86 could occur by chance if ,86% of the changes in gene expression observed in our data corresponded to downregulation; therefore, p-values were determined by comparing the observed c/n (e.g. 18/21<0.86 for P0RP2) with c/n for 10,000 randomly selected sets of 21 genes. These results are summarized as histograms. The p-value of p = 0.0002 (or with a Bonferonni correction; p = 0.0008) for the P0RP2 transition indicates strong concordance of clinical brain-metastatic breast cancer and the P0RP2 transition in our experimental model. This similarity was not seen in the transitions P0RP1, P0RP3, or P0RP4, although the transition P0RP1 nearly reaches a significant p-value. The red vertical lines represent the observed value for c/n and the p-values are indicated. (TIF) Figure S6 Biological replicate of transplantation of N202 tumor cell spheroids to brain tissue. This biological replicate was performed exactly as the initial experiment. There is no reason to anticipate that the rate of evolution in two independent experiments would be identical. Therefore, we sought to match the two sequential passages repP1 and repP2 from the replicate experiment to the appropriate position in the original time course, P0-P4. To achieve this, a spline curve was fit to every gene in the plot of passage vs. RMA-normalized gene expression signal. We used these spline curves to fill a matrix of interpolated intermediate values for each gene at intermediate passages, e.g. P0, P0.002, P0.004,…, P4.000. We then calculated the sum of absolute differences between repP1 and every column in this matrix, and assigned repP1 to the column in the matrix where this statistic was minimum. The position assignment was based on 1000 random bootstrap samples of 100 randomly sampled genes. This was repeated for repP2. Histograms indicated that repP1 and repP2 fit best at P2.65 (blue) and P3.82 (red), respectively. The grey bar corresponds to original P2 as an algorithm check. (TIF) Figure S8 Histogram of PM signals on Affymetrix arrays. Yellow = all PM (perfect match) scores, Red = negative hybridization controls, Blue = positive hybridization ''spike-in'' controls provided by the manufacturer. The vertical red line is placed at the position below which are 95% of the negative controls. Genes were included for which the signal exceeded this threshold for at least one trained or untrained cell line. (TIF) Figure S9 Scatterplot of RMA-normalized untrained control replicates. UT1 and UT2 microarrays were repeated beginning with two independent cell cultures and corresponding RNA preparations, and plotted as a scatterplot. This plot provides a visual demonstration of reproducibility. (TIF) Figure S10 Technical replicate reproducibility. Technical replicates of microarrays for cells trained to grow on brain tissue, beginning with independent cell cultures and corresponding RNA preparations.

(TIF)
Table S1 Significant GO Biological Process Terms related to ''transport'' associated with genes that are UP-REGULATED upon initial contact with Brain, Bone Marrow, and Lung Tissues. (DOCX)