Lineage tracing on transcriptional landscapes links state to fate during differentiation

A simultaneous assay of clonal states and fates

Our strategy for simultaneously capturing transcriptional cell state and fate is to genetically barcode a heterogenous progenitor population, allow cell division, sample some cells immediately for scSeq profiling, and the remainder later (13). This approach provides data for three types of clonal relationships (Fig. 1a): (1) sister cells in the earliest time point may be captured after 1 or 2 rounds of division; (2) clones observed at both early and later time points allow comparing the state of an early cell to the fate outcomes of its sisters; (3) sampling differentiated cells at later time points will reveal clonal relationships between different fates. If recently-divided sister cells (type 1) are transcriptionally similar, then pairs of clonally-related cells sampled both early and late (type 2) should reveal how single cell gene expression changes over time during differentiation. This approach can map the fate of cells from a continuous landscape of starting states and does not require isolation or labeling of specific prospective progenitor populations (2, 14).

We modified a classical strategy for clonal labeling by lentiviral delivery of inherited DNA barcodes (15, 16), to allow barcode detection using scSeq (17). The barcode consists of a random 28-mer in the 3’ UTR of an enhanced green fluorescent protein (eGFP) transgene under control of a ubiquitous EF1α promoter (Fig 1b). Transcripts of eGFP are captured during scSeq, and the barcode is revealed through analysis of sequencing reads. We generated a library of ~0.5 × 10⁶ barcodes, sufficient to label 5000 cells in an experiment with <1% barcode overlap between clones (see Supp Methods 2.3 for estimate of diversity). We refer to the barcoding construct as LARRY (Lineage And RNA RecoverY).

We tested LARRY on mouse embryonic stem (ES) cells and primary HPCs. After profiling by scSeq, one or more barcodes could be robustly detected in 93% of GFP+ cells (Supp Fig 1a–c). Specific barcode sequences overlapped rarely between replicate transduction experiments, at a frequency expected by chance for the library size (0.3% of 5000 barcodes appeared more than once). Therefore, the approach provides an efficient method for simultaneously barcoding large numbers of cells for combined fate and state mapping.

To analyze HSPC fate potential, we applied LARRY to cells cultured in vitro, and to cells transplanted in vivo. For in vitro analysis, we isolated a broad class of oligo-potent (Lin-Sca-Kit+) and multipotent progenitors (Lin-Sca1+Kit+ or LSK) cells (Supp Fig 2a,b) and plated them in media chosen to support broad multi-lineage differentiation (see Methods). Following barcode transduction, cells were cultured for two days to allow lentiviral integration and subsequent division. During this time the cells divided three times on average. We then sampled half the cells (defining the ‘early state’) for scSeq. The other half were re-plated and then sampled after two days (30% of cells) and four days (remaining cells) (Fig 1c). For transplantation, Lin-Sca(hi)Kit+ cells, consisting of mostly short-term and long-term hematopoietic stem cells (HSCs) (Supp Fig 2a,b), were barcoded and placed in culture. After two days, 40% were profiled by scSeq, with the remainder transplanted into ten sublethally irradiated host mice (10) and recovered for scSeq one and two weeks later (Fig 1h). We retrieved 130,887 scSeq transcriptomes from culture and 182,173 single cells after transplantation (see Table 1 and Supp Methods 3 for details of analysis). In these two experiments, 38% and 63% of cells respectively belonged to a clone of two or more cells (5864 and 7751 clones), with 1,816 and 817 clones in total spanning early and late timepoints (Fig 1d,i).

We visualized the cell transcriptomes using force-directed layouts (SPRING plots (18)). In vitro, the cells defined a continuous state map spanning from multipotent progenitors (MPPs) to nine mature cell types that appeared in culture (Fig 1e,f): erythrocytes (Er), megakaryocytes (Mk), basophils (Ba), mast cells (Ma), eosinophils (Eos), neutrophils (Neu), monocytes (Mo), dendritic cells (plasmocytoid pDC; Ccr7+ migratory migDC) and lymphoid precursors (Ly). On this landscape clones exhibited a range of behaviors including uni-lineage and multi-lineage differentiation, and self-renewal of early progenitors (Fig 1g). After transplantation, the cells again defined a continuous landscape spanning from MPPs through several stages of neutrophil maturation, as well as DCs, Mo, Er, B, T and Ba cells. Many of these cell types were internally heterogeneous, with several types of DCs including CD11+, CD8+, migDC and pDC, as well as Ly6C+ classical and Ly6C- non-classical monocytes (Fig 1j; Supp Fig 3). We did not detect megakaryocytes, possibly because they did not survive bone marrow harvest, flow-sorting and single-cell encapsulation intact. Therefore, with these experiments we simultaneously captured single cell state maps and their underlying clonal relationships.

Clonal dynamics identify early transcriptional fate boundaries

With LARRY it is possible to estimate how a single cell changes over time by sampling a clone across multiple time points. Yet the accuracy of this approximation depends critically on the similarity of sister cells at the earliest time point. We found that pairs of sisters profiled on day 2 localized in the SPRING graph, had correlated gene expression (median R=0.846) and that a majority (70%) fell in the same or nearest neighbor cluster (Fig 2a,b; Supp Fig 4a–d; Supp Methods section 5). A minority of cells, however, were more diverged, with 10% falling outside a four-cluster radius (compared to 80% for random cell pairs). We tested and ruled out that similar sister pairs are technical co-encapsulation artifacts (Supp Fig 4e). These tests justified approximating single-cell trajectories by clonal trajectories, though with some loss in resolution of fate boundaries expected due to ~10% diverged sister pairs.

Beginning with the in vitro data, we recorded the clonal fates of each day 2 cell. Visualizing cells from uni-lineage clones revealed well-delineated domains of fate potential (Fig 2c). Where the progenitors for different fates overlapped, we observed bi-potent or oligopotent clones, indicating the location of fate commitment boundaries (Fig 2d). The true number of multipotent clones is likely under-estimated in our data, since some clonal fates were likely missed due to under-sampling (Supp Fig 5) and cell commitment prior to division would result in only one observed fate. Consistent with recent scSeq studies (19), progenitors with different fate potentials did not partition into discrete cell states, but instead formed a structured continuum. Further, bipotent domains formed extended fate boundaries, indicating that differentiation progression can occur independently of fate commitment over some time. Both of these observations differ from the classical model of hematopoiesis represented by discrete, step-wise transitions in state and fate potential.

We interrogated the gene expression heterogeneity defining this continuum and its fate potential. The multipotent progenitor (CD34+) fraction of day 2 cells (Fig 2e) contained several broad domains, including a restricted central domain of stem cell marker (Procr) expression, a wing expressing Gata2 – an erythroid and stem cell marker – and an opposing wing expressing Flt3 indicative of lymphoid priming. Overlaying clonal outcomes (Fig 2f) revealed regions of functional lineage priming consistent with these broad expression domains, but further segregated into subdomains. Megakaryocyte, basophil, mast cell and eosinophil potential were all restricted to the Gata2+ region, yet derived from separate subsets within this region. Testing for differential gene expression, we identified genes enriched within each subdomain of fate potential (Fig 2g) revealing known markers and many that have not been characterized in hematopoiesis (n=447 [391 unique] differentially expressed genes at FDR=0.05; Table 3). For example, Ikaros family zinc finger 2 (Ikzf2) – a myeloid leukemia gene not previously associated with fate choice – was enriched in eosinophil and mast cell progenitors, but not basophil or megakaryocyte.

We similarly identified gene expression correlated with fate outcomes in less differentiated ST-HSCs and LT-HSCs transplanted into irradiated mice. As before the cells spanned a continuous landscape with domains of primed gene expression, including a central domain of stem cell (Procr) and opposing wings of Gata2 and Flt3 expression (Fig 2h), that correlated with output into the nine respective post-transplant fates (Fig 2i). Despite the less mature state of these cells, each fate outcome correlated with unique enriched genes prior to transplantation (Fig 2j; n=190 [173 unique] differentially expressed genes at FDR=0.05; Table 3), indicating a surprising degree of specific priming at this early stage of differentiation. The differentially expressed genes represented a wide range of functional gene categories, from cell adhesion to chromatin regulation as well as intra- and extra-cellular signaling, with cytokine signaling as the major enriched category (p < 10⁻⁵; Table 4). Gene set enrichment analysis for each fate revealed terms associated with the fate’s function, such as ‘lymphocyte activation’ (𝑝𝑝 = 0.002) for T cell progenitors and ‘response to bacterium’ (𝑝𝑝 = 0.001) for neutrophil progenitors. Surprisingly, the majority of top terms enriched in erythrocyte progenitors related to cell motility (8 out of the top 10 terms; Table 5), possibly indicating these progenitors are primed to undergo cytoskeletal and niche rearrangements. We observed differences in clonal fate of phenotypically-similar progenitors (day 2) in vivo compared to in vitro (Supp Fig 6). Such environmental plasticity acts at sub-clonal resolution, as seen by barcoding HSPCs and culturing them in different cytokines (n=958 clones sampled between conditions; n=1,600 clones across time points within conditions; Supp Fig 7a–d). When split across cytokine conditions, sister cells showed consistent shifts of clone size and observed cell fate (Supp Fig 7e–g).

Overall, these observations support the view that functional lineage priming varies across a continuous hematopoietic progenitor landscape and covaries with the heterogenous expression of genes, including transcription factors and a wide array of other functional gene categories. The observed clonal outcomes reflect both priming and environmental inputs.

How predictable is cell fate from gene expression?

Several factors impinge on the fate choice of a cell, including interactions with the environment, gene expression, chromatin state, and stochastic molecular events. Single-cell RNA-seq provides only a limited view of cell state. So far, we have considered the correlates of future fate choice revealed by this measurement. We now ask: to what extent can fate be predicted from scSeq data?

To estimate the predictability of fate choice from gene expression, we considered the machine learning task of predicting a cell’s dominant fate outcome (Fig 3a,b) based on its present scSeq profile (Supp Methods 9.1). We used two machine learning methods: logistic regression a neural network [multi-layer perceptron]. We applied these methods to several sets of genes, including all highly variable genes, genes that are differentially expressed between progenitors (Table 3), and a genome-wide set of transcription factors (n=1811). Transcription factors were only marginally more informative than random size-matched gene sets (10% more informative in vitro; 3% more informative in vivo), whereas differentially expressed genes were substantially more informative (38% more informative in vitro; 20% more informative in vivo). Augmenting the differentially expressed genes with all highly variables genes, which increased the number of genes used by 12-fold in vitro and 28-fold in vivo, did not lead to a significant change in accuracy (1% change in vivo, −4% change in vitro). These results suggest that the predictive content of our gene expression measurements in HSPCs is almost entirely contained within several hundred differentially expressed genes, and only marginally enriched in transcription factors. The poor performance of transcription factors may be due to their low and noisy expression levels, or to the comparable influence of other functional gene categories. These results were recapitulated when predicting the full distribution of fate outcomes rather the dominant one (Supp Fig 8g–j). Viewing predictive accuracy at the single-cell level revealed greater accuracy for increasingly mature cells (Supp Fig 8k–n; Supp Methods 9.2). Across all conditions, the highest overall predictive accuracy from transcriptional state was 60% in vitro and 51% in vivo. These figures provide a lower bound for the cell-autonomous influence of transcriptional state on cell fate.

Functional purity of scSeq-defined cell states

Although fate prediction accuracy could be limited by stochastic fluctuations in cells or their environment, it is also possible that stable cellular properties influence fate choice but are not detected by scSeq. If such ‘hidden variables’ (4) exist, they would challenge the view that scSeq can define functionally pure populations. We tested for the presence of hidden variables by comparing ‘early’ and ‘late’ modes of cell fate prediction. If there were no hidden variables, we reasoned that the information shared between separated sister cells could only decrease as time passes. Conversely, if there are stable properties that influence cell fate but are hidden from scSeq, then the mutual information between sisters could increase over time as these properties manifest in cell fate. This reasoning reflects a formal result known as the Data Processing Inequality (20) (Supp Methods 10.1).

To compare the accuracy of early vs late prediction, we applied a panel of machine learning algorithms to guess the dominant fate of a clone using either the transcriptomes of its day 2 sisters (as in Fig 3a,b), or the transcriptomes of its sisters separated four days in culture (n=502 clones) or one week post-transplantation (n=69 clones) (Fig. 3c–h). We found that late prediction was more informative for all algorithms tested (Fig 3e,h) with the most accurate algorithms achieving late prediction accuracy of 76% in vitro and 70% in vivo, compared to 60% and 52% respectively for early prediction.

These improvements in accuracy for late prediction reflect the high rate of concordance between sisters cell fates, and hold for clones of all potencies (Fig 3d), consistent with recent observations of clonal fate restriction among HSPCs (10). Clones in separate wells produced identical combinations of fates 70% of the time, compared to 22% by chance. One week post-transplantation, sister cells in separate mice also showed highly concordant fate outcomes (Fig 3g): although they only shared the exact same combination of fates 29% of the time (compared to 10% by chance) they shared the same dominant fate 71% of the time (23% by chance). Together, these results imply that, both in culture and during transplantation, there are heritable properties of cell physiology that influence cell fate but are not evident in our scSeq measurements. We cannot tell whether information on cell fate is restricted simply because scSeq data is noisy, or because cell fate depends on cellular properties that are not reflected in the transcriptome, such as chromatin state, protein abundances, cell organization, or the microenvironment.

If scSeq states are not functionally pure then phenotypically-similar progenitors should be primed towards different fates. We tested this prediction by analyzing clones that were detected in three separate samples from our in vitro dataset: at day 2, and in two wells separated until day 6 (n=408 clones; Fig 3i). Without hidden variables, the two fates observed at day 6 should be statistically independent after conditioning on the day 2 state. In this case, the expected frequency of different fate outcome in the separate wells (‘mixed clones’) can be calculated (Fig 3i, left; Supp Methods 10.4). As a result of fate priming, however, we predicted that the frequency of mixed clones rooted in phenotypically-similar day 2 cells would fall below this expectation. For each of three fate choices (Neu vs. Mo, [Neu/Mo] vs. [Er/Mk/Ma/Ba], and [Ly/DC] vs. [all myeloid]), and across different day 2 progenitor states, the proportion of mixed clones was significantly below the expectation for pure bi-potency (Fig 3j–l; Supp Fig 9a,b). This analysis supports the previous conclusion that cell-autonomous fate biases can indeed coexist in the same measured scSeq state.

The above evidence for hidden variables suggests limits to the use of scSeq in building atlases that resolve the functional complexity of HSPCs. For years, cytometry (FACS) has been used to dissect the hematopoietic hierarchy with increasing precision, with the ultimate goal of defining functionally pure subsets of progenitors. Recent studies showing that many commonly used FACS gates are heterogenous in fate and transcriptional state have raised the possibility that genome-wide assays such as scSeq might be required to achieve the necessary resolution. These results indicate that scSeq, while informative, may still be insufficient for defining functionally pure progenitor states.

Distinct routes of monocyte differentiation

Clonal analysis can reveal differentiation paths that may not be apparent by scSeq alone. In the data, monocytes appeared to form a spectrum from neutrophil-like to DC-like, expressing alternatively neutrophil elastase (Elane) and other neutrophil markers, or MHC class II components (Cd74 and H2-Aa) (Fig 4a). No similar overlap occurs with other cell types (Supp Fig 10a). We investigated whether this phenotypic spectrum might result from distinct differentiation trajectories of monocytes (21).

To determine monocyte ontogenies, we scored their clonal relatedness with mature neutrophils and DCs. The monocytes were not uniformly coupled to either cell type (Fig 4b): those with increased expression of neutrophilic markers were clonally related to neutrophils (Supp Fig 10b; 𝑝𝑝 < 10⁻⁷, Mann-Whitney U test), whereas those with DC-like gene expression were clonally related to DCs and lymphoid cells (𝑝𝑝 < 10⁻¹⁷). We did not observe a comparable phenomenon for any other cell type in our data. Thus, monocytes appear unique in showing a phenotypic spectrum that correlated with distinct clonal histories.

The distinct clonal origins of monocytes suggested that they arise from progenitors with different fate potentials, and possibly different gene expression. To define their progenitors, we classified the differentiating monocytes (4–6 days) as either DC-like or neutrophil-like (Supp Methods 11.1), and then examined their early sisters (2 days). Indeed, the predecessors of DC-like and neutrophil-like monocytes segregated by gene expression (Fig 4c,d), with respective expression of early DC and lymphoid markers (Flt3, Bcl11a and Cd74) or early neutrophil markers (Elane, Mpo and Gfi1; see Table 6 for a full list of differentially expressed genes). These early differences were mostly distinct from those distinguishing mature (4–6 day) DC-like and neutrophil-like monocytes (Fig 4e; Table 7). Our data therefore contains two different pathways of monocyte differentiation with distinct clonal relationships and gene expression dynamics.

These results are consistent with a recent finding that immunophenotypically-defined monocyte-dendritic progenitors (MDPs) and granulocyte-monocyte progenitors (GMPs) give rise to monocytes with DC-like and neutrophil-like characteristics, respectively (21). To test whether our observations represent MDP/GMP outputs, we performed scSeq on fresh MDPs and GMPs sorted from adult mouse bone marrow and found that they co-localized with the day 2 progenitors of DC-like and neutrophil-like monocytes. Similarly, scSeq analysis of MDPs and GMPs cultured for 4 days in vitro co-localized with mature DC-like and neutrophil-like monocytes (Supp Fig 10c). Thus, the DC-like and neutrophil-like trajectories observed here likely represent MDP and GMP pathways of monocyte differentiation, and they clarify the location of these states in a gene expression continuum.

Several lines of evidence support the existence of distinct monocyte-neutrophil and monocyte-DC clonal couplings in vivo, and not only in culture: (i) clonal and gene expression relationships following transplantation; (ii) persistent heterogeneity in freshly-isolated mouse and human monocytes; and (iii) results from non-perturbative in-vivo clonal analysis. We present these results in turn.

First, one week after transplantation, monocytes showed distinct clonal relationships to neutrophils and DCs (Fig 4f). As in vitro, the DC-related monocytes were enriched for DC marker genes, whereas neutrophil-related monocytes were enriched for neutrophil markers (Fig 4g,h). Second, we analyzed classical monocytes (Supp Fig 10g) and human peripheral blood monocytes (Supp Fig 10h) by scSeq. Principal component analysis shows that in both cases a spectrum exists of neutrophil-like to DC-like gene expression (see Table 8 for differentially expressed genes), which is also evident in the expression of marker genes (Fig 4i). This analysis agrees with earlier observations (21). Third, in native hematopoiesis we examined the clonal co-occurrence of monocytes with DCs and neutrophils after genetically barcoding HPC clones in a non-perturbative manner using a transposase-based strategy (22) (Supp Methods 12.3). If monocyte heterogeneity correlates with distinct clonal coupling to neutrophils versus DCs, we would expect an anti-correlation between neutrophil and DC relatedness among monocytes (Fig 4j). After a 12 week chase (Fig 4k), we indeed found significantly fewer neutrophil-monocyte-DC tags than would be expected if clonal co-occurrence were independent (2.5-fold reduction; p < 0.001 by binomial test of proportion; Fig 4l,m). Overall, our results support the existence of multiple monocyte ontogenies in native hematopoiesis as well as in culture and during transplantation.

A benchmark for fate prediction in hematopoiesis

To understand hematopoietic fate control, we and others have been interested in developing data-driven models of gene expression dynamics constrained by scSeq data (3, 5, 7, 4, 23). Computational models could identify cellular components driving fate choice, and the sequence of gene expression changes that accompany cell maturation. Due to a lack of ground truth data, existing methods have been hard to compare and validate. Here, we asked whether common approaches for modeling cell state dynamics are consistent with our clonal tracking data.

scSeq-based models do not fully predict fate choice

We first asked how well existing computational models, using only scSeq data, predict cell fate probabilities. We tested three recent approaches, Population Balance Analysis (PBA) (4), WaddingtonOT (WOT) (5), and FateID (7) for their ability to predict the fate of a cell choosing between neutrophil and monocyte fates in culture. We calculated for each cell at day 2 the fraction of its clonal relatives that became a neutrophil or a monocyte (Fig 5a), and then attempted to predict this fraction from transcriptomes alone (Fig 5b; Supp Methods 13.2–4). All three methods were broadly consistent with clonal fate bias as cells began to mature, but in the early progenitor (Cd34+) region, clonal tracking revealed a bifurcation of monocyte and neutrophil potential that was generally not detected by the prediction algorithms, although FateID performed slightly better (Fig 5c,d; R < 0.26 for all methods). All fell considerably below fate predictions obtained from held-out clonal data (R = 0.5; Supp Methods 13.5; correlation is low overall because of noise in the fate outcomes of single cells). These results show that in the absence of lineage information, computational methods may mis-identify fate decision boundaries. It is therefore significant that when genes are ranked by their ability to predict cell fate bias, the top ten genes easily outperformed the prediction algorithms (Fig 5d), including known fate regulators such as Gata2 and Mef2c (Fig 5d,e). The selection of the correct genes to use for prediction, however, required clonal information. These results provide a framework for comparing computational models of differentiation, and may serve as a useful benchmark for improving them.

Temporal progression is captured by pseudotime

A common goal of scSeq is to order gene expression along dynamic trajectories by defining a ‘pseudotime’ coordinate that orders transcriptomes (24). At present, it is unknown how single cells traverse these trajectories, including whether they progress at different rates or even reverse their dynamics (4). Focusing on neutrophil differentiation as a test case, we asked how well ‘pseudotime’ describes the kinetics of differentiation as revealed by clonal tracking. We ordered cells from MPPs to GMPs, to promyelocytes (PMy), to myelocytes (My) (n=63,149 cells Fig. 5f; Supp Fig 11a; Supp Methods 14.2), and compared the pseudotemporal progression of clones sampled at two consecutive days (Fig. 5g). This analysis showed a strikingly consistent forward velocity along differentiation pseudotime. By integrating the velocity across the trajectory, we were able to calculate pseudotime progression as a function of real time for a typical cell (Fig. 5h; Supp Methods 14.3). The time for an MPP to differentiate into a myelocyte was 10 days, consistent with prior literature (25). Pseudotime analysis of sister cells differentiated in separate wells also showed a consistent pace of differentiation both shortly after cell division (day 2), and remaining so four days later (R≥0.89; Fig 5i). Pseudotime velocity was most variable among MPPs (Figs. 5j), which could be explained by cells remaining in the MPP state for a variable duration before initiating neutrophil differentiation. These results support the use of pseudotime methods for mapping differentiation progression.

Agreement of state and clonal differentiation hierarchies

For cells undergoing multi-lineage fate choice, scSeq has been used to estimate lineage hierarchy based on the assumption that cell types with transcriptionally similar differentiation pathways are clonally related (3–5, 7). Yet this assumption may not always hold: similar end states could also arise from non-overlapping clones (26) and distant end states could share lineage through asymmetric division.

To compare fate hierarchies constructed using lineage and state information, for each pair of differentiated states we quantified the number of shared clones as well as the similarity of cell states for each pair of differentiated fates both in vivo and in culture (Fig. 5k,l,p,q; Supp Methods 15.1–2). We found that measures of state distance and clonal coupling are closely correlated in vitro (r=0.93, p < 10⁻³⁵; Fig 5m). When we constructed candidate cell type hierarchies from state distance and clonal distance respectively (Fig 5n,o), they were almost identical, with only one difference in the differentiation path assigned to mast cells. These results held for a broad range of parameters and for different distance metrics (Supp Fig. 12a–h). In vivo, however, the same analysis revealed a weaker correlation between state and fate distance (r=0.58; p = 0.065; Fig 5r) with considerable differences between the resulting cell type hierarchies (Fig. 5s,t). Several factors might explain the weaker relationship between state and fate hierarchy in vivo, such as the longer interval between samples (1 week, compared to every 2 days in vitro), or the complex differentiation environment. These results suggest a set of experimental parameters – operant in our in vitro experiment – that may be favorable for inferring clonal relationships from gene expression topology: dense sampling over time, uniformity of the differentiation environment, and a spectrum of the maturity in the initially barcoded cells.