HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Citation Map Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Reprints/Permissions Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Madras, N. Articles by Aubin, J. E. Search for Related Content PubMed PubMed Citation Articles by Madras, N. Articles by Aubin, J. E.

Modeling Stem Cell Development by Retrospective Analysis of Gene Expression Profiles in Single Progenitor-Derived Colonies

^a Department of Mathematics and Statistics, York University, Toronto, Ontario, Canada;
^b Department of Statistics and Actuarial Science, University of Waterloo, Waterloo, Ontario, Canada;
^c Institute of Biomaterials and Biomedical Engineering, and
^d Department of Anatomy and Cell Biology, University of Toronto, Toronto, Ontario, Canada

Key Words. Stochastic model • Lineage analysis • Osteoprogenitor • Global amplification PCR • Stem cells

Jane E. Aubin, Ph.D., Department of Anatomy and Cell Biology, Faculty of Medicine, University of Toronto, Room 6255 Medical Sciences Building, One King's College Circle, Toronto, Ontario M5S 1A8, Canada. Telephone: 416-978-4220; Fax: 416-978-3954; e-mail: jane.aubin{at}utoronto.ca

	ABSTRACT

Top Abstract Introduction Materials and Methods Results Discussion Conclusions Appendix References

 
The process of development of various cell types is often based on a linear or deterministic paradigm. This is true, for example, for osteoblast development, a process that occurs through the differentiation of a subset of primitive fibroblast progenitors called colony-forming unit-osteoblasts (CFU-Os). CFU-O differentiation has been subdivided into three stages: proliferation, extracellular matrix development and maturation, and mineralization, with characteristic changes in gene expression at each stage. Few analyses have asked whether CFU-O differentiation, or indeed stem cell differentiation in general, may follow more complex and nondeterministic paths, a possibility that may underlie the substantial number of discrepancies in published reports of progenitor cell developmental sequences. We analyzed 99 single colonies of osteoblast stem/primitive progenitor cells cultured under identical conditions. The colonies were analyzed by global amplification poly(A) polymerase chain reaction to determine which of nine genes had been expressed. We used the expression profiles to develop a statistically rigorous map of the cell fate decisions that occur during osteoprogenitor differentiation and show that different developmental routes can be taken to achieve the same end point phenotype. These routes appear to involve both developmental "dead ends" (leading to the expression of genes not correlated with osteoblast-associated genes or the mature osteoblast phenotype) and developmental flexibility (the existence of multiple gene expression routes to the same developmental end point). Our results provide new insight into the biology of primitive progenitor cell differentiation and introduce a powerful new quantitative method for stem cell lineage analysis that should be applicable to a wide variety of stem cell systems.

	INTRODUCTION

Top Abstract Introduction Materials and Methods Results Discussion Conclusions Appendix References

 
Most current models of mesenchymal stem cell and osteoprogenitor differentiation assume a hierarchical arrangement, with an expanding series of cells leading from a small number of stem cells through a larger but finite number of committed progenitors and a still larger number of mature cells with a finite lifespan [1–4]. Although recently the role of plasticity or transdifferentiation events of cells usually considered terminally differentiated, has been the subject of intense interest [5, 6], there is an apparently restricted osteoprogenitor detectable in colony assays performed on cultures of bone or bone marrow cell isolates under particular standardized growth conditions. The ability to identify and isolate osteoprogenitors as pure populations at discrete and identifiable stages has been difficult and lags behind what is possible in certain other lineages (e.g., hematopoietic). Thus, a retrospective assay based on differentiated cell function/phenotype—a colony assay—combined with other biochemical and molecular assays has become widely used and has contributed substantially to the understanding of osteoblast development [5]. Briefly, bone nodule or bone colony formation represents the end-stage of a proliferation-differentiation sequence of a low frequency of cells within unfractionated primary populations of bone-(e.g., rat calvaria [RC]) or bone marrow-derived cell cultures [7–12].

Morphologically and molecularly recognizable osteoblasts and bone nodules form in such cultures at predictable and reproducible times after plating [5, 13]. The process is often subdivided into three stages: A) proliferation; B) extracellular matrix development and maturation, and C) mineralization, with characteristic changes in gene expression at each stage [14]. Expression of osteoblast-associated genes (e.g., type I collagen [COLL I], alkaline phosphatase [ALP], osteopontin [OPN], bone sialoprotein [BSP], osteocalcin [OCN], parathyroid hormone/parathyroid hormone-related protein [PTH/PTHrP] receptor [PTH1R] and others) is asynchronously acquired and/or lost as the progenitor cells differentiate and matrix matures and mineralizes. In general, ALP is thought to increase then decrease when mineralization is well progressed, OPN is high during both proliferation and differentiation stages and is upregulated prior to certain other matrix proteins, including BSP and OCN, with OCN first appearing in postproliferative osteoblasts approximately concomitant with mineralization [9, 15] (Fig. 1A). However, many of these conclusions are based, in large part, on models in which differentiation and bone nodule formation are not synchronous throughout the population, and many of the available osteoblast markers used are not restricted to osteoblasts but may be expressed at substantial levels by other cells (e.g., fibroblasts) in the population, leading to discrepancies in the reported osteoblast developmental sequence in the published literature [5, 16]. In addition, although limiting dilution analysis of bone-derived cell populations, such as those used here (i.e., the RC cell model), has suggested that osteoprogenitor cell differentiation/bone colony formation is a single-hit phenomenon with only one cell type (likely the osteoprogenitor cell) limiting [7], we have also reported evidence that osteoblast differentiation in this model appears to have a stochastic component [9, 15, 17].

View larger version (25K): [in this window] [in a new window]   Figure 1. Models of the osteoblast developmental sequence. A) The established model of osteoblast development follows a linear deterministic sequence in which expression of osteoblast-associated genes (numbered 1-9; see text for definitions) is acquired sequentially as differentiation proceeds. B) A new model based on our analysis. Each shaded circle represents an event in the development of a typical osteoblast colony. A circle containing an unbracketed number denotes the activation of the gene with that label (Table 1). An arrow from event i to event j indicates that event i always precedes event j. The arrow is dashed if there is some evidence for precedence, but it is statistically inconclusive; shorter dashes correspond to less statistical evidence. Bracketed numbers denote genes that must have been activated previously. A circle's horizontal (left to right) position roughly corresponds to the average time at which its event occurs. The circle marked with x represents an early event, perhaps related to expression of gene 1. Genes 3, 4, and 6 are activated independently (and "exchangeably") after this early event. Genes 2 and 5 are activated independently of all other genes analyzed.

	MATERIALS AND METHODS

Top Abstract Introduction Materials and Methods Results Discussion Conclusions Appendix References

 
Cell Culture, Osteoblastic Colony Isolation, and cDNA Preparation
Cells were enzymatically isolated from the calvariae of 21-day Wistar rat fetuses, by sequential digestion with collagenase, and plated at low density as described [8]. Osteoblast colonies for analysis were chosen either based on their morphology [8, 20] or from a replica plating technique [15, Liu and Aubin, submitted]. In either case, dishes were rinsed with phosphate-buffered saline, a cloning ring was placed around selected colonies, and the cells were released with 0.01% trypsin (when matrix lacked mineral) or a 1:1 mixture of 0.01% trypsin and collagenase (when osteoid was mineralizing); the enzyme(s) were neutralized by adding minimal essential medium (-MEM) containing 15% fetal bovine serum. Total RNA was extracted using a miniguanidine thiocyanate method, and cDNA was synthesized by oligo(dT) priming, poly(A)-tailed, and amplified by PCR with oligo(dT) primer as described [8, 21]. The colonies analyzed were collected from 12 independent cell isolations (each resulting from a pool of at least 25 calvariae) and replica plating experiments.

Southern Blots and Hybridization
Amplified cDNA (5 µl) was run on 1.5% agarose gels, transferred onto 0.2 µm-pore-size nylon membrane (ICN; Costa Mesa, CA), and immobilized by baking at 80°C for 2 hours. Prehybridization, hybridization with labeled probes for mRNAs of interest (Table 1), and probe signal quantification were performed as described [8]. The signal intensity for each probe was standardized against total cDNA [8, 20]. For the analysis reported here, each colony is reported as expressing (+) or not expressing (–) a particular gene of interest relative to background; the actual level of expression is not considered (Fig. 2).

View larger version (11K): [in this window] [in a new window]   Figure 2. A schematic of the Southern lineage blot describing expression of nine genes in 99 colonies on a + (expressed; black boxes; value 1) or – (not expressed; white boxes; value 0) scale.

	RESULTS

Top Abstract Introduction Materials and Methods Results Discussion Conclusions Appendix References

 
Development of the Analytical Framework
Expression levels of the nine genes listed and numerically labeled in Table 1 were measured from 99 osteoprogenitor/osteoblast colonies. Our analysis focused on genes expressed above a background level (Fig. 2); we assumed that it was harder to measure magnitudes of expression levels accurately than to decide which genes were expressed at all. This information is encoded in an array of binary variables B_k⁽ⁱ⁾ (i = 1,...,99 , k = 1,...,9), defined to be 1 (respectively, 0) if expression of gene k is (respectively, is not) observed in the i-th colony. Time (days) is indexed so that the experiment begins at time 0 and ends at time . We define the variable T_k⁽ⁱ⁾ to be the time at which gene k is first expressed in the i-th colony (i = 1,...,99 , k = 1,...,9). These times, T_k⁽ⁱ⁾, are not observed directly; we can only observe whether or not they occur before the end of the experiment. That is, if there were no errors in detecting gene expression, then for each i and k we could say that T_k⁽ⁱ⁾ if and only if B_k⁽ⁱ⁾ = 1. To allow for detection errors, we define q_k to be the "false negative" probability that we do not detect the expression of gene k in a particular colony when in fact that gene is expressed there. In terms of conditional probabilities, this says that q_k = Pr{B_k⁽ⁱ⁾ = 0 | T_k⁽ⁱ⁾ } = 1 - Pr{B_k⁽ⁱ⁾ = 1 | T_k⁽ⁱ⁾ }. We assume that there are no "false positive" errors, i.e., Pr{B_k⁽ⁱ⁾ = 0 | T_k⁽ⁱ⁾ > } = 1.

The overriding tenet of our model is that the 99 colonies are stochastically independent and governed by the same probability law. Formally, for each i = 1,...,99, we define the vector V⁽ⁱ⁾ = (B₁⁽ⁱ⁾,...,B₉⁽ⁱ⁾,T₁⁽ⁱ⁾,...,T₉⁽ⁱ⁾) consisting of the i-th colony's 18 associated variables. Then V⁽¹⁾,...,V⁽⁹⁹⁾ is a sequence of 99 independent random vectors drawn from some 18-dimensional multivariate probability distribution. Our goal is to describe the salient features of this multivariate distribution. When discussing the distribution (as opposed to the 99 observations drawn from it) we omit the superscripts and use the random vector V = (B₁,...,B₉,T₁,...,T₉).

A key aspect of our data is displayed in the 9 x 9 matrix E shown in Figure 3. For each pair of genes k and l, let E_kl be the number of colonies in which we detect the expression of gene k but do not detect the expression of the gene l. Formally, E_kl is the number of i's such that B_k⁽ⁱ⁾ = 1 and B_l⁽ⁱ⁾ = 0. Observe that E_kl is 0 and E_1k is large for every k = 2,3,...,9. This strongly suggests that gene 1 (OPN) is expressed before every other gene, i.e., Pr{T₁ < T_k} = 1 for every k = 2,...,9. In contrast, E₃₂ = 19 and E₂₃ = 20, so it seems unlikely that one of genes 2 or 3 must always precede the other. This paper shows how to analyze such claims statistically. For the most part, our methods are novel adaptations of standard techniques in statistics and probability, which may be found in [22] and many other general statistics references.

View larger version (28K): [in this window] [in a new window]   Figure 3. The matrix (Ekl). The 9 x 9 matrix E summarizes observed precedences: each entry, Ekl is the number of colonies in which expression is observed for gene k but not for gene l (genes are numbered as in Table 1).

It is reasonable to assume that q_k increases as the average expression level decreases, and that genes with similar expression levels should have similar q_k's. Therefore, from Table 1, we obtain

(Equation 1)

In Section A.1 of the Appendix, we generalize the method used above for bounding q₁ to obtain the bounds

(Equation 2)

each with (at least) 95% confidence. Importantly, these error bounds were obtained by purely statistical means. We emphasize that they are (conservative) bounds rather than estimates; the true q_k values could be significantly smaller.

Can Bone Development Be Explained by a Linear/ Deterministic Expression of Genes?
Although some aspects of the order of expression events are clear from the data (e.g., OPN is expressed early and OCN late), many others are not (e.g., does platelet-derived growth factor receptor alpha [PDGFR] precede fibroblast growth factor receptor 1 [FGFR1] or vice versa?). We first asked whether the data are consistent with a deterministic linear order of the expression events, with the apparent inconsistency being due to observational errors.

To answer this, we suppose that the events T₁,...,T₉ always occur in the same order. For each k = 1,...,9, write (k) for the label of the gene (as listed in Table 1) that is always k th in this order. For example, if COLL always occurs second, then (2) would be 6. Then, we would have

(Equation 3)

The sequence, (1),...(9), is a permutation of 1,2,...,9. We ask if there is any permutation such that the data are statistically consistent with Equation 3.

Each permutation would imply certain false negatives in the data. For example, suppose that Equation 3 holds, and that in colony 14 we observe B₍₃₎⁽¹⁴⁾ = 0 and B₍₇₎⁽¹⁴⁾= 1. Since gene (3) is always expressed before gene (7), and since there are no false positives, the observation B₍₃₎⁽¹⁴⁾ = 0 must be a false negative error. For each permutation, we computed the minimum number of false negatives that would be necessary in order to explain the data if Equation 3 were true (this is the number of values of i and k such that B⁽ⁱ⁾_(k) = 0 and there exists a value l > k such that B⁽ⁱ⁾_(l) = 1). The fewest such errors is 85, attained by two different permutations: 1,2,5,4,3,6,8,7,9 (i.e., OPN, PDGFR, PTHrP, PTH-R, FGFR1, COLL, BSP, ALP, OCN) and 1,5,2,4,3,6,8,7,9. In Appendix A.2, we show that this many errors is very unlikely, which in turn provides strong evidence against Equation 3 and leads us to reject the possibility that there is a linear sequence of genes that accounts for the development of a bone nodule colony.

Is the Expression of One Gene Independent of the Expression of Other Genes?
We next determined, for each pair of genes k and m, whether their expressions B_k and B_m were independent. (The error probabilities q_k and q_m do not play a role here, since all errors are assumed to be independent.) We excluded gene 1 from this analysis, because B₁⁽ⁱ⁾ = 1 for every colony i, which prevents any meaningful testing of independence. To illustrate the method, consider the data for PDGFR and FGFR1 (k = 2 and m = 3) in Table 2. Testing it for independence using Fisher's exact test (two-tailed) [23] gives a p value of 0.788. Table 3 lists the results for all pairs. Any p value below 0.05 is sufficient to reject the independence hypothesis on any single test. However, since we are testing 28 different pairs, we can reasonably expect to see one or two p values below 0.05 purely by chance, even when the corresponding pair of variables is independent. Thus, we need to be more cautious unless the p value is very small indeed. See Section A.3 of the Appendix for further discussion of this matter. In particular, that section shows that any p value below 0.0018 can be understood as a clear case for rejecting independence, while values between 0.0018 and 0.05 should be regarded as inconclusive indications of dependence. However, of the five "inconclusive" values that appear in Table 3, we expect all but at most one or two to be truly caused by dependence.

View this table: [in this window] [in a new window]   Table 2. The number of events where both FGFR1 and PDGFR were detected (1,1), either FGFR1 or PDGFR were detected (0,1 and 1,0), and neither were detected (0,0).

The following model illustrates how positive correlations could arise. Consider three events (C, D, and E) which occur at the random times T_C, T_D, and T_E. Event C appears first. Once C occurs, then D and E each occur after an additional random (and independent) amount of time. Letting X_CD be the duration of time between events C and D, and X_CE the time between C and E, then T_D = T_C + X_CD and T_E = T_C + X_CE. Assume that X_CD, X_CE, and T_C are independent. Then T_C, T_D, and T_E are mutually positively correlated (in fact, Cov(T_i, T_j) = Var(T_C) > 0 when i j). Moreover, consider the binary variables B_C, B_D,and B_E (where B_k is 1 if T_k and is 0 otherwise). Then each of Cov(B_C, B_D), Cov(B_C, B_E), and Cov(B_D, B_E) is nonnegative, by the FKG Inequality [24].

Our analysis rules out situations in which negative correlations could arise. One such scenario has two alternative developmental pathways, with gene D occurring only on one pathway and gene E occurring only on the other. Then no colony would have B_D = 1 = B_E, and hence Cov(B_D, B_E) < 0 (since E(B_DB_E) = 0).

Can Genes Replace Each Other in the Same Developmental Program (Exchangeability)?
We showed above that genes 3, 4, and 6 are mutually dependent, in fact, positively correlated. We can determine if there is a preferred order in which these genes need to be expressed to result in bone nodule development. In Figure 3, observe that the values of E_ij for i, j {3,4,6} are all approximately equal (all between 11 and 15). This suggests no preferred order among these three genes. We formalize this as follows.

Hypothesis A. The genes 3, 4, and 6 occur in random order. More precisely, as soon as some (early) event occurs (say at a random time T_*), each of genes 3, 4, and 6 waits for a random amount of time before it is expressed. Moreover, these three waiting times are i.i.d. (independent and identically distributed).

That is, there exist i.i.d. random variables X₃, X₄, and X₆, all independent of T_*, such that T_k = T_* + X_k for k = 3,4,6. The event at time T_* could be the expression of gene 1 (i.e., T_* = T₁), the expression of some other gene not examined in this experiment, or some other event.

Under Hypothesis A, T₃, T₄, and T₆ are not independent, but they are identically distributed. In fact, T₃, T₄, and T₆ are exchangeable random variables under Hypothesis A. This means that the probability law of the vector (T₃,T₄,T₆) is the same as that of (T₄,T₆,T₃) and of its four other permutations. For example, exchangeability implies that Pr{B₃ = 1, B₄ = 0, B₆ = 1} = Pr{B₄ = 1, B₆ = 0, B₃ = 1} and similarly for every other permutation.

In Section A.4 of the Appendix, we use this property to devise a statistical test of Hypothesis A. As described there, the tests show that the data are clearly consistent with the hypothesis of exchangeability implied by Hypothesis A. Due to the nature of statistical hypothesis testing, we cannot say that Hypothesis A has been accepted, but at most we can say that our data does not refute it. (Anything more would require information about the power of the test, which is beyond the scope of this paper; indeed, there is no obvious choice of alternative distributions for which one should do power computations.) That is, during osteoblast development, we hypothesize that there is a set of genes necessary for phenotypic maturation, but within this set there is flexibility in the order in which the genes are expressed.

In Section A.4, we also test exchangeability of the random variables {B₂, B₃, B₄, B₅, B₆ } corresponding to the five "early" genes. That is, we tested whether we could extend Hypothesis A to the group of five genes 2,3,4,5,6. As explained in Section A.4, this was rejected.

Other Precedence Constraints
In this section, we consider the three "late" genes: ALP, BSP, and OCN. First, for gene 7 (ALP), observe from Figure 3 that E_7k is much smaller than E_k7 for each of k = 2,3,4,5,6. Therefore, ALP typically is expressed after genes 2 through 6. In particular, E₇₆ = 0, which suggests

Hypothesis B. Gene 6 (COLL) is always expressed before gene 7 (ALP). Moreover, none of genes 2, 3, 4, or 5 are a prerequisite for the expression of gene 7.

By our tests of independence, genes 2 and 5 are clearly not prerequisites for gene 7. Recall that we accepted Hypothesis A, that B₃, B₄, and B₆ are exchangeable among themselves. If gene 6 has a special relationship with gene 7 that genes 3 and 4 do not (in particular, if Hypothesis B is true), then we would expect the following to be false:

Hypothesis B*. The variables B₃, B₄, and B₆ are exchangeable with respect to the joint distribution of genes 3, 4, 6, and 7.

As described in Section A.5 of the Appendix, our tests indicate that Hypothesis B* should be rejected, but we do not reject the exchangeability of B₃ and B₆ with respect to the joint distribution of genes 3, 4, and 7. This is consistent with Hypothesis B. Further data are needed before one can conclusively accept Hypothesis B, however.

Finally, we considered genes 8 (BSP) and 9 (OCN). In particular, we looked at the possibilities that Pr{T₈ < T₇} = 1 and that Pr{T_k < T₉} = 1 for every k 8, as suggested by the matrix E of Figure 3. We found that the statistical evidence was not strong enough to either confirm or reject these possibilities. See Section A.5 of the Appendix for details.

Summary of Expression Analysis
Combining the results of the statistical analyses, we have constructed a new model of osteoblast differentiation that is consistent with the observed relationships of expression among particular genes or subsets of genes (Fig. 1B). Our analysis indicates that the expression of genes 2 (PDGFR) and 5 (PTHrP) may occur independently of other aspects of osteoblast development. Conversely, gene 1 (OPN), which precedes the expression of all of the other bone-related genes, may initiate the expression of genes 3 (FGFR1), 4 (PTH-R), or 6 (COLL). Once all the genes in this group are expressed, further maturation depends upon the expression of gene 7 (ALP), which may be contingent upon gene 8 (BSP) being independently activated. Finally, when genes 6, 7, and 8 are expressed, gene 9 (OCN) is expressed and the gene cascade associated with the clonal development of osteoblast phenotype is accomplished.

	DISCUSSION

Top Abstract Introduction Materials and Methods Results Discussion Conclusions Appendix References

 
In this paper, we readdressed the established linear/deterministic view of osteoblast development (Fig. 1A) by a retrospective analysis of gene expression profiles in individual osteoprogenitor clones that ultimately gave rise to mature osteoblasts making mineralized osteoid, i.e., the bone nodule colony-forming assay. The analysis provides an analytical framework of osteoprogenitor differentiation and is unique because it is based solely on establishing expression relationships among genes at different stages of development and not on any preconceived notion of their expression order. A primary advantage of this approach is that it allows identification of developmental gene expression patterns for genes with both known and unknown expression profiles. It also allows elucidation of the developmental history of a colony without an exhaustive examination of each cell fate division that takes place during its development [24]. The principles established, while tested for osteoblasts, should be widely applicable to other lineages.

The calculated gene expression profiles of individual clones were compared with those based on established osteoprogenitor developmental paradigms (Figs. 1A and 1B). We first asked whether there was a predetermined order to the expression of a set of nine genes associated with osteoblast development. The analysis, which depended on a reduction of sequential expression errors to levels below those due to observational error, did not support a linear or deterministic sequence. Instead, although all the tested colonies eventually developed into fully mineralized bone nodules, the range of expression profiles at different stages of development suggested that it may be possible for osteoprogenitors to differentiate into mature osteoblasts using different developmental routes (Fig. 1B). To confirm and extend this possibility, we examined the correlation between the expression of a particular gene and any other individual gene. In some cases, strong positive correlations were found between particular pairs of genes, e.g., expression of gene 7 (ALP) required expression of gene 6 (COLL). In other cases, such as genes 2 (PDGFR) and 5 (PTHrP), expression appeared independent of any other gene tested. One hypothesis is that acquisition of expression of these latter genes occurs at very early times in osteoblast development, i.e., at a stage prior to the earliest replica plating data points we have available. Alternatively, the expression profiles may be indicative of, or a prefacing of, potential development along pathways other than osteoblastic (e.g., chondrocytic or adipocytic) as may occur from multipotent progenitors (i.e., bipotential cells or mesenchymal stem cells) present in the population [1, 5], but these may not progress under the culture conditions in which osteoblast development predominates (and was thus selected for analysis).

Our revised hierarchy of clonogenic osteoprogenitor differentiation (Fig. 1B) is consistent with two hypotheses of developmental history underlying osteoblast formation. One is that there is more than one pattern of sequential gene expression profiles that results in the same type of mature cell (developmental flexibility hypothesis). The other is that, as each colony matures from an undifferentiated progenitor cell, some divisions result in pathways that are either not supported by the environment (leading to lineage extinction) or give rise to other cell types that may exist in the final colony (i.e., cells expressing gene 2 (PDGFR) or gene 5 (PTHrP)) but do not notably contribute to its overall phenotype under conditions where osteoblast development predominates (also, lineage extinction). At present, our data do not allow us to select between these two possibilities and it is possible that both occur. There is, for example, some evidence for mixed lineages in a small proportion of osteoblast colonies (e.g., adipocyte-osteoblast colonies, [25, 26]). It is also interesting to note that colony size appears to have no (detectable) impact on expression profiles obtained, i.e., as raised above, there is a stochastic component underlying CFU-O differentiation, such that a small bone colony can be as "mature" as a big colony [9, 15]. In addition, the variability seen in the rate of colony development is such that the gene expression "snapshot" taken during replica plating should give us an adequate description of all potential cell fate decisions during colony development (with the possible exception of very early genes such as gene 1, OPN). Taken together, our results argue against the existence of "rogue" differentiation routes, i.e., the expression of genes 2 and 5, or expression of genes 3 and 4, as deadends in the differentiation process, and for the possibility that their expression is required earlier during osteoblast development than was measured using our replica plating strategy. However, our results also suggest that some gene combinations are exchangeable and that higher order interactions among genes exist during osteoblast development. They support the hypothesis that great complexity underlies the interactions among groups of genes, and that further analysis will require larger gene (e.g., microarray technologies) and cellular (e.g., substantially more progenitor colonies) data sets. Interestingly, in contrast to the gene expression programs that appear to underlie early events in osteoblast development, global late gene expression profiles appear to be more fixed (Fig. 1B). This is consistent with other well-defined gene expression sequences during the maturation of committed progenitor cells (including, e.g., erythrocytes) and may be usual for cells that have committed to a particular mature state and lost their developmental flexibility. They also support the view that different developmental routes may contribute to and/or account for the substantial molecular heterogeneity seen amongst mature osteoblasts in vitro [20] and in vivo [27] (as predicted in [12]).

Although a stochastic component defines their development, osteoblast differentiation appears to be a default pathway for progenitor cells in RC populations cultured under the conditions we describe here [5]. Whether survival of particular cell types results from a microenvironment effect, or commitment to particular lineages is initiated by the exposure of an uncommitted cell to a particular signal that, upon exit of the cell from the stem cell state, stabilizes the subsequent identity of the progeny, is not yet known. This so-called directive mechanism clearly depends on the interactions between the cellular sensing systems—typically cell surface receptors—and the cellular microenvironment. Substantial evidence exists that both these mechanisms can occur and that the particular mechanism used may be dependent on the tissue system/cellular microenvironment. For example, in hematopoiesis, several studies suggest that exposure to growth factors may not be obligatory for the differentiation of primitive cells and that, at least under certain conditions, the identity of the differentiated cell population may be intrinsically determined [28]. Particularly interesting, in this regard, is the recent demonstration that coexpression of multiple lineage-restricted genes precedes commitment in multipotent progenitors [29–31]. This multilineage "priming" process [30] is consistent with the flexibility in the gene expression profiles we have seen during osteoprogenitor development and implies that the commitment of a multipotent cell to a particular pathway may reflect the stabilization of a particular subset of a group of expressed molecules. The stabilization process may occur in a stochastic manner in the absence of a particular instructive signal. Conversely, upon commitment of an undifferentiated cell, an instructive signal may stabilize a particular set or subset of expressed transcription factors, resulting in the production of specific cell types.

	CONCLUSIONS

Top Abstract Introduction Materials and Methods Results Discussion Conclusions Appendix References

 
We have presented a novel and powerful way to analyze gene expression relationships during clonal stem or progenitor cell differentiation. This approach has not only given us insight into the cascade of genes expressed during osteoprogenitor development, but has, for the first time, provided quantitative support to observed developmental flexibility during stem and primitive progenitor cell differentiation. Our approach addresses the general need for an analytical/quantitative way to examine stem cell differentiation programs, should be useful to address questions of stem or progenitor cell commitment in response to changes in microenvironmental stimuli, and is a powerful way of analyzing the large amount of information created in genetic databases of stem cell developmental programs [32].

	APPENDIX

Top Abstract Introduction Materials and Methods Results Discussion Conclusions Appendix References

 
A.1 Derivation of Bounds on Experimental Error
To get bounds on q₂,...,q₉, we look at the 0 entries in matrix E (Fig. 3). First, consider E₉₄, which is 0. It turns out that there were 30 colonies in which both genes 4 and 9 were detected. Let A₉₄ be the number of colonies in which genes 4 and 9 were both expressed. Since there are no false positives, we know that A₉₄ 30. The probability of detecting 9 and missing 4 in one of these A₉₄ colonies is q₄ (1 – q₉). But there were no misclassifications of this kind (since E₉₄ = 0), and the probability of this (given A₉₄) is (1 – q₄(1 – q₉))^A94. As above, we are 95% confident that

, i.e., that

(Equation A1)

We believe q₉ q₄, so we deduce that q₉ (1 – q₉) 0.095, and thus q₉ 0.106 with 95% confidence. Inserting this bound into Equation A1 gives q₄ 0.106.

We repeat these calculations for E₉₆, E₉₇, and E₉₈, which are all 0. Each of these cases has 30 colonies in which both genes are detected, implying that q_k 0.106 with high confidence for k = 6,7,8. Alternatively, we could apply q₄ 0.106 to Equation 1 and deduce q_k 0.106 for every k, with high confidence. An analogous calculation for E₇₆ yields q₆ 0.064. By Equation 1, the same bound should hold for q₈. In summary, we have shown

(Equation A2)

each with (at least) 95% confidence.

A.2 Linear/Deterministic Expression of Genes
As described in part 3 of the Results section, we ask if there can be a permutation such that the data are statistically consistent with Equation 3, i.e.,

Recall that each permutation implies certain false negatives in the data. The fewest such errors is 85, attained by two different permutations: 1,2,5,4,3,6,8,7,9 and 1,5,2,4,3,6,8,7,9. We shall show that this many errors is very unlikely for any permutation, which in turn provides strong evidence against Equation 3. We split the 99 x 9 observations B_k⁽ⁱ⁾ into three groups: Group A contains B₁⁽ⁱ⁾, i = 1,...,99 (i.e., all OPN data); Group B contains B₆⁽ⁱ⁾, B₈⁽ⁱ⁾, i = 1,...,99 (all COLL and BSP data); and Group C contains B_k⁽ⁱ⁾, k = 2,3,4,5,7,9, i = 1,...,99 (all other data). Recall that our worst-case bounds on the error probabilities for groups A, B, and C, respectively, are q_A = 0.030, q_B = 0.064, and q_C = 0.106. Let M_A, M_B, and M_C denote the true number of expression events {T_k⁽ⁱ⁾ } that occurred (whether observed or not) in Groups A, B, and C, respectively. We observed 99 of these events in Group A, 139 in Group B, and 386 in Group C. Hence, M_A = 99, M_B 139, and M_C 386. Now let F_A, F_B, and F_C be the number of false negatives in the three groups, i.e., F_A = M_A – 99, F_B = M_B – 139, F_C = M_C – 386. If Equation 3 holds, then 85 F_A + F_B + F_C = M_A + M_B + M_C – 624. For each y = A,B,C, let G_y, be a binomially distributed random variable with parameters n = M_y and p = q_y. Then G_y is stochastically greater than F_y. It turns out that Pr{G_A + G_B + G_C M_A + M_B + M_C – 624} is about 10^–3 or smaller for all values of M_A, M_B, and M_C such that M_A = 99, M_B 139, and M_C 386, and M_A + M_B + M_C 624 + 85. (This can be calculated using the normal approximation for (G_A + G_B + G_C – µ) / , where µ = q_AM_A + q_BM_B + q_CM_C and ² = q_A (1 - q_A)M_A + q_B (1 – q_B)M_B + q_C (1 – q_C) M_C). That is, the probability of seeing enough errors is negligible. Hence we reject the possibility that there is a linear sequence of genes that accounts for the development of a bone nodule colony.

A.3 Testing Independence
In part 4 of the Results section, we performed Fisher's exact test on each pair B_k and B_m (for 2 k < m 9), which is 28 tests. For a single pair, one would normally reject independence when the p value is less than 0.05. However, in 28 tests, there is a reasonable chance of seeing one or more p values below 0.05 by chance alone, even when independence holds. To be cautious, we can choose a threshold, z, such that Pr{min{p₁, K, p₂₈} < z} is at most 0.05, where p_i is the p value from the i-th test under the assumption that the null hypothesis is true. Since Pr{p_i < z} = z (up to discretization error), we can choose z = 0.05/28 = 0.0018 (by Bonferroni's inequality). Therefore, any p value between 0.05 and 0.0018 should be regarded with some caution as being inconclusive. Any p value below 0.0018 may be safely used to reject independence. There are nine such values in Table 3. There are five values between 0.0018 and 0.05 in Table 3 that must be regarded with caution. However, it seems highly unlikely that all five of these values are due to chance resulting from five null hypotheses being true. Indeed, there are at most 28 – 9 = 19 cases where the null hypothesis is not conclusively rejected. If we assume that the null hypothesis holds in all 19 of these cases, and if we further assume that the p values are approximately independent, then we can say that the number of p values below 0.05 in these 19 cases has an approximate Poisson distribution with mean 19 x 0.05 = 0.95. This implies that the approximate probability of seeing j or more p values below 0.05 is 0.246 for j = 1, 0.071 for j = 2, and 0.016 for j = 3. We conclude that all but at most two of the inconclusive p values are likely due to the falsity of the independence null hypothesis.

A.4 Exchangeability
We now develop the statistical test of Hypothesis A based on exchangeability, as described in part 5 of the Results section. First, for each colony i = 1,K, 99, let S₃₄₆⁽ⁱ⁾= {r {3,4,6}: B_r⁽ⁱ⁾ = 1}, which is the random subset of {3,4,6} corresponding to the genes that have been detected in the i-th colony by time . Next, let N₃₄₆⁽ⁱ⁾ = B₃⁽ⁱ⁾ + B₄⁽ⁱ⁾ + B₆⁽ⁱ⁾, which is the size of S₃₄₆⁽ⁱ⁾. The number of colonies in which N₃₄₆⁽ⁱ⁾ was 0 (respectively, 1, 2, 3) was 8 (respectively, 10, 28, 53). For each subset C of {3,4,6}, let n[C] be the number of colonies i such that S₃₄₆⁽ⁱ⁾ = C.

For each k, exchangeability implies that, conditioned on the event that N₃₄₆ = k, the random set S₃₄₆ is uniformly distributed among all k-element subsets of {3,4,6}. We test this uniformity using the classical goodness-of-fit statistic GF = ³_r=1 (n[C_r] – E(n[C_r]))² / E(n[C_r]), where the C_rs are the k-element subsets of {3,4,6} and E(·) denotes expected value. For k = 2, the observed values (n[{3,4}] = 8, n[{3,6}] = 10, and n[{4,6}] = 10) give E(n[C_r]) = 28/3 and GF = 0.29. Under Hypothesis A, the approximate distribution of GF is x²₂ (chi-squared with 2 degrees of freedom). A chi-square table shows that 0.29 has a p value near 0.10. For k = 1, the data n[{3}] = 3, n[{4}] = 5, and n[{6}] = 2 give E(n[C_r]) = 10/3 and GF = 1.40. In this case, there are too few observations for the chi-square approximation to be valid, so we used a simple Monte Carlo simulation to compute Pr{GF 1.40} = 0.63 under Hypothesis A. That is, we simulated a million realizations of the test statistic GF under Hypothesis A and found that 63% of simulated values were greater than 1.40.

The k = 1 and k = 2 tests show that the data are clearly consistent with the hypothesis of exchangeability implied by Hypothesis A. We also used this method to test exchangeability of the random variables {B₂, B₃, B₄, B₅, B₆}, corresponding to the five "early" genes. We rejected uniformity of S₂₃₄₅₆ for k = 3 (with p value 0.007), marginally accepted it for k = 1, and accepted it for k = 2 and k = 4. Overall, then, the data are not consistent with exchangeability of all five of {B₂, B₃, B₄, B₅, B₆}.

A.5 Precedence Constraints for Later Genes
Here we describe the methods referred to in part 6 of the Results section. Hypothesis B* from that section implies that (conditioned on the event {N₃₄₆ = 2}) the expression of gene 7 is independent of the set S₃₄₆. We therefore test independence of B₇ and S₃₄₆ conditioned on {N₃₄₆ = 2}, using the data in Table 4. Fisher's exact test gives the p value 0.045, which marginally rejects Hypothesis B*. The analogous test for independence of B₇ and S₃₄₆ conditioned on {N₃₄₆ = 1} gives little information since the row sum of the B₇ = 1 row is 1.

We now turn to gene 8 (BSP). For each of k = 2,3,4,5,6, E_8k is considerably smaller than E_k8, but the reverse is true for k = 7. This suggests that BSP usually gets expressed after genes 2 through 6 and before gene 7. We observed 47 colonies in which both 7 and 8 were expressed, and four colonies in which 7 was expressed but 8 was not. Could the latter four all be errors? The number of errors of this type would have a binomial distribution that is approximately Poisson with mean 50q₈(1 - q₇). The probability that such a distribution has a value of four (or more) is less than 0.05 only if q₈ < 0.03. Our bound, q₈ 0.064, is too weak for this, so we conclude that there is no strong evidence that Pr{T₈ < T₇} = 1.

Finally, we consider gene 9 (OCN). For every k 8, E_k9 is large and E_9k is small. Apparently OCN is usually expressed after the other genes. In fact, E_9k = 0 for k = 1,4,6,7,8, but it is not clear whether these zeroes are due to strict precedence constraints or simply to late expression times of OCN. Nor can we reject the possibility that Pr{T_k < T₉} = 1 for every k 8, with the nonzero values E₉₂, E₉₃, and E₉₅ due to observational errors. Although the case for strict precedence of gene 9 by the two later genes (7 and 8) is stronger than for the earlier genes, we cannot confirm this statistically without more data or smaller bounds on the q_ks.

	ACKNOWLEDGMENT

Top Abstract Introduction Materials and Methods Results Discussion Conclusions Appendix References

 
Supported by operating grants (P.Z. and N.M.), a Postdoctoral Fellowship (A.G.), and an Undergraduate Research Award (Y.Z.) from the Natural Sciences and Engineering Research Council of Canada and an operating grant from the Canadian Institutes of Health Research (J.E.A.; MT-12390).

	REFERENCES

Top Abstract Introduction Materials and Methods Results Discussion Conclusions Appendix References

 

1. Grigoriadis AE, Heersche JNM, Aubin JE. Differentiation of muscle, fat, cartilage and bone from progenitor cells present in a bone-derived clonal cell population: effect of dexamethasone. J Cell Biol 1988;106:2139–2151.[Abstract/Free Full Text]

2. Grigoriadis AE, Heersche JNM, Aubin JE. Continuously growing bipotential and monopotential myogenic, adipogenic, and chondrogenic subclones isolated from the multipotential RCJ 3.1 clonal cell line. Dev Biol 1990;142:313–318.[CrossRef][Medline]

3. Pittenger MF, Mackay AM, Beck SC et al. Multilineage potential of adult human mesenchymal stem cells. Science 1999;284:143–147.[Abstract/Free Full Text]

4. Aubin JE, Triffitt J. Mesenchymal stem cells and the osteoblast lineage. In: Rodan GA, ed. Principles of Bone Biology, Second Edition. New York: Academic Press, 2002:59-81.

5. Aubin JE. Bone stem cells. J Cell Biochem Suppl 1998;30-31:73–82.

6. Nuttall ME, Patton AJ, Olivera DL et al. Human trabecular bone cells are able to express both osteoblastic and adipocytic phenotype: implications for osteopenic disorders. J Bone Miner Res 1998;13:371–382.[CrossRef][Medline]

7. Bellows CG, Aubin JE. Determination of numbers of osteoprogenitors present in isolated fetal rat calvaria cells in vitro. Dev Biol 1989;133:8–13.[CrossRef][Medline]

8. Liu F, Malaval L, Gupta AK et al. Simultaneous detection of multiple bone-related mRNAs and protein expression during osteoblast differentiation: polymerase chain reaction and immunocytochemical studies at the single cell level. Dev Biol 1994;166:220–234.[CrossRef][Medline]

9. Malaval L, Liu F, Roche P et al. Kinetics of osteoprogenitor proliferation and osteoblast differentiation in vitro. J Cell Biochem 1999;74:616–627.[CrossRef][Medline]

10. McCulloch CAG, Strugurescu M, Hughes F et al. Osteogenic progenitor cells in rat bone marrow stromal populations exhibit self-renewal in culture. Blood 1991;77:1906–1911.[Abstract/Free Full Text]

11. Malaval L, Modrowski D, Gupta AK et al. Cellular expression of bone-related proteins during in vitro osteogenesis in rat bone marrow stromal cell cultures. J Cell Physiol 1994;158:555–572.[CrossRef][Medline]

12. Aubin JE. Osteoprogenitor cell frequency in rat bone marrow stromal cell populations: role for heterotypic cell-cell interactions in osteoblast differentiation. J Cell Biochem 1999;72:396–410.[CrossRef][Medline]

13. Bellows CG, Aubin JE, Heersche JNM et al. Mineralized bone nodules formed in vitro from enzymatically released rat calvaria cell populations. Calcif Tissue Int 1986;38:143–154.[Medline]

14. Stein GS, Lian JB, Stein JL et al. Mechanisms regulating osteoblast proliferation and differentiation. In: Rodan GA, ed. Principles of Bone Biology. San Diego: Academic Press, 1996;69-86.

15. Aubin JE. Molecular fingerprinting of osteoblast differentiation: from primitive osteoprogenitor to mature osteoblast. In: Kronenberg HM, ed. The Molecular and Cell Biology of Bone. Tokyo: International Bone and Calcium Institute Incorporation, 1996;54-59.

16. Aubin JE. Regulation of osteoblast formation and function. Rev Endocr Metab Disord 2001;2:81–94.[CrossRef][Medline]

17. Aubin JE, Liu F. The osteoblast lineage. In: Rodan GA, ed. Principles of Bone Biology. San Diego: Academic Press, 1996;51-67.

18. Raetz CR, Wermuth MM, McIntyre TM et al. Somatic cell cloning in polyester stacks. Proc Natl Acad Sci USA 1982;79:3223–3227.[Abstract/Free Full Text]

19. Brady G, Billia F, Knox J et al. Analysis of gene expression in a complex differentiation hierarchy by global amplification of cDNA from single cells. Curr Biol 1995;5:909–922.[CrossRef][Medline]

20. Liu F, Malaval L, Aubin JE. The mature osteoblast phenotype is characterized by extensive plasticity. Exp Cell Res 1997;232:97–105.[CrossRef][Medline]

21. Brady G, Iscove NN. Construction of cDNA libraries from single cells. Methods Enzymol 1993;225:611–623.[Medline]

22. Cox DR, Hinkley DV. Theoretical Statistics. London: Chapman and Hall, 1974.

23. Grimmett G. Percolation. Berlin: Springer-Verlag, 1989.

24. Ogawa M. Differentiation and proliferation of hematopoietic stem cells. Blood 1993;81:2844–2853.[Abstract/Free Full Text]

25. Rickard DJ, Kassem M, Hefferan TE et al. Isolation and characterization of osteoblast precursor cells from human bone marrow. J Bone Miner Res 1996;11:312–324.[Medline]

26. Oizumi K, Yoshiko Y, Bonnelye E et al. Activation of PPAR and , but not PPAR alone, converts osteoprogenitor cells to mature adipocytes in primary rat calvaria cell cultures. J Bone Miner Res 2000;15:S553.

27. Candeliere GA, Liu F, Aubin JE. Individual osteoblasts in the developing calvaria express different gene repertoires. Bone 2001;28:351–361.[Medline]

28. Mayani H, Dragowska W, Lansdorp PM. Cytokine-induced selective expansion and maturation of erythroid versus myeloid progenitors from purified cord blood precursor cells. Blood 1993;81:3252–3258.[Abstract/Free Full Text]

29. Fairbairn LJ, Cowling GJ, Reipert BM et al. Suppression of apoptosis allows differentiation and development of a multipotent hemopoietic cell line in the absence of added growth factors. Cell 1993;74:823–832.[CrossRef][Medline]

30. Hu M, Krause D, Greaves M et al. Multilineage gene expression precedes commitment in the hemopoietic system. Genes Dev 1997;11:774–785.[Abstract/Free Full Text]

31. Liu L, Roberts RM. Silencing of the gene for the ß subunit of human chorionic gonadotropin by the embyronic transcription factor Oct-3/4. J Biol Chem 1996;271:16683–16689.[Abstract/Free Full Text]

32. Phillips RL, Ernst RE, Brunk B et al. The genetic program of hematopoietic stem cells. Science 2000;288:1635–1640.[Abstract/Free Full Text]

This article has been cited by other articles:

				J. E. Aubin Close Encounters of the Bone-Blood Kind IBMS BoneKEy, January 1, 2008; 5(1): 25 - 29. [Full Text] [PDF]

				D. Visnjic, Z. Kalajzic, D. W. Rowe, V. Katavic, J. Lorenzo, and H. L. Aguila Hematopoiesis is severely altered in mice with an induced osteoblast deficiency Blood, May 1, 2004; 103(9): 3258 - 3264. [Abstract] [Full Text] [PDF]

				F. Liu, L. Malaval, and J. E. Aubin Global amplification polymerase chain reaction reveals novel transitional stages during osteoprogenitor differentiation J. Cell Sci., May 1, 2003; 116(9): 1787 - 1796. [Abstract] [Full Text] [PDF]

This Article Abstract