The influence of tumor size and environment on gene expression in commonly used human tumor lines

Background The expression profiles of solid tumor models in rodents have been only minimally studied despite their extensive use to develop anticancer agents. We have applied RNA expression profiling using Affymetrix U95A GeneChips to address fundamental biological questions about human tumor lines. Methods To determine whether gene expression changed significantly as a tumor increased in size, we analyzed samples from two human colon carcinoma lines (Colo205 and HCT-116) at three different sizes (200 mg, 500 mg and 1000 mg). To investigate whether gene expression was influenced by the strain of mouse, tumor samples isolated from C.B-17 SCID and Nu/Nu mice were also compared. Finally, the gene expression differences between tissue culture and in vivo samples were investigated by comparing profiles from lines grown in both environments. Results Multidimensional scaling and analysis of variance demonstrated that the tumor lines were dramatically different from each other and that gene expression remained constant as the tumors increased in size. Statistical analysis revealed that 63 genes were differentially expressed due to the strain of mouse the tumor was grown in but the function of the encoded proteins did not link to any distinct biological pathways. Hierarchical clustering of tissue culture and xenograft samples demonstrated that for each individual tumor line, the in vivo and in vitro profiles were more similar to each other than any other profile. We identified 36 genes with a pattern of high expression in xenograft samples that encoded proteins involved in extracellular matrix, cell surface receptors and transcription factors. An additional 17 genes were identified with a pattern of high expression in tissue culture samples and encoded proteins involved in cell division, cell cycle and RNA production. Conclusions The environment a tumor line is grown in can have a significant effect on gene expression but tumor size has little or no effect for subcutaneously grown solid tumors. Furthermore, an individual tumor line has an RNA expression pattern that clearly defines it from other lines even when grown in different environments. This could be used as a quality control tool for preclinical oncology studies.


Background
Preclinical animal models of human tumors represent a major tool for the selection and development of effective anticancer agents. There is a considerable number of well characterized human cancer cell lines, many of which can be grown as solid tumors (either subcutaneously or orthotopicaly) in immunodeficient mice. The ease of use and low cost of these models make them desirable for screening in vivo activity of anticancer compounds as compared to induced or transgenic rodent tumor models. Rodent models also allow pharmacodynamic and pharmacokinetic parameters to be directly measured and related to antitumor efficacy. Commonly used human cancer cell lines, such as the panel of 60 lines (NCI60) used by the Developmental Therapeutics Program at the National Cancer Institute, have been extensively studied at the molecular level in vitro and to a lesser extent in vivo. Recent advances in genomic technology allow molecular characterization of these models to an extent never before possible using RNA expression profiling, comparative genomic hybridization, proteomic and metabonomic profiling. Ross et al. [1] looked at the in vitro gene expression in the NCI60 panel. The gene expression pattern for many lines indicated a relationship with the tumor tissue of origin and also correlated with doubling time and drug metabolism. Virtanen et al. [2] have analyzed the expression patterns of 85 lung tumor samples from both clinical samples and established tumor lines. The fresh tumors clustered according to pathological subtype with many of the cell lines also clustering within the same groups. Pedersen et al. [3] performed a similar study with human small cell lung cancer cell lines and compared them to ressected tissue samples. Dan et al. [4] recently performed RNA expression profiling on 39 human cancer cell lines and related the gene expression patterns to chemosensitivity of 55 anticancer agents. Zembutsu et al. [5] performed a similar study but profiled 85 human cancer xenografts.
What was apparent in all these studies was that each tumor model was distinctly different and could be distinguished from related models using routine data analysis methods. Furthermore, it appears that it is possible to identify models that are no longer true to their origins. Ross et al. [1] identified MDA-MB-435 as a possible melanoma derived line, which is at odds with its supposed origins of a metastatic breast carcinoma [6]. However, such identification also raises the possibility that the sample that Ross et al. obtained was not the true "MDA-MB-435." Clearly there is the opportunity to use profiling technology as an advanced quality control method that could not only identify mislabeled tumor lines but possibly genetic drift.
While these reports have tremendous value they do not address some basic questions about human tumor mod-els that could impact the design of in vivo drug development studies. For example, most, if not all, xenografts demonstrate Gompertzian growth kinetics with continuously increasing doubling times as they grow larger [7]. Despite this obvious change in biology, it is unknown whether xenografts alter their expression as they increase in size or whether expression differences exist for the same tumor line grown both as a solid tumor in a mouse or in tissue culture. To investigate these questions we grew human colon tumors (HCT-116 and Colo205) in Nu/Nu mice and harvested tumors at three sizes (200 mg, 500 mg, 1000 mg). Isolated RNA was analyzed on Affymetrix GeneChips. The effect of mouse strain on tumor gene expression was also investigated by comparing tumors at 500 mg grown in Nu/Nu and C.B-17 SCID mouse strains. Lastly we profiled additional tumor models to investigate the changes in gene expression that occur when a tumor line is grown in vivo or in vitro.

Xenograft and tissue culture methods
Cell lines were grown in DMEM/F12 supplemented with 10% fetal bovine serum (Invitrogen, Carlsbad, California). They were passaged in 75-cm 2 tissue culture flasks in an atmosphere of 5% CO 2 in air and were subcultured weekly, using 0.05% trypsin EDTA (Invitrogen, Carlsbad, California). RNA from tissue culture samples was harvested at mid-log phase. Human tumor xenografts were grown in either Nu/Nu mice or C.B-17 SCID female mice obtained from Charles River Laboratories (Wilmington, Massachusetts) between four and five weeks of age. Mice were housed five to a cage in animal rooms maintained at between 21-25°C with a 12 h alternating dark/light cycle. All animal studies were conducted under Veterinary Use Protocols approved by the Institutional Animal Care and Use Committee. Tumors were maintained by serial passage of 30 mg tumor fragments between animals, implanted subcutaneously into the right axillary region using a trocar needle aseptically. Tumors were passed when the primary had reached between 500 and 1000 mg and were never passed more than ten times. Tumor growth was followed by caliper measurements of perpendicular measures of the tumor. The weight in mg was estimated by the formula: tumor weight = a(b 2 )/2, where a and b are the tumor length and width respectively in mm. Tumor tissue was harvested immediately following animal sacrifice by excising the tumor and powdering it in a liquid nitrogen cooled crucible and pestle. Tumor powder was stored at -80°C until RNA isolation.

Affymetrix GeneChip
RNA was extracted using TRIZOL ® reagent (Invitrogen, Carlsbad, California) according to the manufacturer's protocol. RNA integrity was monitored using denaturing agarose gel electrophoresis in 1X MOPS. Biotinylated target RNA was prepared from 15 µg of total RNA using the Affymetrix protocol. Briefly, double-stranded cDNA was prepared from the RNA template using a modified oligo-dT primer containing a 5' T7 RNA polymerase promoter sequence and the Superscript Choice System for cDNA Synthesis (Invitrogen, Carlsbad, California). Following phenol-chloroform extraction and ethanol precipitation, one-half of the cDNA reaction (0.5 -1.0 µg) was used as the template in an in vitro transcription reaction containing T7 RNA polymerase, a mixture of unlabeled ATP, CTP, GTP, and UTP, and biotin-11-CTP and biotin-16-UTP (BioArray High Yield Kit, ENZO, Farmingdale, New York). The resulting biotinylated-cRNA "target" was purified on an affinity resin (RNeasy, Qiagen, Valencia, California) and quantified using the convention that 1 O.D. 260 nm corresponds to 40 µg/mL of RNA. Typical yields ranged from 50 -100 µg with transcript sizes between 3.0 to 0.25 kilonucleotides as determined by denaturing gel electrophoresis. Fifteen micrograms of biotinylated cRNA was randomly fragmented to an average size of 50 nucleotides by incubating at 94°C for 35 minutes in 40 mM TRIS-acetate, pH 8.1, 100 mM potassium acetate, and 30 mM magnesium acetate. The fragmented cRNA was hybridized in a solution containing 100 mM MES, 1 M [Na + ], 20 mM EDTA, 0.01% TWEEN 20, 50 pM of Control Oligonucleotide B2, 0.1 mg/mL of sonicated herring sperm DNA, and 0.5 mg/mL BSA for 16 hours at 45°C on either the Human U95A or the Mouse U74A Affymetrix GeneChips ® (Affymetrix, Santa Clara, California). Each hybridization included a mixture of four bacterial biotinylated-RNA transcripts (BioB, BioC, BioD, and cre) spiked at 1.5, 5, 25, and 100 pM, respectively. The hybridization reactions were processed and scanned according to the standard Affymetrix protocols.

Data analysis
All arrays were global scaled to a target intensity value of 600 using the standard Affymetrix protocol. Calculation of the scaling factor, background, noise and percent present, was performed according to Affymetrix protocols using the Data Mining Tool (Affymetrix Santa Clara, California). All resulting data sets were filtered using the absolute call metric (present or absent) using Microsoft Access (Microsoft Corporation, Redmond, Washington). Genes selected had expression levels classified as present at least once in the samples selected for the particular analysis.
To determine the relationship between tumor samples harvested at different sizes the filtered RNA profiling data was analyzed with classic multidimensional scaling (MDS), implemented in R [8,9]. MDS is an unsupervised learning technique that attempts to preserve the relationship between points from high dimensional space at lower dimensional spaces. The program R is a free integrated suite of software facilities for data manipulation, calculation and graphical display http://www.rproject.org.
Analysis of variance (ANOVA) was used to test the effects of tumor size and tumor line, using the following model: Expression of gene i ~ tumor.size + tumor.line + tumor.line*tumor.size. To test the effect due to mouse background, we used the following model: Expression of gene i m ouse background. RNA profiling data filtered on the absolute call metric was used for this analysis and ANOVA was implemented in R. ANOVA is a statistical linear modeling procedure that partitions the total variance into parts corresponding to various sources in the model [10,11]. It has been used previously in microarray data analysis [12][13][14][15].
To rigorously select genes with expression differences between samples, ANOVA p-values were adjusted using multiple comparison procedures. Multiple comparison procedures are tools to adjust p-values that might be inflated as a result of performing multiple hypothesis tests. The Benjamini and Hochberg procedure controls the false discovery rate, which is the expected fraction of false discoveries in all rejected hypothesis [16]. This procedure is less stringent than methods controlling the family wise error rate (e.g. the Bonferroni correction); hence it is more powerful.
Hierarchical clustering was performed in GeneSpring 5.0 (Silicon Genetics, Redwood City, California). The distinction calculation from Spotfire DecisionSite 6.2 (Spotfire Inc. Somerville, Massachusetts) was used to select genes differentially expressed in xenograft samples or tissue culture samples. All data from the tissue culture samples that had an in vivo pair (8 samples) were selected into one group and all data from the xenograft samples (8 samples) were selected into a second group. Genes were prefiltered using the absolute call metric by selecting genes that were present at least once in the selected samples. A distinction value score and p-value was calculated for each gene. The score (≥1) and p-value (≤0.001) was then used to select genes that were differentially expressed between xenograft samples and tissue culture samples. To functionally classify gene lists, web resources such as NCBI (National Center for Biotechnology Information, http:// www.ncbi.nlm.nih.gov/) were searched and the data compiled. Further searching for gene associations in PubMed (NCBI) was also performed.

Variation in tumor xenograft gene expression due to size
We focused on two human colon carcinoma xenografts (HCT-116 and Colo205) to investigate the effects of tumor size and mouse strain on gene expression. Samples were harvested in quadruplicate at three different tumor sizes (200 mg, 500 mg and 1000 mg) for both tumor models grown in Nu/Nu mice (except for the 500 mg sample of Colo205 where five samples were harvested). These sizes were selected as they represent the range at which sensitivity to anticancer agents are traditionally tested and because most models approximate log-linear growth at these sizes. RNA expression profiling data was obtained from Affymetrix U95A GeneChips containing approximately 12600 genes. Genes present (above background) once or more across all samples were selected for further analysis (approximately 7600 genes).
Initial analysis of the expression data with multi-dimensional scaling (MDS) showed that samples from the same tumor line clustered together and that there was clear separation between samples from HCT-116 and Colo205 ( Figure 1). Compared to the profound effect due to tumor Multidimensional scaling plot of Colo205 and HCT-116 samples line, there was no clear separation among samples of different sizes in the MDS plot, suggesting that there was little alteration in gene expression due to differences in tumor size (tumor size effect).
The result from MDS was further confirmed by analysis of variance (ANOVA). Using ANOVA we modeled the effects of tumor line and tumor size on gene expression. Since in the ANOVA we conducted approximately 7000 statistical tests (on the selected genes), with a p-value cutoff of 0.01, we would expect approximately 70 genes (1% of 7000) scored as significantly changed due to chance alone. Indeed, the observed number of significantly changed genes (p ≤ 0.01) due to tumor size effect was 154, approximately twice what was predicted by chance alone ( Figure  2). In contrast, the number of significantly changed genes due to tumor line effect (4731 genes p-value < 0.01) was far greater than chance alone, indicating the two tumor lines were extremely different as suggested by MDS. The distribution of p-values confirmed the profound effect of tumor line on gene expression with genes affected at all expression levels ( Figure 3).
To rigorously identify genes that may suggest functionally significant changes as a tumor increases in size from 200 mg to 1000 mg the ANOVA p-values were adjusted using multiple comparison procedures. Following this analysis none of the previously identified 154 genes had p-values < 0.05. This indicated that there was no change in gene expression as a human tumor xenograft increased in size from 200 mg to 1000 mg. This was true for both the HCT-116 and the Colo205 tumor lines.

Variation in tumor xenograft gene expression due to mouse strain
Additional RNA samples were prepared from 500 mg HCT-116 tumors grown in C.B-17 SCID mice and hybridized to Affymetrix U95A GeneChips. The expression data were compared to data from the same line grown in Nu/ Nu mice harvested at 500 mg. The MDS analysis showed there was no distinct separation or grouping for samples from the two different mouse strains (Figure 1). This suggests that mouse strain played only a small role in altering the expression profiles of HCT-116 tumor samples. However, ANOVA did reveal that considerably more genes Expected verses observed number of significantly changed genes Figure 2 Expected verses observed number of significantly changed genes. The graph shows the number of genes (y-axis) that fall into specific categories (x-axis) based on the ANOVA calculated p-values. The categories are as follows: Predicted, number of genes expected by chance alone; Tumor line effect, number of genes that fall within the specified ranges due to differences in the tumor line; Tumor size effect, number of genes that fall within the specified ranges due to changes in the tumor size; Mouse strain effect, number of genes that fall within the specified range due to changes in the mouse strain the tumors were grown in.
were altered in their expression due to the mouse background (493 genes) compared to chance alone (p < 0.01). To identify genes that may suggest biologically significant differences due to the mouse background effect, the ANOVA p-values were adjusted using multiple comparison procedures. Using the Benjamini and Hockberg procedure [16] 63 genes were found to have a p-value < 0.05 with 32 genes increased in the C.B17-SCID strain and 31 increased in the Nu/Nu strain (Table 1 and 2). Functional classification of these genes using a gene ontology approach did not identify functions that could be linked to known biological pathways. Therefore a biological understanding of the changes in gene expression due to mouse strain remains elusive. P-value distribution/volcano plot of line and size effect Figure 3 P-value distribution/volcano plot of line and size effect. Scatter plot with mean expression level (log 10) on the x-axis and -p-value (log 10) on y-axis. Each gene is plotted twice, once for the p-value resulting from the tumor-line effect ANOVA (black) and once for the p-value resulting from the tumor-size effect ANOVA (red). The blue line represents a p-value of 0.01. Genes with a lower p-value (more significant) have a higher -p-value log10. There are far more genes with significant p-values due to differences in the tumor lines than due to changes in tumor size.

Comparison of lines grown in both tissue culture and solid tumor xenografts
To explore the relatedness of different tumor lines to each other and to investigate the gene expression differences between growth in tissue culture and growth as a subcutaneous solid tumor, RNA samples were prepared from the 13 human tumor lines grown in tissue culture to mid-log phase (Table 3). An additional eight RNA samples were obtained of the same lines grown as xenografts in Nu/Nu mice. The gene expression profiling data resulting from hybridizing the 21 samples to Affymetrix U95A Gene-Chips was filtered for genes classified as present at least once across all samples.
Hierarchical clustering analysis using Spearman rank for samples and Pearson correlation for genes was used to determine the relationship between 13 different human tumor lines (Figure 4). The clustering analysis revealed that for each individual tumor line the xenograft and tissue culture profiles were more similar to each other than any other profile (with one exception). That is, tumor lines clustered together based on their genotypes rather than their growth conditions, suggesting "nature" (genotype) was more influential than "nurture" (growth conditions). This result was confirmed with MDS and principle component analysis (data not shown). The one exception from this pattern was the xenograft sample of ZR-75-1, which did not cluster on the same node as the ZR-75-1 sample grown in tissue culture. However, subsequently we have found that our xenograft ZR-75-1 line has aberrant biological characteristics that are inconsistent with the known phenotypes of this line.
Although each tissue culture-xenograft sample did cluster together there were many gene changes between the paired samples (Table 4). When fold change values were calculated for each tissue culture xenograft pair and the number of genes with ≥2 fold-change tabulated there was a median of 425.5 genes increased in xenografts and 387 decreased. Again the ZR-75-1 sample was aberrant with many more gene changes than the other lines.

Selection and functional assessment of differential expressed genes in xenograft or tissue culture samples
Since genes that are co-expressed in a particular environment may provide information about the adaptive changes to the environment, we identified genes that showed a pattern of high expression in xenograft samples or tissue culture samples using Spotfire's DecisionSite calculation. This analysis identified a small number of genes that matched either pattern. There were 36 genes with increased expression in the xenograft samples relative to the tissue cultures samples and 17 genes that showed high expression in the tissue culture samples relative to the xenograft samples ( Figure 5). Functional classification of these genes showed that they separated into distinct functional groups (see Table 5 and 6). Many of the genes consistently expressed in tissue culture samples encoded proteins involved in cell division, cell cycle, transcription and translation. In contrast, genes expressed in xenograft samples encoded proteins involved in extracellular matrix, cell adhesion, cell surface receptors transcription and translation.

Discussion
The first part of our study was designed to address whether gene expression patterns in xenografts varied with tumor size. No significant changes (p-value < 0.05) were found in any genes after statistical analysis of the RNA profiling data. This was a somewhat surprising result as there is evidence to suggest that tumors show size-dependent biological variation in both a clinical setting and as implanted solid tumors in rodents. It is generally assumed clinically that the larger the tumor mass the lower the likelihood of curing the patient, regardless of the treatment [17,18]. Pathological and genetic heterogeneity increase with increasing tumor cell number and contribute to the emergence of clinical drug resistance [19]. As tumors grow larger their vascular surface decreases, intercapillary distance increases, interstial pressure increases and necrotic foci develop [20]. Tumor doubling times and cell loss also increase with increasing tumor size [17].
Implanted animal tumor models most likely represent clinical end-stage disease and may not necessarily recapitulate clinical behavior seen during tumor development in humans. There are few reports on the molecular and biochemical changes that occur as implanted solid tumors increase in size. Massaad et al. [21] found that some drug metabolizing enzyme systems, including glutathione-Stransferases (GST) did change both activity and expression with increasing tumor size in the mouse colon adenocarcinoma Co38.
Preclinical tumor models show variability in response to treatment with a given anti-tumor agent despite being derived from a single tumor and being implanted into inbred mice [22]. This has been attributed to heterogeneity of the host or the tumor cell populations [23,24] along with tumor-to-tumor variations in perfusion and drug distribution [25]. We found there was no significant difference in gene expression between tumors of either the same size or different size as assessed in four replicate animals. If differences in tumor population, host metabolism or perfusion exist between animals implanted with the same tumors, our data suggests they are not detectable using gene expression analysis of the whole tumor.
Our results suggest that once tumor xenografts have grown to 200 mg, physiological influences on tumor transcription (such as the interaction with host stroma, nutrient supply and oxygenation) have reached steady state and do not change appreciably as the tumor grows to 1000 mg. It is also conceivable that tumors smaller and larger than the sizes investigated in this study may show gene expression changes. When a tumor is first implanted the hypoxic environment is likely to drive the expression of pro-angiogenic cytokines that will recruit new vessels to the tumor. As a tumor becomes very large (>2 g) angiogenic factors may again become highly expressed due to the inability for the existing vasculature to adequately support the large tumor mass. We have not excluded the possibility that there may be changes in protein activities For any two samples, the vertical distance from the sample roots to the first node joining them is a measure of their similarity; the shorter the distance the more similar. The color in each cell of the table represents the median adjusted expression value of each gene. The color scale used to represent the expression ratios is shown on the right, with yellow indicating increased expression relative to the median and blue decreased.  Heat-map plot of differential expressed genes in xenograft or tissue culture samples By comparing the same tumor line grown in different mouse strains we identified 63 genes with a significant change in gene expression with biological functions ranging from cellular signaling, RNA processing and translation. Although different genes were differentially expressed in different mouse strains, particular gene functions did not associate with one strain. Therefore, the biological significance of these changes remains unclear. It should be noted that although the changes were statistically significant, in general the magnitude of the change was small. Only 12% of the significant genes displayed a change greater or equal to 2-fold. This raises the question of whether statistical significance equates to biological significance and how relevant thresholds can be determined for analyzing expression profiling data.

Hierarchical clustering analysis of 13 tumor lines grown in different environments
In the second part of the study we compared the gene expression patterns of tumor lines grown as xenografts with the same lines grown in tissue culture. We found that each xenograft -tissue culture pair was more similar to one another than any other line. This has been previously observed with RNA expression profiling of small cell lung cancer lines grown as both xenografts and tissue culture [3]. In our data, the one exception to this clustering pattern, ZR-75-1, exhibited aberrant growth characteristics; it had a rapid doubling time (2-3 days) and was not estrogen dependent as has been reported [26,27].
There was a gene expression pattern consistent with growth either in a xenograft or tissue culture environment, but remarkably it consisted of only a small number of genes. We identified 36 genes with high relative expression in xenograft samples and 17 genes with high relative expression in tissue culture samples. Their biological functions were consistent with the differences between a tissue culture sample and a xenograft sample. Genes expressed at increased levels in tissue culture samples were primarily involved in cell division, cell cycle, transcription and translation and were consistent with a greater percentage of cycling cells in the tissue culture samples, which we have previously observed with flow cytometry (data not shown). Many of the genes expressed in xenograft samples encoded for proteins involved in extracellular matrix, cell adhesion, cell surface receptors transcription and translation and suggest increased interactions with an in vivo stromal environment, including the development of a 3dimensional matrix.
There are several reports in the literature of tumor line expression profiling studies demonstrating that lines derived from a common tissue of origin group together in most cases [1, 4,5]. We observed that most of the breast tumor lines clustered together but the other lines generally did not show clustering patterns based upon tissue of origin. However, our study evaluated fewer samples and consisted predominantly of breast cancer tumor lines with only two representatives of non-small cell lung cancer and colon cancer.
Since each tumor line was substantially different to every other line tesed it suggests that it would be possible to identify a group of genes that could be used to distinguish individual lines. Also, it has been our experience that it is possible to identify aberrant samples resulting from labeling error by expression profiling that were missed by other methods. What is less clear is whether genetic drift causes changes in RNA expression that can be identified using RNA profiling technology. Recently a report has been published showing that MCF-7 sublines were substantially different to each other both at the genetic and RNA expression levels [28] suggesting that variation within individual lines can be identified. Given this, it seems reasonable that RNA expression profiling could be used as a comprehensive methodology for identifying aberrant or incorrectly labeled samples. This would provide an additional quality control tool to standardize tumor models and the in vivo testing of therapeutic agents.

Conclusions
Our data suggest that the environment a tumor line is grown in can have a significant effect on gene expression even though tumor size has little or no effect for a subcutaneously grown solid tumor. Furthermore, an individual tumor line has an RNA expression pattern that clearly defines it from other lines even when grown in different environments such as tissue culture or in vivo. Routine RNA expression profiling of a selected set of genes could be used as a quality control tool for preclinical oncology