Densely interconnected transcriptional circuits control cell states in human hematopoiesis (original) (raw)

Cell. Author manuscript; available in PMC 2012 Jan 21.

Published in final edited form as:

PMCID: PMC3049864

NIHMSID: NIHMS265972

Noa Novershtern,1,2,3,* Aravind Subramanian,1,* Lee N. Lawton,4 Raymond H. Mak,1 W. Nicholas Haining,5 Marie E. McConkey,6 Naomi Habib,3 Nir Yosef,1 Cindy Y. Chang,1,6 Tal Shay,1 Garrett M. Frampton,2,4 Adam C. B. Drake,2,7 Ilya Leskov,2,7 Bjorn Nilsson,1,6 Fred Preffer,8 David Dombkowski,8 John W. Evans,5 Ted Liefeld,1 John S. Smutko,9 Jianzhu Chen,2,7 Nir Friedman,3 Richard A. Young,2,4 Todd R. Golub,1,5,10 Aviv Regev,1,2,10,‡ and Benjamin L. Ebert1,5,6,‡

Noa Novershtern

1Broad Institute, 7 Cambridge Center, Cambridge MA, 02142

2Department of Biology, Massachusetts Institute of Technology, Cambridge MA, 02140

3School of Computer Science, Hebrew University, Jerusalem, Israel

Aravind Subramanian

1Broad Institute, 7 Cambridge Center, Cambridge MA, 02142

Lee N. Lawton

4Whitehead Institute for Biomedical Research, 9 Cambridge Center, Cambridge, MA 02142

Raymond H. Mak

1Broad Institute, 7 Cambridge Center, Cambridge MA, 02142

W. Nicholas Haining

5Dana-Farber Cancer Institute, Boston, MA 02115

Marie E. McConkey

6Brigham and Women's Hospital, Boston, MA 02115

Naomi Habib

3School of Computer Science, Hebrew University, Jerusalem, Israel

Nir Yosef

1Broad Institute, 7 Cambridge Center, Cambridge MA, 02142

Cindy Y. Chang

1Broad Institute, 7 Cambridge Center, Cambridge MA, 02142

6Brigham and Women's Hospital, Boston, MA 02115

Tal Shay

1Broad Institute, 7 Cambridge Center, Cambridge MA, 02142

Garrett M. Frampton

2Department of Biology, Massachusetts Institute of Technology, Cambridge MA, 02140

4Whitehead Institute for Biomedical Research, 9 Cambridge Center, Cambridge, MA 02142

Adam C. B. Drake

2Department of Biology, Massachusetts Institute of Technology, Cambridge MA, 02140

7Koch Institute for Integrative Cancer Research, Massachusetts Institute of Technology, Cambridge, MA 02139

Ilya Leskov

2Department of Biology, Massachusetts Institute of Technology, Cambridge MA, 02140

7Koch Institute for Integrative Cancer Research, Massachusetts Institute of Technology, Cambridge, MA 02139

Bjorn Nilsson

1Broad Institute, 7 Cambridge Center, Cambridge MA, 02142

6Brigham and Women's Hospital, Boston, MA 02115

Fred Preffer

8Massachusetts General Hospital, Boston, MA 02114

David Dombkowski

8Massachusetts General Hospital, Boston, MA 02114

John W. Evans

5Dana-Farber Cancer Institute, Boston, MA 02115

Ted Liefeld

1Broad Institute, 7 Cambridge Center, Cambridge MA, 02142

John S. Smutko

9Nugen Technologies, San Carlos, CA 94070

Jianzhu Chen

2Department of Biology, Massachusetts Institute of Technology, Cambridge MA, 02140

7Koch Institute for Integrative Cancer Research, Massachusetts Institute of Technology, Cambridge, MA 02139

Nir Friedman

3School of Computer Science, Hebrew University, Jerusalem, Israel

Richard A. Young

2Department of Biology, Massachusetts Institute of Technology, Cambridge MA, 02140

4Whitehead Institute for Biomedical Research, 9 Cambridge Center, Cambridge, MA 02142

Todd R. Golub

1Broad Institute, 7 Cambridge Center, Cambridge MA, 02142

5Dana-Farber Cancer Institute, Boston, MA 02115

10Howard Hughes Medical Institute

Aviv Regev

1Broad Institute, 7 Cambridge Center, Cambridge MA, 02142

2Department of Biology, Massachusetts Institute of Technology, Cambridge MA, 02140

10Howard Hughes Medical Institute

Benjamin L. Ebert

1Broad Institute, 7 Cambridge Center, Cambridge MA, 02142

5Dana-Farber Cancer Institute, Boston, MA 02115

6Brigham and Women's Hospital, Boston, MA 02115

1Broad Institute, 7 Cambridge Center, Cambridge MA, 02142

2Department of Biology, Massachusetts Institute of Technology, Cambridge MA, 02140

3School of Computer Science, Hebrew University, Jerusalem, Israel

4Whitehead Institute for Biomedical Research, 9 Cambridge Center, Cambridge, MA 02142

5Dana-Farber Cancer Institute, Boston, MA 02115

6Brigham and Women's Hospital, Boston, MA 02115

7Koch Institute for Integrative Cancer Research, Massachusetts Institute of Technology, Cambridge, MA 02139

8Massachusetts General Hospital, Boston, MA 02114

9Nugen Technologies, San Carlos, CA 94070

10Howard Hughes Medical Institute

*These authors contributed equally to this work

Supplementary Materials

01.

GUID: 75AE693C-4987-4650-A858-BFED1D01C5B1

02.

GUID: 3513A877-63D5-4CEE-8128-CAA62B038300

03.

GUID: CB63BBA7-51C0-41BE-90FD-A269BF2894AD

04.

GUID: 22CB47EF-7178-416F-9A58-29C97C2926C2

05.

GUID: 123DFECD-8D1B-4FD7-AB53-EE1A51FEA367

06.

GUID: EF4D0AFD-7521-4593-BFBD-5E06252641DA

07.

GUID: 4C4C2EB0-9781-44C4-8E6C-292B606FF24B

08.

GUID: A2550CC5-F029-4538-8AB1-D9ED1AE956C2

Abstract

While many individual transcription factors are known to regulate hematopoietic differentiation, major aspects of the global architecture of hematopoiesis remain unknown. Here, we profiled gene expression in 38 distinct purified populations of human hematopoietic cells and used probabilistic models of gene expression and analysis of _cis_-elements in gene promoters to decipher the general organization of their regulatory circuitry. We identified modules of highly co-expressed genes, some of which are restricted to a single lineage, but most are expressed at variable levels across multiple lineages. We found densely interconnected _cis_-regulatory circuits and a large number of transcription factors that are differentially expressed across hematopoietic states. These findings suggest a more complex regulatory system for hematopoiesis than previously assumed.

Introduction

Hematopoiesis is an ideal model for the study of multi-lineage differentiation in humans. More than 2 × 1011 hematopoietic cells from at least 11 lineages are produced daily in humans from a small pool of self-renewing adult stem cells (Quesenberry and Colvin, 2005). Production of each cell type is highly regulated and responsive to environmental stimuli. Mutations or aberrant expression of regulatory proteins cause both benign and malignant hematologic disorders.

The hematopoietic system is also well-suited for an analysis of the global architecture of the molecular circuits controlling human cellular differentiation. Hematopoietic stem cells, progenitor cells, and terminally differentiated cells can be isolated using flow cytometry. Moreover, many aspects of hematopoietic differentiation can be recapitulated in vitro. Finally, high-speed multiparameter flow cytometry and cDNA amplification procedures allow us to purify and profile gene expression from rare subpopulations (Ebert and Golub, 2004).

A dominant model of hematopoiesis posits that it is controlled by a hierarchy of a relatively small number of critical transcription factors (TFs) that are sequentially expressed, are largely restricted to a specific lineage, and can interact directly to mediate and reinforce cell-fate decisions (Iwasaki and Akashi, 2007). Genetically engineered mice have been used to map the maturation stage at which key TFs are essential (Orkin and Zon, 2008)

Recent genome-wide studies suggest a more complex architecture in regulatory circuits, involving larger numbers of TFs that control different combinations of modules of co-expressed genes (Amit et al., 2009; Suzuki et al., 2009). Complex circuits with a larger number of TFs than previously assumed, each with a major regulatory effect, are emerging from studies in immune cell types (Amit et al., 2009; Suzuki et al., 2009), stem cell populations (Muller et al., 2008) and cell differentiation in invertebrates (Davidson, 2001).

These two views leave open several key questions in understanding the regulatory architecture of human hematopoiesis. (1) Are distinct hematopoietic cell states characterized mostly by induction of lineage-specific genes, or by a unique combination of modules, where the distinct capacities of each cell type are largely determined through the re-use of modules? (2) Is hematopoiesis determined solely by a few master regulators, or does it involve a more complex network with a larger number of factors? (3) What are the regulatory mechanisms that maintain cell state in the hematopoietic system and how do they change as cells differentiate?

Here, we measured mRNA profiles in 38 prospectively purified cell populations, from hematopoietic stem cells, through multiple progenitor and intermediate maturation states, to 12 terminally differentiated cell types (Figure 1). We found distinct, tightly integrated, regulatory circuits in hematopoietic stem cells and differentiated cells, implicated dozens of new regulators in hematopoiesis, and demonstrated a substantial re-use of gene modules and their regulatory programs in distinct lineages. We validated our findings by experimentally determining the binding sites of four TFs in hematopoietic stem cells, by examining the expression of a set of 33 TFs in erythroid and myelomonocytic differentiation in vitro, and by investigating the function of 17 of these TFs using RNA interference. Our data provides strong evidence for the role of complex interconnected circuits in hematopoiesis and for ‘anticipatory binding’ to the promoters of their target genes in hematopoietic stem cells. Our dataset and analyses will serve as a comprehensive resource for the study of gene regulation in hematopoiesis and differentiation.

An external file that holds a picture, illustration, etc. Object name is nihms265972f1.jpg

Hematopoietic differentiation

The 38 hematopoietic cell populations purified by flow sorting and analyzed by gene expression profiling are illustrated in their respective positions in hematopoiesis. Grey - hematopoietic stem cell (HSC1,2), common myeloid progenitor (CMP), megakaryocyte/erythroid progenitor (MEP); Orange - erythroid cells (ERY1-5); Red – CFU-MK (MEGA1) and megakaryocyte (MEGA2); Purple - granulocyte/monocyte progenitor (GMP), CFU-G (GRAN1), neutrophilic metamyelocyte (GRAN2), neutrophil (GRAN3), CFU-M (MONO1), monocytes (MONO2), eosinophil (EOS), and basophil (BASO); Blue - myeloid dendritic cell (DENDa2) and plasmacytoid dendritic cell (DENDa1); Light green - early B-cell (Pre-BCELL2), pro-B-cell (Pre-BCELL3), naïve B-cell (BCELLa1), mature B-cell, class able to switch (BCELLa2), mature B-cell (BCELLa3), and mature B-cell, class switched (BCELLa4); Dark green - mature NK cell (NK1-4); Turquoise - naïve CD8+ T-cell (TCELL2), CD8+ effector memory RA (TCELL1), CD8+ effector memory (TCELL3), CD8+ central memory (TCELL4), naïve CD4+ T-cell (TCELL6), CD4+ effector memory (TCELL7), and CD4+ central memory (TCELL8). See Table S1 for markers information.

Results

An expression map of hematopoiesis reveals cell state-specific profiles

We defined 38 distinct cell states based on cell surface marker expression, representing hematopoietic stem and progenitor cells, terminally differentiated cells, and intermediate states (Figure 1, Table S1). For each state, we purified samples separately from 4 to 7 independent donors by multiparameter flow cytometry (Experimental Procedures), yielding 211 samples. Cells from all stem and progenitor populations were purified from umbilical cord blood. Terminally differentiated lymphocyte populations were purified from peripheral blood, as terminal differentiation is completed in these cells upon exposure to antigens after birth (Table S1). In all cases, cells were harvested fresh, and were processed and sorted immediately. We isolated mRNA from each cell type and measured expression profiles using Affymetrix microarrays (Experimental Procedures).

The global transcriptional profiles are consistent with the established topology of hematopoietic differentiation. Replicate samples from a single state but different donors and samples from multiple states within a lineage are highly correlated with each other, and profiles from related lineages are also similar (Figure 2A). Notably, hematopoietic stem cell (HSC) samples do not form a separate cluster, but are highly similar to early progenitors in the megakaryocyte/erythrocytes lineage (MEGA/ERY), suggesting that their transcriptional state is largely maintained in some of the early progenitors. These findings are also apparent in a systematic unsupervised analysis, using non-negative matrix factorization (NMF) (Brunet et al., 2004) (Figure S1A), and hierarchical clustering (Figure S1B). We further validated our dataset by confirming that previously published lineage-specific gene signatures are significantly enriched in the expected lineage compared to other lineages (FDR < 0.25, Figure S1C, Supplementary Experimental Procedures).

An external file that holds a picture, illustration, etc. Object name is nihms265972f2.jpg

A transcriptional map of hematopoietic differentiation identifies lineage specific transcription

(A) Similarity in global expression profiles between proximate differentiation states. The heat map shows the pair-wise Pearson correlation coefficients between all 211 samples, ordered according to the differentiation tree (right and top). A positive correlation is portrayed in yellow and a negative correlation in purple. (B) Signature genes characterizing the five main lineages. Expression levels are shown for the top 50 marker genes (rows) for each of four major lineages plus hematopoietic stem and progenitor cells. High relative expression is shown in red, and low relative expression in blue; the expression of each gene is normalized to a mean expression of zero across all the samples; labels as in Figure 1. Genes were selected by high expression in one lineage compared to the others (_t_-test). (C) The number of genes that are differentially expressed, according to an outlier statistic, was calculated for all hematopoietic cell states profiled (red), a compendium of 79 tissues in the GNF atlas (Su et al., 2004) (blue) and datasets of lymphomas (Monti et al., 2005) (turquoise), lung cancers (Bhattacharjee et al., 2001) (purple), and breast cancers (Chin et al., 2006) (green). See also Figure S1.

Unique and complex gene signatures characterize distinct hematopoietic lineages

In a supervised analysis we found that each of the five dominant states in our dataset – HSPCs, differentiated erythroid cells, granulocytes/monocytes B cells, and T cells – is distinguished by a set of significantly differentially expressed genes specific to each lineage as compared to the others (Figure 2B, Table S2). Some of these genes are expressed more than 100-fold higher in one cell type (e.g. granzyme genes in certain T-cell and NK-cell populations, PROM1 (CD133 antigen) and HOXA9 in stem and progenitor cells).

The signature genes are enriched for molecular functions and biological processes consistent with the functional differences between lineages (Figure S1D, Table S2). Notably, a set of 16 genes comprised of the 5′ partners of known translocations in leukemias (Mitelman et al., 2010) is enriched in the HSPC population (P<0.013). This suggests that the 5′ partners of leukemia-causing translocations, containing the promoters of the fusion genes, tend to be selectively expressed in stem and progenitor cell populations.

The diversity of gene expression across hematopoietic lineages is comparable to the diversity in gene expression observed across a host of human tissue types. The number of genes that are differentially expressed throughout our hematopoiesis dataset (Outlier Analysis (Tibshirani and Hastie, 2007), Supplemental Experimental Procedures), is comparable to that determined for an atlas of 79 different human tissues (Su et al., 2004), and far higher than in lymphomas (Monti et al., 2005), lung cancers (Bhattacharjee et al., 2001), or breast cancers (Chin et al., 2006) (Figure 2C).

Coherent functional modules of co-expressed genes are re-used across lineages

To dissect the architecture of the gene expression program, we used the Module Networks (Segal et al., 2003) algorithm (Experimental Procedures) to find modules of strongly co-expressed genes, and associate them with candidate regulatory programs that (computationally) predict their expression pattern. We identified 80 gene modules (Figure 3A; modules are numbered arbitrarily by the algorithm), covering the 8968 genes that are expressed in the majority of the samples of at least one cell population. The genes in each of the modules are tightly co-expressed (Figure S2), the 80 modules have largely distinct expression patterns (Figures 3A and S2), and are enriched for genes with distinct biological functions (Figure 3B, Table S3).

An external file that holds a picture, illustration, etc. Object name is nihms265972f3.jpg

Expression pattern and functional enrichment of 80 transcriptional modules

(A) Average expression levels of 80 gene modules. Shown is the average expression pattern of the gene members in each of the 80 modules (rows) across all 211 samples (columns). Colors and normalization as in Figure 2B. The samples are organized according to the differentiation tree topology (top) with abbreviations as in Figure 1. The number of genes in each module is shown in the bar graph (left). The expression profiles of a few example modules discussed in the text are highlighted by vertical yellow lines. The expression of individual genes in each module is shown in Figure S2. (B) Functional enrichment in gene modules. Functional categories with enriched representation (FDR < 5%) in at least one module are portrayed. Categories were selected for broad representation. The complete list appears in Table S3. See also Figures S2 and S7.

A small number of modules are expressed in very specific cell states and reflect the unique functional capacities of a single lineage. For example, Module 889 is expressed in terminal erythroid differentiation and is enriched for genes encoding blood group antigens and organic cation transporters; Module 691 is expressed in B lymphocytes and is enriched for genes encoding immunoglobulins and BCR signaling pathway components; and Module 721 is expressed in granulocytes and monocytes and includes genes encoding enzymes and cytokine receptors essential for inflammatory responses.

Conversely, most modules are expressed at varying levels across multiple lineages, suggesting re-use of their genes in multiple hematopoietic contexts. These include modules expressed in both HSC and progenitor populations (e.g. #865, 679, 805), in both B and T cells (e.g. #673, 703), in both granulocyte/monocyte populations and lymphocytes (e.g. #817, 799, 649), and across all myeloid (e.g. #583) or all lymphoid cells (e.g. #931).

Re-use of modules reflects the differential functional requirements for specific biochemical programs in the various cell states. For example, mitochondrial and oxidative phosphorylation modules (e.g. #847, 583, 883) are induced in erythroid progenitors that produce high levels of heme and are affected most by mitochondrial mutations (Chen et al., 2009; Fontenay et al., 2006), as well as in granulocytes and monocytes, which are capable of a respiratory burst following phagocytosis.

Module states persist through multiple differentiation steps

To delineate the relation between gene expression and differentiation, we projected each module's expression pattern onto the known topology of the differentiation tree (Figures 4 and S4). For example, consider Module 865 (Figures 4A and S3), which is strongly induced in hematopoietic stem and progenitor cells and contains genes encoding key HSPC cell surface markers (CD34 and CD117) and transcriptional regulators (GATA2, HOXA9, HOXA10, MEIS1, and N-MYC). By projecting the module on the differentiation tree, we observe that its induced state in HSCs persists through several consecutive differentiation steps and is repressed at three main points (arrowheads, Figure 4A): (1) after the granulocyte/monocyte progenitor, (2) after erythroid progenitors, and (3) in the differentiation of HSCs towards the lymphocyte lineage.

An external file that holds a picture, illustration, etc. Object name is nihms265972f4.jpg

Propagation and transitions in modules' expression along hematopoiesis

Shown are the mean expression levels of the module's genes in each cell state (colored squares), and selected changes in the predicted regulators, as highlighted in the text (upward arrowhead – regulator induced; downward arrowhead – regulator repressed). Member genes (rather than regulators) in each module encoding TFs are noted below each module, as these may reflect alternative regulators at the same differentiation points. TFs that were validated as regulators of erythroid or granulocyte/monocyte differentiation in a functional assay (Figure 7) are highlighted in bold. The color bar at the bottom of each tree denotes the key lineages, as in Figure 1. (A) HSC and progenitor expression in Module 865; (B) Lineage specific induction in late erythrocytes in Module 727; (C) Lineage specific induction in granulocytes and monocytes in Module 721; (D) Lineage specific induction in B cells in Module 589; (E) One-sided propagation of induced state from HSC to the erythroid lineage in Module 655. (F) Re-use of Module 817, which is inactive in HSCs and independently induced in both lymphoid cells and granulocytes. See also Figures S3 and S4.

We identified a host of such differentiation-associated patterns in gene regulation. One major pattern (31 modules) is HSC persistent states: such modules are active in the HSC state and persist in an active state in several progenitor populations, on either the erythroid/myeloid branch (Figure 4A,E), the lymphoid branch (Figure S4A), or both (Figure S4B,H). The HSC state changes gradually at different points in different modules. Indeed, only Module 631 (Figure S4C) is primarily HSC specific and includes the known stem cell-specific TFs NANOG and SMAD1 (Xu et al., 2008). In other patterns, modules have low or inactive expression in HSCs, but are activated in a single lineage (10 modules), on either the erythroid/myeloid branch (Figure 4B,C, S4D), or the lymphoid branch (Figure 4D). In most cases (39 modules), modules are inactive in HSPCs but are activated in multiple independent lineages (Figure 4F, S4F).

A sequence-based model of the regulatory code

The high degree of co-expression of genes within modules suggests that they may be co-regulated by common transcriptional circuits. We therefore examined each module for enrichment of known and novel _cis_-regulatory elements in their promoters (Supplemental Experimental Procedures). We used six motif-finding methods and a motif clustering pipeline to identify a non-redundant library of enriched elements. We scored each module for the enrichment of each of the novel sites or of known elements or binding events (Sandelin et al., 2004; Subramanian et al., 2005) (Supplemental Experimental Procedures). We identified 156 sequence motifs and 28 binding profiles of 12 TFs (measured by ChIP), that were enriched in the promoters of at least one module (Data available on http://www.broadinstitute.org/dmap/). Of these, 66 are previously unannotated motifs and 118 are associated with 72 TFs (Table S4).

Of these 72 TFs, 11 are known hematopoietic factors (Table S4), and their sites are often enriched in modules consistent with their known functions. For example, the site for the erythroid TF GATA1 (Pevny et al., 1991) is enriched in the late erythroid Module 889, and sites for the lymphocyte regulators Helios and NFATC (Aramburu et al., 1995; Hahm et al., 1998) are enriched in the T and NK Module 559. We also found significant enrichments for TFs with roles in other differentiation processes, which were not previously implicated in hematopoiesis, such as HNF4 alpha (in the HSPC Module 865) and HNF6 (in the lymphoid Modules 859 and 961).

Tightly integrated _cis_-regulatory circuits govern differentiation states

To explore how these _cis_-regulatory associations can give rise to stable cell states, we assembled the regulatory circuits connecting the 276 TFs that were enriched in any gene-set with each other (Figure 5). We connected an edge from each factor with a known motif to all the factors that harbor this motif in their gene promoters (Supplemental Experimental Procedures), and focused only on those factors that were expressed in a given cell state. For example, the circuit of HSC-expressed TFs with known binding sites (Figure 5A) includes many major known regulators of the HSC state (Orkin and Zon, 2008), which are densely interconnected through auto-regulatory (12 of 23 active factors), feedback (15 and 39 loops of size 2 and 3) and feed-forward (206 loops of size 3) loops. Abnormal expression of many of the circuit's TFs is known to cause hematologic malignancies (Look, 1997). This integrated circuitry can give rise to a robust transcriptional network in terminally differentiated cells and HSCs. Notably, since the sequence of the binding site for most TFs is unknown, including 66 of the putative enriched binding sites, the density of regulation is likely even greater than we observed.

An external file that holds a picture, illustration, etc. Object name is nihms265972f5.jpg

Dynamic organization of tightly integrated _cis_-regulatory circuits in HSCs and erythroid cells

(A,B) Shown are _cis_-regulatory networks between TFs (nodes) that are enriched in at least one gene-set, and are expressed (fold-change > 1.5) in (A) HSCs or (B) late erythroid cells. Nodes represent TFs that are active (purple) or non-active (grey) in each of the 4 phases of the erythroid lineage (HSC, MEP, Early ERY, Late ERY). An edge from node a to node b indicates that the promoter of the gene in node b has a binding site for the TF encoded by the gene in node a. Edge colors indicate the Pearson correlation between the expression profiles of the TFs in the connected nodes: red – positive correlation (coefficient > 0.4), black – no correlation (absolute Pearson ≤ 0.4), grey – non-active edge (at least one of the two connected nodes was not expressed in that phase). See Table S4 for enriched motif information.

During the course of differentiation, the HSC circuit gradually disappears along multiple lineages due to loss of expression of the relevant TFs (Figure 5A and data available on http://www.broadinstitute.org/dmap/). Conversely, in terminally differentiated cells, other dense circuits emerge, through the induction of other TFs. For example, the 14 factors in the erythroid circuit include many of the known major regulators of erythroid differentiation (Cantor and Orkin, 2002), including GATA1, LMO2, FOXO4, NFE2 and RXRA (Figure 5B). We find similarly distinct networks in the granulocyte lineage, T cells, and B cells.

Hundreds of transcription factors are differentially expressed across lineages in coherent modules

The dense regulatory circuits between TFs in our sequence-based model suggest that the expression of TF genes is likely to be highly regulated in hematopoiesis. Indeed, supervised analysis finds that many TF genes are strongly differentially expressed in each primary lineage (Figures 6A and S5A), and that the diversity of TF gene expression is comparable between hematopoiesis and the tissue compendium (Su et al., 2004) (Figure S5B).

An external file that holds a picture, illustration, etc. Object name is nihms265972f6.jpg

Lineage specific regulation of TF expression

Signature TF genes with lineage-specific expression in the five main lineages. Shown are the expression levels of the top 50 marker TF genes (rows) selected for each of four major lineages plus hematopoietic stem and progenitor cells (labels as in Figure 1). Genes were selected by high expression in one lineage compared to the others (_t_-test). High expression is shown in red, and low expression in blue; the expression of each gene is normalized to a mean expression of zero across all the samples. See also Figures S5 and S6.

Some TFs are expressed predominantly in a single lineage, including well-studied TFs that are known to be essential for differentiation in HSCs or a particular lineages (Figure S6). However, the expression of those factors often increases gradually along differentiation (Figure S6D,H,I), similar to the gradations observed in gene modules (Figure 4 and S4).

Many other TFs are ‘re-used’ across lineages, either through persistent expression from a single progenitor population or by independent activation in multiple lineages (Figures 4 and S4). For example, Module 793 (Figure S4F), which is induced in both B cells and late erythroid cells, includes several TFs and chromatin regulators. Among these, KLF3 has a reported role in erythroid cells (Funnell et al., 2007), while NFAT5 has a demonstrated function in B cells (Kino et al., 2009).

Many TFs – not previously associated with these lineages – are expressed similarly to known factors, and belong to the same modules, suggesting that the transcriptional circuit consists of a greater number of TFs than previously assumed. For example, the late erythroid Module 727 (Figure 4B) contains four TFs: two are known erythroid TFs (GATA1 and FOXO3A) (Bakker et al., 2007), whereas the others (NFIX1, MYT1) were not previously linked to erythropoiesis. Similarly, the granulocytes/monocytes Module 721 (Figure 4C) contains eight TFs, only two with known roles in the lineage (CEBPA, PU.1/SPI1).

An expression-based model of the regulatory code of hematopoiesis identifies putative regulators controlling changes in differentiation

To identify the potential regulatory role of differentially expressed TFs, we examined the combinations of TFs (regulatory program) which the Module Networks algorithm (Segal et al., 2003) used in order to ‘explain’ the expression of each of the 80 modules (Experimental Procedures). For example, the algorithm associated Module 865 (Figure S3, bottom) with five regulators, most prominently PBX1 (‘top regulator’), SOX4 (‘2nd level regulator’) (Figure S3, top). It predicts that when both PBX1 and SOX4 are induced (in HSCs, CMPs, MEPs, GMPs, early ERY and early MEGA cells) the module's genes are induced too. PBX1 is an established regulator of HSPCs, and SOX4 has recently been shown to be a direct target of HOXB4, a known HSC regulator (Lee et al., 2010), supporting the algorithm's result. The regulators were chosen by their expression alone, and while the model chooses one combination of ‘representative’ regulators, there may be several highly similar TFs that could fulfill the role.

We next interpreted these regulatory connections within the context of the lineage tree., We associated each regulator with the tree positions (Figure 4 and S4, arrowheads), in which a change in the regulator's expression is associated with a change in the module's expression. For example, there are four such positions for PBX1 and SOX4 in Module 865 (Figure 4A, arrowheads), including the association between the repression of PBX1 and the repression of the module in differentiation towards lymphoid lineages (downward arrows, labeled PBX1, Figure 4A). In this way, we predict the roles of distinct TFs at distinct differentiation points, such as MNDA at the granulocyte/monocytes progenitor (Figure 4C, S4G) or NCOA4 and KLF1 at late erythrocytes (Figure S4D).

Overall, the algorithm associated 220 TFs (Table S3) with at least one regulatory program, and 63 TFs as ‘top’ regulators (e.g. Figure S3, top) of at least one module. These include 15 TFs previously associated with hematopoiesis (e.g. TAL1, KLF1, BCL11b, LMO2 and MYB), and 7 associated with differentiation in other systems (e.g. CREG1, MEF2A and NHLH2). For example, we correctly found HOXA9 associated with HSPCs and early erythroid induction (Module 679); NFE2, RXRA, KLF1, and FOXO3 associated with late erythroid induction (Modules 727, 895, 889, and 739, Figure 4B and S5D), HIVEP2 and BCL11b associated with T cell induction (Modules 859, 949, and 667); and HOXC4 and POU2AF1 associated with B cell induction (Module 589, Figure 4D). In addition, the algorithm predicted a regulatory role for proteins that were not previously associated with regulating hematopoietic differentiation (e.g. MNDA, NCOA4). The selected regulators are enriched for TFs that are known to participate as 3′ partners in fusions in hematologic cancers (Mitelman et al., 2010) (25 of the regulators, P<0.028), consistent with a regulatory role in hematopoiesis.

Finally, we compared the predictions of the expression- and sequence-based models. The two models were different, due to two reasons. First, 85% of the TFs chosen as regulators in the expression model (187 of 220) do not have a characterized binding motif in current databases and cannot be identified in the sequence model. Second, 29 of 41 TFs (70%) whose known sites are incorporated in the sequence model, and appear in the expression model, show little or no correlation in expression (absolute Pearson < 0.4) to the module with which they are associated in the sequence model (Data available on http://www.broadinstitute.org/dmap/). Thus, the two models are likely complementary, each capturing a substantial but distinct number of known regulators in the relevant states. To gain confidence in their predictions, we next pursued experimental approaches.

Direct targets of MEIS1, TAL1, IKAROS, and PU.1 in HSPCs reveal dense circuits and anticipatory binding

To validate and further investigate the gene modules and _cis_-circuits, we examined the direct binding of TFs across the genome using chromatin immunoprecipitation followed by sequencing (ChIP-Seq) in HSPCs. We analyzed the binding of MEIS1, TAL1, PU.1/SPI1 and IKAROS/IKZF1, four key regulators of the specification, maintenance or differentiation of HSCs (Argiropoulos et al., 2007; Lecuyer and Hoang, 2004; Ng et al., 2007; Singh et al., 1999), in two replicates, often in independently expanded populations of primary human HSPCs (Supplemental Experimental Procedures). We scored each experiment for statistically significant binding (Supplemental Experimental Procedures, Table S5), and tested each of our expression modules for enrichment in binding events (Table S5).

In modules whose genes are highly induced in terminal differentiation, we found enrichment of binding by corresponding lineage specific factors in HSPCs, suggesting ‘anticipatory’ regulation. For example, Module 727 (Figure 4B), expressed in terminally differentiated erythroid cells, was enriched with target genes bound by in HSPCs by TAL1, an erythroid transcription factor (Table S5). Similarly, genes in the granulocyte/monocyte Module 763 were enriched for targets bound by PU.1 in HSPCs (Table S5); and genes in the lymphoid Module 949 were enriched for target genes bound by IKAROS in HSPCs (Table S5). In many (but not all) cases, expression of the target module is already moderate in HSCs, and increases with differentiation. This strongly supports an ‘anticipatory’ regulation, where relevant differentiation TFs are bound at target promoters in HSPCs, resulting in mild expression of targets which persists and further increases upon differentiation.

Some of our expression-based model's predictions for HSPCs are supported by the ChIP-Seq data. For example, the two modules that are induced in HSPCs and are associated in our model with either MEIS1 (Module 961) or its known binding partner PBX1 (Module 865, Figure 4A) are enriched in target genes bound by MEIS1. MEIS1 and HOXA9 are members of Module 865 consistent with MEIS1's auto-regulatory binding (Table S5). The ChIP-Seq data also support module re-use. For example, several of the modules enriched with PU.1 are ‘re-used’ in granulocytes and B lymphoid cells (e.g., Modules 853. 649, 979, 769 and 817), consistent with an established role for PU.1 in both lineages. In other cases module re-use may be mediated by combinatorial binding of two factors (e.g. by both PU.1 and IKAROS in Module 607 which is expressed in granulocytes, monocytes and some lymphoid cells).

The individual binding events in our profiles also support the overall organization observed in the _cis_-circuits in the sequence model. First, three of the factors bind their own promoter (IKAROS and MEIS1) or enhancer (PU.1), forming auto-regulatory loops, as observed for many known master regulators (Boyer et al., 2005) and in our sequence model. Second, PU.1, IKAROS and MEIS1 are integrated in a feed-forward loop. Third, there is a significant overlap between the targets of any pair of factors (Table S5). Finally, in aggregate, the factors bind 13 of the 23 other TFs in our HSC circuit, further increasing its density.

Differential expression of candidate transcription factors during in vitro differentiation

We confirmed the lineage-specific expression of 33 TFs in primary human hematopoietic progenitor cells induced to differentiate in vitro. We focused on the erythroid and myelomonocytic lineages, as differentiation of primary human hematopoietic progenitor cells can be faithfully recapitulated and genetically manipulated along these lineages in vitro. We selected a set of 33 TFs identified in either the sequence or gene expression-based models as candidate regulators of these two lineages.

We developed a quantitative, multiplexed assay to detect the expression of the signature genes in a single well using ligation mediated amplification (LMA) followed by amplicon detection on fluorescent beads (Peck et al., 2006). We cultured primary human CD34+ cells from adult bone marrow in vitro in cytokine conditions promoting either erythroid or myelomonocytic differentiation. We harvested cells at 12 time points between days 3 and 10 of erythroid and myelomonocytic differentiation and determined TF gene expression using the multiplexed bead-based assay. We confirmed that the 33 TFs are differentially expressed between the two lineages, providing a robust expression signature that can distinguish between the two states independent of profiling platform, in cells derived from adult bone marrow or umbilical cord blood, and in cells that differentiated in vivo or in vitro (Figure 7A).

An external file that holds a picture, illustration, etc. Object name is nihms265972f7.jpg

Experimental validation of 33 TFs

(A) The expression of 33 TFs was detected in primary human bone marrow CD34+ progenitor cells undergoing differentiation in vitro, harvested at 12 time points between days 3 to 10 of differentiation, and detected by a multiplexed assay using LMA followed by fluorescent bead-based detection (left heat map). In the heat map in the right panel, the expression of the same TFs in the original Affymetrix dataset is illustrated. The labels at the far left indicate whether the TF was chosen as a regulator in the expression-based model or in the sequence-based model. (B) Differentiation following TF silencing with shRNA. Human bone marrow CD34+ cells expressing shRNAs targeting TFs were induced to differentiate in-vitro for 10 days, and the ratio of erythroid (glycophorin A positive) and myelomonocytic (CD11b positive) cells was measured by flow cytometry. Each black dot represents an individual shRNA (mean of three replicates), and bars indicate their average. The effect of a control shRNA targeting the luciferase gene, which is not expressed in human cells, is indicated with a dashed line. Below the shRNA labels, * or ** indicates p < .05 for one or both shRNAs, respectively. Bottom – classification of the TFs according to their roles in the expression-based and sequence-based models, and to their induction pattern in the LMA profiling. (C) The effects of additional shRNAs targeting candidate TFs expressed in CD34+ cells derived from both umbilical cord blood and adult bone marrow and assayed as in panel B (* indicates p < .01). (D) Representative flow cytometry scatter plots from shRNAs expressed in umbilical cord blood. See additional information in Tables S5, S6 and S7.

Changes in expression levels in transcription factor circuits functionally modulate differentiation in vitro

We next tested whether acute loss of expression of each TF using RNA interference can functionally affect erythroid and myelomonocytic differentiation. We used our multiplexed bead-based assay to identify short hairpin RNAs (shRNAs) that effectively knockdown each TF, and found 17 genes, with at least two different effective shRNAs. Next, we infected primary human adult bone marrow CD34+ cells with the validated lentiviral shRNAs, cultured the cells in cytokine conditions supporting both erythroid and myelomonocytic differentiation, and assessed the number of erythroid (glycophorin A positive) cells relative to myelomonocytic (CD11b positive) cells by flow cytometry (Figure 7B). In most cases, the shRNA perturbation dramatically altered differentiation, with the ratio of erythroid to myeloid cells ranging from less than 1:10 to more than 10:1 with different shRNAs.

The perturbations associated with the lowest fraction of erythroid cells in culture corresponded to the samples expressing shRNAs targeting 9 TFs expressed at higher levels in the erythroid lineage (Table S6). Consistent with our models, six were regulators in either the expression or the sequence model, and the other three members of erythrocyte induced modules (Figure 7B, bottom). These include GATA-1 and KLF1, TFs with well-established roles in erythroid differentiation, (Funnell et al., 2007; Pevny et al., 1991), and TAL1 and FOXO3A, which have been implicated in erythroid differentiation (Aplan et al., 1992; Bakker et al., 2007). The TF YY1 was identified in our sequence-based models, has higher expression in erythroid cells, and was functionally validated by our shRNA screen. A physical association between YY1 and GATA-1 was reported in the chicken α-globin enhancer (Rincon-Arano et al., 2005). Finally, we validated a new role for HIF3A and AFF1 (AF4) in erythroid differentiation based on module membership and perturbation. Notably, AFF1 is a common translocation partner with the MLL gene in leukemia (Li et al., 1998).

Conversely, 8 perturbations resulted in the lowest fraction of myelomonocytic cells and corresponded to samples expressing shRNAs targeting 7 TFs induced in granulocyte/monocyte cells and one (E2F1) with higher expression in erythroid cells. Four TFs were predicted by the expression model to regulate modules induced in granulocytes/monocytes and five were predicted in the sequence network (Figure 7B, bottom). These include the well-established granulocyte/monocyte TFs, PU.1/SPI1 and C/EBP family members (Hirai et al., 2006; Scott et al., 1994), and VDR, a gene that has been implicated in myeloid differentiation (Liu et al., 1996).

We further validated three TFs that had not previously been associated with erythroid differentiation (AFF1, HIF3A, and YY1) alongside a known erythroid regulator (FOXO3A) and a known granulocyte regulator (PU.1/SPI1) (Figure 7C,D). We tested additional shRNAs for each gene by quantitative PCR and identified two shRNAs per gene that decrease expression of their target genes by 63 to 95% in human CD34+ cells derived from both adult bone marrow and umbilical cord blood. Using flow cytometry for lineage-specific markers following 10 days of differentiation, we validated our initial findings that AFF1, HIF3A, and YY1 decrease the relative production of erythroid lineage cells. These results are further supported by profiling mRNA levels following knockdown of these five TFs at 4 days following lentiviral infection. Compared to a control shRNA, gene expression in cells expressing AFF1, HIF3A, and FOXO3A were anti-correlated with erythroid profiles and positively correlated with granulocytes (Figure S5C) and knockdown of PU.1/SPI1 had the inverse pattern, as expected. Knockdown of YY1 caused a transcriptional profile more similar to HSCs, indicating a more substantial block in terminal differentiation. Taken together, our findings indicate that modulating the expression of TF genes can powerfully alter hematopoietic differentiation.

Web-based portal as a research resource

To facilitate interrogation of our hematopoietic gene expression database by the broader scientific community, we have created a Web-based portal (http://www.broadinstitute.org/dmap), providing access to the primary data, sample information, processed results from both models, and a suite of analytic tools.

Discussion

General principles of transcriptional circuits in differentiation

The changes in gene expression over the course of hematopoietic differentiation are profound. The number of differentially expressed genes is similar within hematopoiesis and across human tissues, suggesting comparable complexity. Our findings reveal several major principles about the organization of this transcriptional program.

Gene expression in hematopoiesis can be decomposed into modules of tightly co-expressed genes, some of which are restricted to specific lineages, whereas most are re-used in multiple lineages. Furthermore, modules' transcriptional state persists through multiple differentiation steps. For example, the transcriptional state of HSCs is not switched off immediately, but instead persists with gradually decreasing expression in progenitor cells.

Many of the TFs with known binding sites can be assembled into densely-interconnected circuits. These can provide a mechanism for robust gene regulation in both terminally differentiated cells and HSCs. Since the binding sites for many factors remain unknown, we therefore expect that the circuit's density and complexity is even higher.

A large number of TFs are differentially expressed across hematopoiesis, often in tightly co-regulated modules, and at comparable complexity to that of the other (non-regulatory) genes. Leveraging this correspondence, we associated TFs to the modules and differentiation states that they may regulate. We automatically re-discovered (without prior knowledge) many of the key known TFs, and predict regulatory functions for numerous additional TFs.

By monitoring the binding of four major TFs in HSPCs, we found that ‘anticipatory’ regulation may be a major feature of these circuits. In such cases, TFs that direct lineage-specific differentiation bind a significant portion of their target genes in HSPCs. These target genes are often moderately expressed in the stem and progenitor cells, with substantial further induction as differentiation progresses in the relevant lineage. This is consistent with the concept of ‘lineage priming’ in HSCs (Akashi, 2005), providing flexibility in cell fate commitments.

Discovering and validating transcriptional regulators in hematopoiesis

Our examination of the global architecture of hematopoietic differentiation offers a complementary strategy to studies of individual genes in murine models. In this approach, gene expression and sequence-based analyses nominate a host of candidate regulators, and point groups of factors that may act together and hence introduce redundancies. The two computational approaches complement each other: the expression model may identify factors whose binding specificity is unknown; the sequence model may help detect those factors whose mRNA levels do not change or do not correspond to changes in targets (Lu et al., 2009).

We used a perturbation-based approach to validate TFs derived from the sequence and expression models in an in vitro differentiation system. Modulating expression of candidate TFs with RNA interference altered differentiation of hematopoietic progenitor cells in vitro in the direction predicted by our models. We re-confirmed the role of several known factors, and identified several new ones (e.g. YY1, AFF1, HIF3A). In vitro manipulation can be more sensitive than genetic ablation experiments in vivo, where perturbations may be corrected by homeostatic mechanisms, such as cytokine or transcriptional feedback loops.

A transcriptional roadmap for hematologic malignancies

Balanced translocations involving TFs play a major role in the pathogenesis of human leukemias. Of 200 known translocations in AML (Mitelman et al., 2010), there are 53 in which at least one translocation partner is a TF, 16 as a 5′ partner and 43 as a 3′ partner (6 in both). 25 of these 43 known 3′ partners are among the 220 regulators in the expression model (P<0.028), and 5 of the known 3′ partners are among the 72 known TFs in the sequence model. Furthermore, 5′ partners are enriched in genes expressed in HSPCs. These results support the role of chosen regulators in differentiation, and are consistent with a broader paradigm in which lineage-specific promoters can dysregulate key TFs to disrupt differentiation (Rosenbauer and Tenen, 2007; Tomlins et al., 2005).

Impaired or blocked hematopoietic differentiation is a defining characteristic of leukemia, and the gene expression profiles of leukemias cluster strongly into subgroups that correspond to specific molecular subgroups (Bullinger et al., 2004; Tamayo et al., 2007; Valk et al., 2004). Gene signatures induced in various leukemias significantly overlap those induced in normal hematopoiesis (Figure S7A,B). In most cases there is a coherent overlap between the leukemia subtype and the cell type from which it is known to arise. However, human leukemias often have more complex combinations of modules that are not observed in normal samples, including HSPC modules, as has been reported in murine models of leukemia (Krivtsov et al., 2006).

Towards a programming ‘code’ of hematopoietic differentiation

A more complete understanding of hematopoietic differentiation will likely require an integration of gene expression data with other genomic data, including epigenetic analyses, genome-wide ChIP-Seq studies, proteomics, and systematic functional studies. Given the ability to produce high quality measurements from small numbers of cells, gene expression data provides a first draft of the transcriptional program controlling hematopoiesis, opening the way to manipulate and reprogram these circuits, through perturbation and manipulation of each regulatory factor. These can highlight avenues for therapeutic intervention, including ‘reprogramming’ of cells to more desired states. While deriving mechanistic insights from mammalian gene expression profiling data has been challenging, hematopoiesis provides a paradigm for the testing of more advanced algorithms (Kim et al., 2009). Our dataset and analyses provide a resource for further inquiries into normal and pathologic hematopoietic differentiation in humans.

Experimental Procedures

Subjects and samples

Human umbilical cord blood was harvested from postpartum placentas at Brigham and Women's Hospital under an Institution Review Board (IRB)-approved protocol. Peripheral blood samples were obtained from healthy volunteers at the Dana-Farber Cancer Institute with informed consent under an IRB-approved protocol.

The majority of cells were purified from umbilical cord blood, an enriched source of undifferentiated populations. However, terminally differentiated lymphocyte populations including T-cells (TCELL1-8), B-cells (BCELL1-4), natural killer cells (NKa1, NKa2, NKa3, NKT) and dendritic cells (DENDa1, DENDa2) were purified from adult peripheral blood because terminal differentiation in these populations requires exposure to antigens after birth. For each cell population, we purified samples from 4-7 distinct donors. All blood samples were harvested fresh and immediately processed for flow sorting.

Cell sorting strategy and flow cytometry

First, mononuclear cells were isolated by Ficoll-Hypaque sedimentation. For relatively rare populations, including hematopoietic stem cell populations (HSC1 and HSC2), progenitor populations such as common myeloid progenitor (CMP), megakaryocyte/erythroid progenitor (MEP), granulocyte-monocyte progenitor (GMP), and for the erythroid lineage populations (ERY1-5), lineage depletion was performed using antibodies against CD2, CD3, CD4, CD5, CD8, CD11b, CD14C, CD19, and CD56 with a magnetic column (Miltenyi Biotec, Auburn, CA). Positive selection was then performed using flow cytometry for labeled antibodies to the markers described in Table S1. For the more common or terminally differentiated populations including neutrophil populations (GRAN1-3), basophils (BASO1), monocytes (MONO1-2), eosinophils (EOS2), megakaryocytes (MEGA1-2), B-Lymphoid Progenitor (PRE_BCELL1), Pro and Early B-lymphocytes (PRE_BCELL2, PRE_BCELL3), dendritic cells (DENDa1, DENDa2), mature T-cells (TCELL1-8), mature B-cells (BCELL1-4) and natural killer cells (NKa1-3, NKT), cells were positively selected using flow scatter properties and antibodies based on the immunophenotypes described in Table S7. The gene expression profiles for a subset of the lymphoid populations, has been analyzed previously (Haining et al., 2008).

Sorting was performed with Vantage SE Diva or FACSAria flow cytometers (BD (Becton-Dickinson, San Jose, CA). Cell populations of interest were collected into tubes containing PBS in a collection unit at 4°C. The >95% purity of populations was confirmed by performing FACS analysis of the sorted cells. Sorted cells were spun down, immediately re-suspended in TriZol (Invitrogen, San Diego, CA) and stored at -70°C.

Microarray data acquisition

Total RNA was isolated from TriZol. The concentration of RNA was quantified using the RiboGreen© RNA Quantitation Kit (Invitrogen, San Diego, CA). Ten nanograms of total RNA were amplified using the Ovation Biotin RNA Amplification and Labeling System (NuGEN, San Carlos, CA). The cDNA was fragmented, labeled and hybridized to Affymetrix HG_U133AAofAv2 microarrays (Affymetrix, Santa Clara, CA), which contains 22,944 probes. Transcript levels were normalized using RMA, genes where filtered and normalized again such that the mean of each gene is zero. Dataset is available on http://www.ncbi.nlm.nih.gov/geo/, GSE24759.

Expression-based Module Networks model

The modules and their regulation programs were automatically learned using the Module Networks algorithm (Segal et al., 2003). This method detects modules of co-expressed genes and their shared regulation programs. The regulation program is a small set of genes, whose expression is predictive of the expression level of the module genes using a decision (regression) tree structure. Given the expression values and a pool of candidate regulator genes, a set of modules and their associated regulation programs are automatically inferred by an iterative procedure. This procedure searches for the best gene partition into modules and for the regulation program of each module, while optimizing a target function. The target function is the Bayesian score, derived from the posterior probability of the model (see (Segal et al., 2005) for a detailed description of the algorithm).

Supplementary Material

01

02

03

04

05

06

07

08

Acknowledgments

We thank E. Lander, I. Amit and I. Gat-Viks for critical review of the manuscript; D. Scadden for helpful discussions; L. Gaffney and S. Hart for assistance with figure generation; and D. Peck, J. Lamb, R. Onofrio, and the Broad Genetic Analysis Platform for assistance with expression arrays. The work was funded by the NIH (grants R01 HL082945 and P01 CA108631 to BE, and the PIONEER award to AR), the Burroughs-Wellcome Fund (CAMS to BLE and CASI to AR), funds from Landon and Lavinia Clay (RAY) and HHMI (TRG and AR). AR is an investigator of the Merkin Foundation for Stem Cell Research at the Broad Institute.

Footnotes

Further details for data analysis, chromatin immunoprecipitation, and functional validation experiments are described in the Supplemental Experimental Procedures.

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

References