development and differentiation are characterized by genetic circuitry or gene regulatory networks ( grns ) that have inherent forward momentum encoded by a number of regulatory motifs ( davidson , 2010 ) . to self - renew and maintain differentiation potential , stem cells must structure their grns so as to arrest or buffer this forward trajectory . networks at early multipotent stages may bear little relation to those of mature differentiated cells , making comparison between them difficult . detailed time - series data are useful in this regard ( bruno et  al . , 2004 ; may et  al . , 2013 ) , but there is likely substantial asynchrony between individual cells at any given time point . cells also may undergo lineage commitment through different initial gene expression trajectories ( pina et  al . , 2012 ) . together , these factors can confound attempts to infer network architectures and gain molecular insights into commitment and subsequent lineage specification from averaged gene expression profiles . analysis of gene expression in single cells offers a different approach that , in principle , makes use of cellular heterogeneity as a source of variation for establishing gene - gene associations . 2013 , 2015 ) or inferred pairwise gene associations using only co - expressing cells ( sthlberg et  al . , 2011 ) . it has been suggested that levels of expression are better accounted for in cells that co - express both genes , and may be obscured by presence / absence effects when all cells are considered ( rusnakova et  al . , 2013 ) . we have tried to address this constraint in our exploration of gene expression networks around the erythroid versus myelo - monocytic lineage choice of multipotent hematopoietic progenitor cells . additionally , we have focused on capturing networks from closely related cells in the vicinity of the commitment boundary , to gain insight into the evolution of grns relevant to lineage specification . the erythro - myeloid bifurcation is an intensively studied paradigm , and both transcription factors ( tfs ) and regulatory motifs involved in physiological transcriptional programming of these alternative fates have been described ( wolff and humeniuk , 2013 ) . key players include gata factors , the ets family protein pu.1 , and c / ebp family members , whose potency has been demonstrated in the experimental reprogramming of blood lineages ( graf and enver , 2009 ) . we used two distinct cell commitment scenarios to obtain high resolution around the early phase of commitment and lineage specification . first , we identified and prospectively isolated cells spontaneously committing to both lineages under culture conditions that maintain self - renewing ( sr ) cells and lack pro - differentiative cytokines . second , we used inducible variants of the lineage - affiliated tfs gata1 and pu.1 to drive cells into lineage specification , again in the absence of pro - differentiative cytokines . this approach allows timed sampling after instigation of a discrete lineage trigger , and it may provide a more homogeneous molecular entry into commitment and lineage development . we describe state - distinct networks in multipotent and early lineage - committed cells and , in particular , highlight the existence of lineage - conflicting programs at the emergence of lineage choice . we further identify an axis involved in lineage specification that includes gata2 and ddit3 , a c / ebp family member previously implicated in stress response ( zinszner et  al . , 1998 ) and described as a potential target of erythropoietin signaling in erythro - leukemic cell lines ( coutts et  al . , 1999 ) . we have explored commitment in the non - transformed bone marrow ( bm)-derived hematopoietic multipotent cell line fdcpmix . fdcpmix is karyotypically normal , il-3 dependent , and capable of multilineage differentiation in response to the appropriate environmental cues . under maintenance culture conditions , sr and lineage - committed cells ( erythroid- or myeloid - committed progenitors [ ecps or mcps ] ) co - exist in culture ( figure  1a ) and can be isolated on the basis of their surface phenotype ( figure  1b ) . sr cells are kit+gr1 cells ( figure  1b ) with proliferative capacity in bulk and clonal cultures ( figures s1a and s1b ) and are uniquely able to faithfully reconstitute the cellular heterogeneity observed in maintenance cultures ( figure  s1c ) . lineage - committed cells devoid of sr potential are gr1 + mcps , with an early myeloid morphology ( figure  s1d ) and no erythroid differentiation capacity ( figure  s1e ) , and kitgr1 ecps , with accelerated erythroid differentiation ( figure  s1e ) but minimal or no contribution to neutrophil cultures ( figure  s1f ) . the transcriptional signatures of sr , mcp , and ecp compartments are readily distinct ( figure  s1 g ) and confirm their lineage affiliation ( table s1 ) . we next explored cellular heterogeneity within these cell compartments using single - cell multiplex qrt - pcr . the results showed substantial cell - to - cell heterogeneity ( figure  1c ) and overlap ( figure  1f ) within and among all three compartments ; nevertheless , the compartments may be robustly identified using single and dual - gene classifiers ( table s1 ) . the transcriptional heterogeneity observed within individual populations may highlight the capture of multiple contemporaneous molecular programs underlying lineage commitment under self - renewal conditions . since tfs are potent instigators of lineage determination or reprogramming , we next examined cellular and transcriptional heterogeneity in fdcpmix cells driven to erythroid or myeloid lineages through expression of inducible gata1 or pu.1 estrogen receptor fusions , respectively ( figure  s1h ) . we analyzed the transcriptional programs of single cells captured at various time points ( figure  s1i ) after induction . in parallel , we studied the temporal dynamics of lineage commitment in this setting and functionally tested their commitment status by evaluating the following : ( 1 ) their retention of self - renewal potential , i.e. , their capacity to re - initiate maintenance cultures ; and ( 2 ) their lineage potential in response to various cytokine cues . this experimental design affords a dynamic appreciation of cellular and molecular mechanisms employed in lineage specification . activation of gata1 in sr cells led to phenotypic changes at as early as 46  hr ( figure  1d ) , accompanied by the loss of neutrophil potential ( at 6  hr ) and followed by the loss of sr capacity ( at 45  hr ) ( figure  1e ) . enforced pu.1 activity resulted in the loss of clonogenic sr potential and elicited a myeloid differentiation bias in a more extended time frame ( figures s1j and s1k ) . single - cell transcriptional profiling during tf - driven commitment confirmed the lineage identity of the cells obtained , which broadly separated away from the sr state along erythroid ( gata1 ) and myeloid ( pu.1 ) axes ( figure  1f ) . gata1-driven cells co - occupied the same transcriptional space as ecp cells and showed a similar extent of cell - to - cell variation . in contrast , pu.1-driven cells , while similarly heterogeneous to mcps , appeared to occupy a distinct territory ( data not shown ) . this presumably reflects the neutrophilic status of mcps and a monocytic bias of pu.1-ert - differentiated cells , consistent with prior reports of pu.1-driven cell fate ( laslo et  al . , 2006 ) . we exploited the heterogeneity of cells at early stages of spontaneous and tf - driven lineage commitment to explore the transcriptional networks controlling lineage specification . inspection of patterns of expression for individual genes revealed a fraction of cells in which the gene is off , and a fraction of cells expressing the gene to varying levels ( on ) ( figure  2a ) . the on / off status can be described as binary while the distribution of on values represents a continuous component of the data . thus , for any given gene pair , both binary and continuous relationships are possible ; this is exemplified for gata1 and epor in figure  2b . methodologically , we used odds ratio ( or ) to quantify on / off gene - to - gene associations ( figure  2c ) and spearman rank correlation to measure correlations between gene expression levels ( figure  2d ) . we combined gene associations obtained by both methodologies to infer putative regulatory networks characterizing sr states and the different modes of lineage commitment ( figure  2e ) . at coarse grain , the networks revealed increased connectivity in the lineage - committed compared to the sr state . also , commitment appeared associated with a higher frequency of negative associations between genes ( figures 2f and s2a ) , including known lineage - determining factors ( table s1 ) . while this may be , to some extent , a function of the genes analyzed , it also may reflect mechanistically distinct processes governing acquisition of lineage identity versus exit from self - renewal . negative associations are less prominent in the full - activation time courses of tf - driven commitment , as the networks capture not only the early processes of lineage specification , but also the later consolidation of the differentiation program . in contrast , detailed temporal analysis of gata1-ert - driven lineage specification revealed that cross - antagonistic associations between lineage determinants peak at 6  hr ( figure  2 g ) , coincident with early loss of neutrophil differentiation potential ( figure  1e ) en route to lineage commitment , suggesting that resolution of lineage conflicts is an early step in acquisition of lineage identity . in this respect , ddit3 emerges as an interesting candidate in lineage cross - antagonism : it is positively associated with gata2 , both in sr and committed cells ( figure  s2b ) , and negatively associated with the neutrophil determinant cebpa , either directly or through gata2 , in mcps and at early stages of erythroid lineage commitment ( figure  s2b ; table s1 ) . since ddit3 has not previously been tabled as a central regulator of erythro - myeloid lineage specification , we functionally tested its impact in loss- and gain - of - function experiments . knockdown of ddit3 ( figures s3a and s3b ) in fdcpmix cells resulted in the loss of erythroid and mixed - lineage colonies , with no change to myelo - monocytic potential ( figure  3a ) . the same loss of erythroid potential in colony - forming assays was observed in stem and progenitor cells ( kls ) from mouse bm upon knockdown ( figure  3b ) and constitutive knockout ( figure  3c ) of ddit3 expression . the data are compatible with a requirement for ddit3 in the erythroid lineage , while it is dispensable for the development of the myeloid lineage . the negative association observed between ddit3 and cebpa in inferred transcriptional networks from early stages of lineage specification ( figures 2e and s2b ; we tested this in myeloid - committed granulo - monocytic progenitors ( gmps ) by enforcing ddit3 expression ( figure  s3c ) , resulting in a transient re - acquisition of self - renewal potential ( figure  3d ) and a dramatic change in the nature of the colonies obtained ( figure  3e ) , with the predominance of large gm colonies of immature appearance ( figure  3f ) . cells in these colonies expressed immature surface  markers and were predominantly lineage - negative kit+cd34+cd16/32 + , thus presenting an essentially gmp phenotype albeit with variable levels of sca1 expression ; in contrast , cells in control colonies exhibited a differentiated gr1+mac1 + phenotype ( data not shown ) . taken together , the data suggest that ectopic expression of ddit3 in gmps blocks lineage progression and transiently re - activates self - renewal capacity . we used single - cell gene expression profiling of gmps , either wild - type or transduced with a control vector or a ddit3-expressing lentivirus , to interrogate the transcriptional program changes imposed by enforced expression of ddit3 and to inspect its role in remodeling of the transcriptional networks underlying lineage progression and/or identity . enforcement of ddit3 changed the expression of two - thirds of genes ( figures s3d and s3e ) predicted as its neighbors in our inferred transcriptional networks ( figure  s2b ) , attesting to the robustness of our inference approach . principal component analysis ( pca ) of the populations of individual wild - type and transduced gmp cells separated ddit3-expressing cells from controls ( figure  3 g ) . this separation is mostly attributable to the increased expression of early erythroid regulators gata2 , tal1 , zfpm1/fog1 , and gfi1b and the relative loss of m , gm , and g - csf receptors ( csf1r , csf2ra , and csf3r , respectively ) as well as of c / ebp family members ( figure  3h ) . the relative gain in the expression of erythroid - affiliated genes and loss of myeloid csf receptors and c / ebp family tfs further developed with prolonged expression of ddit3 in gmps under differentiation conditions that support multilineage output ( figures s3f and s3 g ) . these data confirm ddit3 as a positive regulator of erythroid lineage specification at the expense of myeloid fate , providing an experimental validation of the predictive power of the networks we derived by analyzing the heterogeneity of single cells undergoing lineage specification . we next asked if the relative gain in importance of erythroid - affiliated regulators in ddit3-transduced gmps was associated with a global remodeling of the transcriptional networks underlying gmp lineage identity . indeed , we observed an overall loss in network connectivity specific to the activity of ddit3 ( figure  4a ) moreover , there was a relative gain in connectivity of gata2 at the expense of myeloid hubs , as quantified in figure  s4a . for a broader appreciation of the transcriptional changes induced by ddit3 , we performed rna sequencing ( rna - seq ) . gmps transduced with ddit3 or control vector ( csiem ) were cultured for up to 5  days under conditions supportive of multilineage output . similarly to cells obtained from colony - forming assays , ddit3-transduced cells retained a gmp - like phenotype , while control cells acquired differentiated myeloid surface markers ( figure  s4b ) . the global transcriptional profiles of cells with enforced expression of ddit3 ( figure  s4c ) were clearly distinct from control - transduced and wild - type gmps . gene set enrichment analysis ( gsea ) showed that ddit3 expression is associated with global loss of gmp programs and concomitant upregulation of meg - e - affiliated signatures ( figure  4b ) . interestingly , signatures representative of pre - gm cells , the developmental precursors of gmps ( figure  4c ) , also were upregulated ( figure  4b ) . these data suggest that ddit3 acts through the development or stabilization of a more primitive precursor with inherent meg - e potential ( figure  4c ) . analysis of the networks derived from wild - type gmps ( figure  s4d ) as well as ddit3-transduced cells exposed to conditions supportive of multilineage output for 2  days ( figure  4a ) revealed increased importance of specific erythroid versus myeloid regulatory nodes . given the association of gata2 and ddit3 seen in our network analysis of lineage commitment ( figure  2e ) and the increase in activity of gata2 as a hub in ddit3-transduced gmps ( figure  4d ) , we explored the behavior of gata2 target genes , previously identified by chromatin immunoprecipitation sequencing ( chip - seq ) in fdcpmix ( may et  al . , 2013 ) , in response to ddit3 expression . gsea on the rna - seq data from ddit3 and control vector - transduced cells provided evidence for a coincidence of gata2 and ddit3-driven gene expression programs ( figure  4e ) . to further explore the ddit3-gata2 axis in lineage specification , we focused on gene expression programs resident in those single cells that co - expressed the two genes , and derived networks from multipotent and lineage - committed compartments ( table s1 ) . within the ecp compartment , this analysis revealed the following : ( 1 ) an association between the two genes that was not seen when cells were used irrespective of ddit3/gata2 status , and ( 2 ) the involvement of both these genes in anti - correlations with myeloid - affiliated regulators . a similar picture now also emerges in multipotent cells , highlighting the gains in network information obtainable from interrogating single cells selected on the basis of specific co - expression patterns . our analysis of cells in the vicinity of commitment highlights the existence of lineage - conflicting programs at the emergence of lineage choice , and it identifies the importance of a ddit3-gata2 axis in this process . in revealing this association , it was instrumental that our network inference approach acknowledged and took account of both binary and continuous components of single - cell gene expression . as such , we captured both on / off relationships and the associations seen between distributions of on expression levels in co - expressing cells . furthermore , having established potentially interesting associations , we zoomed in on single cells where the associations were present to refine their specific network identity . in the case of ddit3-gata2 , this supported the importance of cross - antagonistic interactions in lineage decisions . while in principle this approach can be applied to any gene - gene interaction in an unsupervised manner , it would require profiling of considerably higher numbers of individual cells to ensure statistical robustness and systematically explore all possible interactions . to validate our inferred networks , we elected to experimentally test the predicted role of ddit3 in erythroid lineage specification . functional experiments revealed an early erythroid effect from ddit3 loss of function produced by knockdown or genetic knockout . however , we note that the constitutive ddit3 knockout mice used in this study did not display any discernible erythroid defect ( c.f . , unpublished data ) . thus , lineage - determining networks may be robust to ddit3 deletion or constitutive deletion may be compensated during development . forced ddit3 expression does promote erythroid programming , highlighting a potential role for the ddit3 node in early development of the erythroid lineage . previous reports of ddit3 function in late differentiation of erythro - leukemic cells ( coutts et  al . , 1999 ) interestingly , our inspection of cfu - e potential in mouse bm stem and erythroid progenitor compartments upon ddit3 ablation did not reveal a late erythroid defect ( c.f . , unpublished data ) . the capacity of tfs to re - program cells has been used as a test of their lineage - determining capacity . in the case of gata1 , its introduction into gmps results in the expression of erythroid potential with the appearance of large blast - like multipotential colonies in  vitro ( heyworth et  al . , 2002 ) . molecular analysis suggests that ddit3-enforced expression in gmps leads to expression of erythroid - affiliated genes and an overall transcriptional state similar to that of pre - gm , the precursors of gmps . one may presume that since pre - gm cells lie upstream of gmps their transcriptional programs should have diverged less from multipotent cells that retain both gm and e potential . this would explain the e signatures seen in ddit3-expressing gmps . by which cellular mechanism does ddit3 effect these changes ? one possibility is that rare pre - gms exist within prospectively isolated gmp populations and that these are preferentially selected for by ddit3 . this seems unlikely , given the kinetics of the changes in cells and gene expression observed . alternatively , ddit3 may regulate a subset of the erythroid program that is simply overlaid on the existing gmp program , resulting in a mixed - lineage program that is reminiscent of pre - gm cells . in such a model , ddit3 effects a pre - gm state from gmps , recapitulating its physiological role in lineage programming . recently , nerlov and colleagues have suggested that high levels of erythropoietin ( epo ) may have instructive effects on lineage specification in  vivo ( grover et  al . , 2014 ) , and interestingly ddit3 has been suggested as a target of epo signaling in erythro - leukemic cells ( coutts et  al . , 1999 ) . the strong network association observed between ddit3 and gata2 may indicate that ddit3 acts on an early erythroid signature primed by gata2 in cells of mixed - lineage potential ( may et  al . , 2013 ) and present in the pre - gm state . elements of this signature may be required for erythroid lineage progression and , thus , explain the erythroid defect observed upon ddit3 loss of function . it is interesting to speculate as to the existence of cross - antagonistic interactions between gata2-centered networks and c / ebp - driven myeloid programs putatively effected through ddit3 . ddit3 heterodimerizes with c / ebp family members to form complexes that can not bind dna , thus blocking activation of c / ebp - driven programs , a mechanism that has been described to block differentiation in mesenchymal lineages ( han et  al . , 2013 ; shirakawa et  al . , 2006 ) . a combination of ddit3 structure - function mutant studies and direct investigation of gata2 and c / ebp dna binding in ddit3-expressing gmps will contribute to clarifying dislodgement of c / ebp complexes from their target genes as a putative mechanism of gmp lineage remodeling . b6.129s - ddit3tm1dron / j ( ddit3 ko ) mice ( jackson immunoresearch laboratories ) and c57bl/6 mice were maintained in the john radcliffe hospital and cr - uk london research institute animal facilities in accordance with home office regulations . fdcpmix culture conditions , lentiviral transductions with gata1-ert and pu.1-ert constructs , and tamoxifen activation were performed as described previously ( may et  al . , 2013 ) . lentiviral transductions of fdcpmix cells with ddit3-small hairpin rna ( shrna ) were performed under maintenance culture conditions ; transductions of primary bm cells were performed in serum - free expansion medium ( sfem ) ( stemcell technologies ) supplemented with mouse stem cell factor ( scf ) and flt3l ( 50  ng / ml ) . , gfp gmps were cultured for up to 5  days in iscove s modified dulbecco s medium ( imdm )  + 10% fetal calf  serum ( fcs ) supplemented with mouse scf ( 50  ng / ml ) , mouse il-3 and il-6 ( 10  ng / ml ) , and human recombinant erythropoietin ( eprex , 10  u / ml ) . colony - forming cell ( cfc ) assays used m3234 supplemented with rat scf ( 100  ng / ml ) ; mouse il-3 ( 0.01  ng / ml ) and eprex ( 10  u / ml ) ( fdcpmix ) ; and m3434 or m3234 supplemented with mouse scf ( 100  ng / ml ) , mouse il-3 , il-11 , gm - csf , and tpo ( 10  ng / ml ) , and eprex ( 10  u / ml ) . all mouse cytokines were from peprotech and cfc media were from stemcell technologies . transcriptional profiling of up to 48 genes in individual cells was performed on  a fluidigm platform and the data retrieved and quality - controlled as described previously ( teles et  al . , 2014 ) . the ct values were calculated to  the mean of the three control genes utilized . heatmap representation of z  score - normalized ct values was performed in genesis ; pca plots used the statistical toolbox in matlab ( mathworks ) . single - cell gene expression data were linearly transformed as described previously ( teles et  al . , 2013 ) . logistic regression linear classifiers were used to infer the best predictor genes in the separation between two cell populations ( teles et  al . , 2013 ) . single - cell transcriptional networks were inferred by calculating significant pairwise associations using both continuous ( spearman rank correlations ) and binary ( or ) components of linearly transformed expression data . spearman rank correlations were calculated between all pairs of genes co - expressed by a minimum of ten cells in a given population . correlation coefficients > 0.4 with p  < 0.01 were considered to be significant . or and respective 95% confidence intervals ( cis ) were calculated based on presence / absence patterns of expression for all pairs of genes in a given population . significant positive and negative associations were called when lower95ci > 1 and upper95ci  < 1 , respectively . network representations of significant pairwise associations in both methods were produced using cytoscape ( smoot et  al . t.e . supervised the project . c. pina and t.e . wrote the paper with input from j.t .