A transcriptional profile of multicellular development in Dictyostelium discoideum

The soil amoeba Dictyostelium discoideum has a relatively simple developmental program. It is characterized by a series of highly coordinated cellular, physiological and morphological changes and serves as a paradigm in the research of multicellular development (Kessin, 2001; Loomis, 1975). Upon starvation, the unicellular amoebae stop dividing, aggregate, form a multicellular organism that consists of two cell types and eventually construct a fruiting body in which a ball of spores is carried aloft on a cellular stalk (Fig. 1).

Dictyostelium development is a series of synchronous, coordinated morphological and physiological changes. The first is a transition from growth to development that is induced by starvation (Clarke and Gomer, 1995), but is not accompanied by macroscopic morphological changes (Fig. 1, 0 hours). The second change occurs when the amoebae begin to aggregate and communicate through secretion of cAMP (Fig. 1, 6 hours) (Parent and Devreotes, 1996). 2-4 hours later, cell motility and chemotaxis to cAMP mediate the aggregation of groups of up to 100,000 cells into loose mounds (Fig. 1, 10 hours). Overt cell-type divergence follows aggregation as the amoebae differentiate into prespore and prestalk cells that are distinguishable by molecular markers (Takeuchi, 1991; Williams et al., 1989). Later, an extracellular matrix is secreted and the cells become enveloped in an acellular sheath (Fig. 1, 14 hours). In the following 6-8 hours, the multicellular organism undergoes a series of morphological changes, including a remarkable transition into a slug-shaped structure that migrates towards light and heat (Fig. 1, 18 hours) (Miura and Siegert, 2000; Raper, 1940). This feature illustrates that the slug is a bone fide multicellular organism, capable of sensing its environment and responding by coordinated movement. The last dramatic morphological transition begins after 18-20 hours of development (Loomis, 1975). The prestalk cells undergo terminal differentiation and form a multicellular stalk while the prespore cells encapsulate and become dormant spores. The ultimate structure consists of a stalk that carries a ball of spores about 1 mm away from the substratum (Fig. 1, 24 hours). There are numerous studies that correlate these developmental transitions with changes in gene expression and the number of developmentally regulated genes was estimated as between 300 and 1000 (Alton and Lodish, 1977; Firtel, 1972; Loomis, 1978; Morrissey et al., 1984). We were interested to know whether all of the morphological changes coincide with discernable physiological changes and which developmental events are accompanied by the largest physiological change.

A general way to approach these questions is to apply genomic methods to experimentally tractable developmental systems. One such approach is transcriptional profiling with microarrays. Recent reports on the transcriptional profiles of development in C. elegans and in Drosophila have concentrated on the identity and possible function of the developmentally regulated genes and did not specifically address the above questions (Furlong et al., 2001; Hill et al., 2000; Kim et al., 2001; White et al., 1999). In addition, the complexity of these developmental systems may have presented a challenge when trying to correlate developmental events with specific transcriptional profiles. The relative simplicity of the Dictyostelium multicellular phase and the fact that a large number of cells can be induced to develop with a high degree of synchrony make Dictyostelium an attractive model for the study of genomic control systems in development. At the same time, the molecular and cellular events that lead to terminal differentiation in Dictyostelium are very similar to those described in the development of metazoa, including cell motility and sorting, cell-cell adhesion, morphogenesis, intercellular signaling, signal transduction and coordinated regulation of gene expression (Chung and Firtel, 2000; Ma et al., 2001; Parent and Devreotes, 1999). Consequently, studies in other systems should benefit from the comprehensive studies that are possible in this organism.

There are two complementary approaches to the analysis of expression array data. In one approach, genes are grouped according to their expression pattern under the premise that coregulated genes are likely to share functional characteristics. This approach has been useful in the analysis of the cell cycles of human fibroblasts, yeast and bacteria (Cho et al., 2001; Iyer et al., 1999; Laub et al., 2000; Spellman et al., 1998). However, an observation was made in yeast that raises questions about the general applicability of this approach (Winzeler et al., 1999). In that study, the authors made a systematic comparison between gene function as measured by the fitness of null-mutant strains under certain conditions and gene expression under the same conditions as measured with an expression array. Surprisingly, there was no obvious correlation between gene function and gene expression (Winzeler et al., 1999). This study shows that assigning function to genes based on their pattern of expression alone is not applicable in every experimental system.

The other approach in array data analysis is to consider the expression pattern as a reflection of cell physiology. This approach has been applied successfully to the classification of cancer cells, where the pattern of gene expression provided enough detail to differentiate between tumors that were nearly indistinguishable by other means (Alizadeh et al., 2000; Bittner et al., 2000; Golub et al., 1999). This type of analysis was extended into a comparative approach in yeast (Hughes et al., 2000). In that study, a similarity between the transcriptional profiles of two mutant yeast strains was shown to be sufficient to indicate functional similarities between the mutated genes. Similarly, the transcriptional profile of drug-treated cells was similar to that of mutant cells in which the drug target gene was deleted (Hughes et al., 2000). We applied this type of analysis to the investigation of the physiological changes which occur during Dictyostelium development.

We compared the transcription profiles of cells at 2-hour intervals throughout the 24-hour period of development and found that the largest transition in cell physiology occurs at the same time as the morphological transition from unicellular development to multicellular development. We also found that about a quarter of the genes in the genome were regulated during that period and that this regulation occurred with little or no variation between cells from different strains or different nutritional histories. Comparison of the transcriptional profiles of cell-type-enriched genes revealed unexpected patterns of gene expression that suggest the use of common physiological modules in different developmental stages.

A collection of 7385 clones was used: 5655 cDNA clones from the Dictyostelium cDNA project (Morio et al., 1998) (http://www.csm.biol.tsukuba.ac.jp/cDNAproject.html); 987 cDNA clones were selected from a low redundancy screen of a lambda library (Shaulsky et al., 1995) and a plasmid library (S. Lu and A. Kuspa, unpublished) of cDNA from late developmental stages and from vegetative and early developmental stages of AX4 cells, respectively; 647 genomic DNA clones from the Dictyostelium Genome Project at Baylor College of Medicine (http://dictygenome.bcm.tmc.edu/) were selected as long open reading frames that matched published protein sequences; and 96 clones were from miscellaneous sources. All the clones were sequenced to verify their identity. Sequences were compared to public databases and annotated (see supplemental data). The degree of redundancy is estimated at less than 20%, based on sequence analysis. The array also contained 198 control targets that were made from the Dictyostelium ribosomal 17S RNA gene, histone H1, actin8 and mhcK, as well as control targets from yeast genes and ‘no DNA’ controls. Altogether, the array contained 7744 targets. The entire array was printed in duplicate.

Wild-type Dictyostelium discoideum strains AX2 (Watts and Ashworth, 1970) and AX4 (Knecht et al., 1986) were grown in association with Klebsiella aerogenes on SM plates or in HL-5 medium (Sussman, 1987). Exponentially growing cells were washed from the nutrient source, deposited on nitrocellulose filters and developed in the dark at 22°C (Shaulsky and Loomis, 1993). At each time point, 1×10⁸ cells were collected, resuspended in 1 ml Trizol reagent (Life Technologies) and total RNA extracted according to the manufacturer’s protocol.

AX2 cells were grown in association with K. aerogenes bacteria. Vegetative cells were harvested, resuspended in 20 mM K₂PO₄, pH 6.4, and plated on 1.5% agar plates at a density of 5×10⁶ cells/cm². The plates were incubated at 4°C for 4 hours and then at 22°C for 15-18 hours. Slugs were collected on a 77 μm nylon membrane, resuspended in 0.1% (w/v) Pronase and 0.1% BAL in 50 mM Tris, pH 6.8 (Takeuchi and Yabuno, 1970), disaggregated by trituration, washed with 20 mM K₂PO₄, 20 mM EDTA, pH 6.5, and resuspended at 1×10⁹ cells/ml in the wash buffer. Disaggregated cells were loaded on top of a solution of 50% Percoll in 5 mM MES, 20 mM EDTA, pH 6.8, and centrifuged at 27,000 g for 15 minutes. The two cell types accumulated in two bands; the upper band contained prestalk cells and the lower band contained prespore cells. The bands were collected separately, washed and separated again by centrifugation through Percoll as above (Ozaki et al., 1988). RNA was extracted from the separate cell populations as above and purity was verified by northern blot analysis (data not shown).

Spores were collected from 2- to 4-hour-old fruiting bodies, washed once in 20 mM K₂PO₄, pH 6.4, and centrifuged for 5 minutes at 340 g in a clinical centrifuge. The pellet was resuspended in 20 mM K₂PO₄, 10 mM EDTA, 0.1% NP40, pH 6.5, and triturated through an 18-G syringe needle to eliminate unencapsulated cells. Spores were centrifuged for 7 minutes at 340 g in a clinical centrifuge, washed with 20 mM K₂PO₄, pH 6.5, and ground in a mortar and pestle while in liquid nitrogen until 90% or more of the spores appeared broken by microscopic examination. The spore extracts were collected into Trizol and RNA was extracted as above.

Stalks were collected, resuspended in 20 mM K₂PO₄, pH 6.4, and filtered 5 times through a 77 μm nylon membrane to eliminate spores. Stalks were ground in liquid nitrogen and total RNA was extracted as above. In the cell-type enriched experiments, RNA samples were from 2-4 different preparations and the data were averaged.

DNA targets were amplified from plasmids by polymerase chain reaction (PCR) with common oligonucleotides and their size verified by gel electrophoresis. PCR products were purified by precipitation with 50% isopropyl alcohol, 0.3 M sodium acetate, pH 5.2, washed once with 70% ethanol, dissolved in water and adjusted to 800 mM NaCl, 200 mM Na₂PO₄, pH 10.5.

Glass slides (Gold Seal Products, VWR) were washed by sonication in acetone for 10 minutes, rinsed twice in distilled water, immersed in 0.1 N NaOH for 10 minutes, washed in distilled water, and immersed for 3 minutes in 6% (v/v) 3-glycidoxypropyltrimethoxysilane (Aldrich) made in 95% ethanol/acetic acid, pH 5.0. The slides were washed in 100% ethanol, dried with nitrogen gas and baked for at least 2 hours at 120°C.

Target DNA was printed on the activated glass slides on 200 μm centers with a Cartesian Pixsys5500 robot using Chipmaker II pins (Tele-chem International). Arrays were stored desiccated in the dark.

The DNA primers for cDNA production from total RNA were deoxythymidine octadecamers [dT(18)]. The 5′ terminal nucleotide was modified with Cy3 or with Cy5. All the primers were HPLC-purified by the manufacturer (Operon Technologies). Total RNA (10 μg) was mixed with 0.5 μg of labeled [dT(18)] primer in 13 μl of water, incubated at 70°C for 10 minutes and on ice for 2-5 minutes. Reaction buffer (Gibco-BRL), 0.1M DTT, 0.5 mM of each dNTP and 200U Superscript II (Gibco-BRL) were added and the reaction incubated at 42°C for 2 hours. Reverse transcription reactions were terminated with 0.1 M EDTA. RNA was degraded by adding 0.3 N NaOH and incubating at 60°C for 20 minutes. The reaction was neutralized with 0.4 M Tris-HCl, pH 7.6. Labeled cDNA was precipitated at room temperature with 0.3 M sodium acetate, pH 5.2, and 2 volumes of ethanol, washed once with 70% ethanol and resuspended in 5 μl of water.

For each microarray experiment, an experimental RNA sample was compared to a reference RNA sample that was a pool of equal portions of RNA samples prepared from several developmental stages (0, 3, 6, 12, 17, 24 hours). The experimental sample was labeled with Cy5, the reference sample with Cy3. Both probes were combined and 130 μl of Perfecthyb^TM Plus Hybridization buffer (Sigma) was added. The solution was boiled for 2 minutes, cooled to 65°C and applied to the array. A GeneTAC^TM hybridization station (Genomic Solutions) was used for hybridization for 2 hours at 65°C, followed by three washes in 5× SSC + 0.1% SDS, 2× SSC + 0.1% SDS and 0.1× SSC.

The arrays were scanned with a Scanarray5000 scanner (GSI Lumonics) and images were processed with the GLEAMS software package (NuTec Sciences, Inc.). Quantified data were passed through a single-chip normalization procedure to correct for spatial artifacts, to estimate the variability of replicate log-ratios and to bring the data to a common measurement scale to allow for multi-array comparisons (see supplemental data). The normalization was implemented in seven distinct steps: thresholding, spatial adjustment, averaging of on-chip replicates, by-signal-size variance estimation, identification of outliers, by-signal-size log-ratio adjustment (Callow et al., 2000) and scaling of the final values using the estimated by-signal-size variance. The processed data values were the quantities assessed in all multi-array analyses. One of the most important points to notice is that the expression level values are reported as the gene expression level relative to the average gene level across time.

To identify genes whose expression levels were altered dramatically during the transition from unicellular to multicellular development, we used the function (y=x/12–1). This function can be thought of as describing the expression of a hypothetical gene where x is the developmental time in hours and y is the normalized gene expression level. In the range 0-24 hours, the expression levels of this hypothetical gene increase from –1 to +1. Every gene trajectory in the experimental data was compared to the function. Genes that fit the function well are expressed at a low level in early development and at a high level later in development. These genes receive a high positive score. Genes that are expressed at a high level early in development and at a low level later on receive a high negative score. Genes whose expression is not altered dramatically during the transition receive a score close to 0. Therefore, the absolute value of the score that fits the gene trajectory to the function (y=x/12–1) determines whether the gene is regulated during the transition from unicellular to multicellular development, and the real value of that score can be used to sort the genes into informative groups and to order them within each group. Additional details are provided in the supplemental data.

To assess the reproducibility of the experimental system, AX4 cells were grown on HL5 in three separate experiments and developed separately. RNA samples were collected and analyzed in duplicate from each sample as described above. The function (y=x/12–1) was used to calculate scores for each gene in each experiment at each time point. The values within a time point were averaged within each biological replication. The mean squared error of each gene was calculated as the average (across time points) of the within-time-point variance of the gene. These values were used to calculate a by-gene T-statistic for the gene’s trajectory. Finally, the scores were compared to the scores obtained in the four strain/nutrition experiments to determine the number of genes whose trajectories were reproducible (α=0.05).

A different assessment of reproducibility compared the relative contribution of the biological and the technical variations. We define biological variation as that observed in the three independent RNA preparations and the technical variation as that observed between duplicate array experiments from each one of the RNA preparations. The analysis was performed by fitting a 2-way ANOVA model to each gene, with one factor being time and the other factor being biological preparation. Details are given in the supplement data.

Data from the cell-type enriched samples were subjected to by-gene F tests to identify genes that showed a differential behavior in at least one cell type (α=0.05). Genes that showed significant F tests were further subjected to tests against five linear contrasts, representing particular patterns of cell-type enriched expression. For example, the linear contrast that describes spore-enriched genes is (3, –1, –1, –1) where the positions describe RNA extracted from spores, prespore cells, prestalk cells and stalks, respectively. Experimental data from each gene are multiplied by the respective coefficients and the results are added to give a score. A gene that is expressed at a high level in spores and at a low level in all other cell types would receive a high score and vice versa.

The genes with significant cell-type contrasts were then examined in the time-course data. Two linear contrasts [(–1, 1.6, 1.6, 1.6, –1, –1, –1, –1, –1, –1, –1, 1.6, 1.6) and the function (y=x/12–1)] were used to order the time course data for these cell-type enriched genes as described above. Detailed description of the analysis is provided in the supplemental data.

In order to monitor the transcriptional changes that occur during Dictyostelium development we assembled a DNA microarray that consists of hybridization targets for about 75% of all genes (for details see Materials and Methods and supplemental data). Using these microarrays, we monitored the relative abundance of mRNA in samples that were collected throughout the course of development.

The most distinct morphological states in development are shown in Fig. 1. We wanted to determine whether all of these morphological states are associated with distinct transcriptional profiles and which developmental transition was accompanied by the largest transcriptional change. We therefore developed cells, collected RNA at 2-hour intervals and examined the relationship between the 13 time-point samples using the expression array. The data from each time point were compared to all the other data by calculating the dissimilarity between the expression levels for each of the 7385 genes represented on the microarray. The results are summarized in the form of a dendrogram in Fig. 2A. Each leaf in the dendrogram represents a time-point RNA sample and the height of each join is directly proportional to the dissimilarity between the joined leaves or lower joins. For example, the most similar RNA samples are the ones collected at 0 hours and at 2 hours, as the join that connects them is the lowest. The largest dissimilarity between the time-point samples was found between the 0-6 hour group (red) and the 8-24 hour group (green) (Fig. 2A). This dissimilarity coincides with the transition from unicellular development to multicellular development between 6 and 8 hours, which is one of the most dramatic morphological transitions in Dictyostelium development (Fig. 1). The next large dissimilarity in the data distinguishes the 8-10 hour group from the later time points. The 8-10 hour period is not accompanied by gross morphological changes, but coincides with the beginning of cell-type divergence. These large dissimilarities in the expression data implicate the transition to multicellular development and the initiation of cell type divergence as the two developmental processes that require the largest change in cell physiology. On the other hand, the pre-aggregation stage (0-6 hours) and the finger stage (14-18 hours) coincide with the most coherent groups of transcription profiles in the experiment. These stages are also characterized by the smallest changes in gross morphology.

It was also interesting to find that the dendrogram conserved the temporal relationships between the samples even though temporal information was not included in the dissimilarity calculation. This result serves as an important control since the subtle asynchrony between the organisms in our samples (±1 hour) demands that adjacent time points be most similar. Therefore, conserving the temporal order of the time points indicates that the large-scale transcriptional changes that we observed are a true reflection of cellular physiology.

Dictyostelium cells develop with nearly invariant timing and morphology (Loomis, 1975). We hypothesized that if the invariant developmental timing and morphology result from a common and robust developmental program, it should be accompanied by an invariant pattern of gene expression. To test this idea, we analyzed the transcriptional pattern during the development of two common laboratory strains, AX2 and AX4, which were grown on two different nutritional sources, bacteria or nutrient broth (axenic growth). Vegetative cells growing on bacteria are different from axenically grown cells. The doubling time of bacterially grown cells is 3 hours and that of axenically grown cells is 8 hours. The cell volume and the protein, RNA, DNA and carbohydrate content of bacterially grown cells vary by 1.5- to 5-fold from those measured in axenically grown cells, indicating that the cells are physiologically different (Ashworth and Watts, 1970; Leach and Ashworth, 1972). However, when starved, the cells develop with nearly indistinguishable morphology and timing, regardless of their nutritional history (Loomis, 1975). We used the different nutritional conditions to test the robustness of the developmental transcriptional profile. Cells from the two different strains were grown axenically or in association with bacteria and developed. RNA was extracted at 2-hour intervals and tested with the expression array. The normalized gene expression levels were compared across the four strain and growth conditions. About 4000 genes were found that exhibited a consistent pattern of expression in at least two of the four strain/nutrition conditions and about 3000 of them were regulated at some point during the 24-hour course of development (data not shown). The dissimilarity calculation shown in Fig. 2A indicated that the expression of many genes must be altered during the transition from unicellular to multicellular development. We therefore selected 2021 genes that exhibited a prominent transcriptional change during that transition in two or more of the experiments and defined them as the consensus group of developmentally regulated genes in Dictyostelium (Fig. 2B). These findings indicate that approximately 25% of the estimated 8,000 genes in the genome are regulated during the transition from unicellular to multicellular development.

There are two prominent trends in the data shown in Fig. 2B. The top 711 rows represent genes that display lower-than-average expression during growth and early development and higher-than-average expression in later times. This pattern of gene expression is characteristic of many developmentally induced genes such as the spore-coat genes cotA, cotB and cotC (Fosnaugh and Loomis, 1991; Haberstroh et al., 1991) that were found among the top 711 genes in Fig. 2B. The bottom 1310 rows in Fig. 2B represent genes that are expressed at a higher-than-average level during growth and early development and at a lower-than-average level later in development. Genes that were previously described to have this pattern, e.g. the cysteine protease gene cprD (Souza et al., 1995) and a group of vegetative ribosomal genes (Singleton et al., 1989) were found in this cluster (Fig. 2B). The correlation with published reports of developmental gene regulation indicates that our expression array results provide an authentic measure of gene expression during development.

The 2021 genes that we defined as the consensus group of developmentally regulated genes were selected because they were expressed in a similar manner in at least two out of the four strain/growth experiments. However, this definition could have been skewed if one of the two conditions, strain or nutritional history, had resulted in more dramatic gene regulation profiles than the other condition. To test that possibility, we compared the four individual strain/nutrition data sets as follows: the genes in each data subset were ordered by fitting to the function (y=x/12–1) and the resulting orders were imposed on the two most different data subsets, axenically grown AX4 cells and bacterially grown AX2 cells (Fig. 3). We found that regardless of the experimental condition, the set of 2021 genes conserved its expression pattern and was not sensitive to the strain or the nutritional history. These findings support the notion that Dictyostelium development is a robust process, which is accompanied by an invariant physiological process as reflected in its transcriptional profile.

To test the reproducibility of the experimental system we repeated the developmental time course analysis of AX4 cells grown in liquid broth three independent times. The data were used to assess the biological reproducibility by calculating the variation between the three different experiments and the technical reproducibility by calculating the variation between the two replicate arrays within each biological sample. We found that out of the 2021 genes that exhibited a conserved expression pattern across the four different strain/nutrition experiments, only 84 genes (about 4%) showed an altered trajectory relative to their consensus pattern. This finding illustrates that the expression trajectories are highly reproducible. Comparing the technical variation within each sample to the biological variation within each time point we found that about 20% of all the genes on the array exhibited a biological variation that was equal to or greater than the technical variation and that number was much smaller (9%) in the consensus group of 2021 genes. This finding indicates that the largest component of the variability can be attributed to technical sources.

Following the transition from unicellular to multicellular development, Dictyostelium cells differentiate into two types of cells, prestalk and prespore, which express specific genes. The prestalk marker gene ecmA was first described as a gene whose expression was induced by the stalk differentiation-inducing factor DIF-1 (Williams et al., 1987). The gene is expressed exclusively in prestalk cells after about 12 hours of development and it encodes an extracellular-matrix protein. Expression of ecmA is controlled by a number of factors, including STAT proteins, and its regulatory element contains sub-domains that have been useful in defining various prestalk cell subpopulations (Williams, 1997). Prespore cells are defined by the expression of pspA and the coordinately regulated spore coat genes cotA, B and C (Early et al., 1988; Fosnaugh and Loomis, 1991; Haberstroh et al., 1991). To begin exploring cell-type divergence from a global transcriptional perspective we performed expression array analyses on RNA from separated prespore and prestalk cells and from separated spores and stalks. We found 873 targets that showed a strong cell-type preference. Of those, 328 targets were enriched in spores, 335 were enriched in prespore cells or in spores, 150 were enriched in prestalk cells and 60 were enriched in stalks (Fig. 4A). These numbers exceed previous estimates of the number of cell-type specific genes (Iranfar et al., 2001; Morrissey et al., 1984), so it was important to validate them by comparison to published findings. We found that nearly every gene previously described as prespore-specific (Iranfar et al., 2001) was found in our list of prespore enriched genes (supplemental data) and the ecmA gene showed a stalk-specific pattern (Fig. 4B). The correlation of these and other cell-type specific targets (supplemental data) with the published data validates our definition of cell-type-enriched genes. Our genome-wide survey of cell-type-enriched genes indicates that Dictyostelium cell-type differentiation is a more complex process than suggested by previous studies.

Expression of the spore coat genes is induced several hours before the expression of ecmA, suggesting that prespore cell differentiation precedes prestalk cell differentiation (Fosnaugh and Loomis, 1991; Morrissey et al., 1984; Williams et al., 1989). That observation initiated the idea that cell-type proportioning in Dictyostelium is governed by a mechanism of lateral inhibition in which the nascent prespore cells secrete a factor that imposes a prestalk cell fate on the remaining cells (Loomis, 1993). An expectation of this model is that the induction of many prespore specific genes should precede the expression of most prestalk-specific genes. To test that hypothesis, we selected the cell-type-enriched targets from Fig. 4A and traced them in the time-course data (Fig. 4B).

Most of the cell-type-specific genes described previously are not expressed at early stages of development and are induced after the transition to multicellular development. In fact, this pattern was the basis for defining the time of cell-type divergence to 10-12 hours after starvation (Morrissey et al., 1984). This predicted pattern of gene expression was prominent in the prespore- and prestalk-enriched RNA samples (Fig. 4B), but did not support the hypothesis that prespore gene expression precedes prestalk gene expression. However, we found an unexpected pattern that may be an antecedent to the prespore and prestalk cell differentiation. In the charts for PSP and PST (Fig. 4B), the expression of many prespore- and prestalk-enriched genes was higher-than-average early in development and lower-than-average later in development. The color scheme in Fig. 4B represents the level of gene expression relative to the mean expression level of that gene over the entire time course. Therefore, lower-than-average expression levels (blue) at the late time points do not necessarily mean that the gene is not expressed at those time points. With that in mind, notice that Fig. 4B describes groups of cell-type-enriched genes that are expressed at a higher-than-average level during growth and during the early stages of development, before the presumed time of cell-type divergence. The mRNA levels of these genes are reduced dramatically after the transition to multicellular development and they become enriched in one of the cell types.

The temporal analysis of the spore- and the stalk-enriched genes also revealed two patterns of expression. Genes with the expected pattern of lower-than-average expression during early development and higher-than-average expression during late development are clustered in the charts for SP and ST (Fig. 4B). These clusters consist of genes with an unexpected expression pattern: two peaks of higher-than-average expression, one at 2-6 hours of development and the other at the end of development. This pattern of gene expression may be an indication of a common physiological state that is shared by the pre-aggregation phase and the fruiting body phase.

The data presented here provide a new perspective on the complexity of the developmental program, the portion of the genome that is regulated during development and the timing and possible mechanisms of cell-type divergence in Dictyostelium. We found that almost 40% of the genes in the genome are regulated at some point in development in a strain- and nutrition-independent manner. The largest group of developmentally regulated genes, about 25% of the genome, exhibited a marked change in relative abundance of mRNA during the transition to multicellular development. Therefore, the differentiating cells in the multicellular organism appear to be very different physiologically from the vegetative unicellular amoebae.

Dictyostelium development is initiated by depletion of nutrients. In yeast, starvation induces a stress-related transcription pattern within 15 minutes and most of the additional responses occur within the first 4 hours (Jia et al., 2000; Natarajan et al., 2001). Our findings indicate that multicellular development in Dictyostelium is more than a starvation response. The finding that the largest transition in the transcription pattern occurs between 6-8 hours after starvation, coincident with the onset of the multicellular state, suggests that multicellular development is a distinct and specialized phase in the life cycle of Dictyostelium.

Comparing the expression array profiles of all the developmental time points also revealed two stages in which the transcriptional profile was relatively unchanged. The least amount of change was found during the first 6 hours of development (Fig. 2A). At that time, the cells sense starvation and begin to regulate the expression of genes that are necessary for cAMP signaling and chemotactic aggregation, but there is very little change in their gross morphology (Aubry and Firtel, 1999; Clarke and Gomer, 1995; Parent and Devreotes, 1999). The findings presented in Fig. 2A indicate that this period is not accompanied by vast changes in the transcriptional profile. This result is somewhat surprising in light of the complex signaling and motility mechanisms that are being established early in development (Aubry and Firtel, 1999; Parent and Devreotes, 1999). However, it is clear that the array detected transcriptional changes during the preaggregation stage. For example, many of the stalk- and spore-enriched genes shown in Fig. 4B are subject to regulation during the first 6 hours of development. Therefore, the data on the transcriptional transitions during the first 6 hours of development indicate a low level of change relative to the other periods in development. Another consideration is that many of the genes that encode regulatory elements may not be subject to dramatic changes and their effect on cell physiology may not be immediately reflected in the transcriptional profile.

The finger stage of development, between 14-18 hours after starvation, was the second least variable stage (Fig. 2A). At that stage the cells have already differentiated into the two major cell types and their proportions and spatial distribution are largely invariant (Kessin, 2001; Loomis, 1975; Williams et al., 1989). The low level of change in the transcriptional profile during that period is therefore consistent with the morphological findings.

Previous estimates of the number of developmentally regulated genes range between 300 and 1000, whereas our studies suggest that 2000-3000 genes are regulated during development. The previous estimates were based on mutagenesis studies, RNA hybridization studies and protein studies with two-dimensional gel electrophoresis (Alton and Lodish, 1977; Firtel, 1972; Loomis, 1978; Morrissey et al., 1984). The mutagenesis studies must have underestimated the number of developmentally regulated genes because they discovered only genes that are essential for proper developmental morphology (Loomis, 1978). The other studies relied on more subjective estimates of significance and therefore must have been fairly conservative. In addition, our definition of developmental regulation includes some genes that are expressed during all stages of growth and development as long as there is a sharp contrast between the early and the late expression levels, as well as genes that are downregulated during development. Such genes were not included in previous studies, but we included them in our definition because they are coordinately regulated with the developmentally induced genes and because they show consistent expression patterns regardless of the strain or the nutritional history of the developing cells (Figs 2B, 3; supplemental data).

The idea that transcriptional profiles are a reflection of cell physiology stems from the successful application of expression array data to the characterization of mutants in yeast (Hughes et al., 2000) and to the diagnosis of human tumor samples (Alizadeh et al., 2000; Bittner et al., 2000; Golub et al., 1999). Our data support that notion and extend it to the analysis of a developmental time-course. First, we found that the application of a dissimilarity calculation to the time-course expression array data conserved the correct temporal order of the samples despite the fact that temporal information was not included in the calculation (Fig. 2A). This finding indicates that the expression array data properly reflect the temporal progression of physiological change in the cells. Second, we found a conserved pattern of gene expression, which is insensitive to strain or nutritional history (Figs 2B, 3). This finding is consistent with the robust morphological and physiological aspects of the developmental program (Loomis, 1975). Third, we found that the degree of change in the transcriptional pattern was consistent with the degree of change in developmental morphology or physiology. We therefore conclude that transcriptional profiling is a useful tool in the analysis of development and we propose that it will be applicable to the analysis of mutations that affect development.

The data from the developmental time course can be analyzed in a variety of ways, one of which is to follow the pattern of expression of selected groups of genes. This approach was applied to the analysis of cell-type-enriched RNA (Fig. 4). It was somewhat surprising to find that many of the genes determined to be enriched in one cell type or another are expressed at high levels during growth or very early in development, much earlier than the expression of cell-type-specific structural genes and earlier than the generally accepted time of cell-type specification (Fosnaugh and Loomis, 1991; Haberstroh et al., 1991; Williams et al., 1987; Williams et al., 1989). We consider two opposing interpretations of this observation: (1) that the earlier expression of a gene has no bearing on its later function and no direct influence on the cell-type divergence event, or (2) that the earlier expression of the cell-type-enriched genes is determinative for cell-type specification. The first possibility indicates only that there may be a common physiological state between the pre-aggregation cells and the fully differentiated cells. For example, components of the cAMP signaling mechanism are utilized in different physiological contexts throughout development (Aubry and Firtel, 1999; Wang et al., 1999). However, if the early expression patterns are instructive then the early cohort of cell-type-enriched genes is either expressed in all cells and helps to dictate downstream events that determine cell differentiation, or is already expressed in a cell-type-specific manner. In either case, this finding suggests that cell-type divergence is affected by events that occur during growth and early in development.

The notion that growth-phase factors influence cell-type choice has a firm experimental basis. Cells grown on rich medium have a higher propensity to differentiate as spores than cells grown on poor medium (Blaschke et al., 1986; Thompson and Kay, 2000). Cells in the late G₂ or S-phase of the cell cycle have a higher propensity to differentiate as stalk cells and cells in the mid-G₂-phase have a higher propensity to differentiate as spores (Araki et al., 1997; Araki and Maeda, 1998; Gomer and Firtel, 1987; McDonald and Durston, 1984). These studies indicate that the vegetative cell physiology affects the subsequent cell-type choice. We propose that genes which are expressed at higher-than-average levels in the early time points and at cell-type-enriched, lower-than-average levels late in development may play a role in that developmental predisposition. For example, the rtoA gene is required for the coupling between cell-cycle phase and cell-type choice and is expressed at higher-than-average levels early in development and at lower-than-average levels late in development (Wood et al., 1996) and the tagA gene, whose expression peaks at 2-4 hours of development, is required for limiting the proportion of prestalk cells in slugs at 16 hours of development (J. Good, M. Cabral and A. Kuspa, unpublished observations). In light of these findings, our results suggest that differential gene expression during growth and early development presages or may help to direct cell-type divergence. One possibility is that the early expressed, cell-type-enriched genes described in Fig. 4B reflect mechanisms that inhibit cell-type specific differentiation. Such mechanisms may be present in all cells during the unicellular stage to prevent premature cell-type differentiation. After cell type divergence, the prespore specific mechanisms may inhibit prestalk cell differentiation and vice versa.

Our analysis of cell-type-enriched genes confirms and extends the published lists of cell-type specific genes (Iranfar et al., 2001; Kessin, 2001). The relatively large number of cell-type enriched genes is somewhat unexpected and it may be of interest to investigate their functional roles by more direct experimental approaches. Our analysis also suggests the presence of molecular mechanisms that function during growth and early development and have an effect on subsequent development and cell-type differentiation. Such mechanisms have been proposed before (Araki and Maeda, 1998; Blaschke et al., 1986; Gomer and Firtel, 1987; McDonald and Durston, 1984; Wood et al., 1996), but our analysis demonstrates that they involve a much larger portion of the genome than previous studies suggested.

By comparing the biological and the technical reproducibility of the experimental system we found that, in general, the biological variation was smaller than the technical variation and that the error rate in assigning a gene to a specific category was smaller than 5%. This variation should be considered when selecting individual genes as candidates for further studies, but it has little impact on the determination of the consensus group of 2021 developmental genes. The main reason for this is that our analysis is focused on the similarities between the four different strain and nutrition conditions. We therefore may have underestimated the number of developmentally regulated genes by 5-10%. Analysis of the transcriptional differences between the two strains and analysis of the developmental consequences of the nutritional history will require additional replication of all the experiments.

We consider the pattern of gene expression as a reflection of cellular physiology rather than as an indicator of the function of individual genes. This idea is supported by the finding that stress-induced gene expression and stress-related gene function were not correlated in the yeast Saccharomyces cerevisiae (Winzeler et al., 1999) and by the application of expression array data to the characterization of mutants (Hughes et al., 2000) and to the diagnosis of tumor samples (Alizadeh et al., 2000; Bittner et al., 2000; Golub et al., 1999). Based on this idea, our findings in Fig. 4B suggest the existence of overlapping physiological states between the 2-6 hour cells and the fully differentiated cells. These shared states may reflect the utilization of common molecular mechanisms to carry out different developmental functions.

In summary, our analysis has uncovered a large group of genes that exhibit an invariant pattern of expression despite the variable strain and nutritional histories of the cells. This finding indicates that Dictyostelium development is controlled by a robust physiological program and illustrates the power of applying transcriptional profiling to the analysis of development in multicellular organisms.

We thank Kirsten Kibler, Sijie Lu and Guokai Chen for unpublished cDNA and genomic DNA clones, Chunliu Su for technical assistance, Brian Desany for informatics assistance, the Human Genome Sequencing Center at Baylor College of Medicine for genomic DNA clones, the Department of Molecular and Human Genetics, Department of Biochemistry and the Office of Research at Baylor College of Medicine for financial support, and Stephen Elledge and Victoria Lundblad for useful comments on the manuscript. This work was supported by a Basil O’Connor research award from the March of Dimes Birth Defects Foundation (5-FY99-735; G.S.) and a grant from the National Institute of Child Health and Human Development, P01 HD39691-01 (A.K.). Work in Japan was supported by JSPS-RFTF (No. 96L00105 and No. 00L01412) and Grants-in-Aid for Scientific Research on Priority Areas A (No. 08283107) and C (No. 12206001).