Multiplexed Division Tracking Dyes for Proliferation-Based Clonal Lineage Tracing

Determining the contribution of asymmetric cell division (ACD), intercellular communication, quorum sensing, lineage priming, and autonomous programming to clonal cell fate is a key focus of immunology and many other fields of biology (1–6). However, progress has been impeded by the low throughput and laborious nature of common lineage-tracing and fate-mapping approaches such as time lapse microscopy. Recently introduced technologies such as retroviral barcoding (4, 7), CRISPR-induced heritable genetic lesions (8), and the development of fluorescent lineage reporters (2, 9, 10) have contributed to improved throughput in lineage-tracing experiments and thus revealed important discoveries into the emergence of diverse cell types across multiple systems. Despite this success, such methods remain highly resource dependent and time consuming. Furthermore, these methods typically lack information regarding clonal division progression, an important source of information in understanding the mechanisms that drive cell fate decisions (11–17). Although many cellular processes across multiple systems have demonstrated an association between cell state transitions and division (12, 15, 17–20), thorough examination of these associations across the progeny of expanding single-cell lineages has, to date, been limited. A fast, easy, high-throughput method that, for individual clones, simultaneously measured division progression as well as cell state in the form of marker and/or fluorescent reporter expression could therefore significantly contribute to progress in this field. In this study, we introduce such a method, using multiplexed division tracking dyes in combination with flow cytometry–based phenotyping. Earlier variants of the dye-multiplexing approach have been applied to high-throughput cytotoxicity assays (21), analysis of clonal division progression (22), and identification of distinct cocultured cell populations (23). In this study, we demonstrate that the utility of this method can be significantly extended by integrating phenotypic information with proliferation-based lineage tracing and by provision of the statistical tools necessary for data interrogation.

The three murine strains, wild-type C57BL/6, OVA-specific Bim-deficient (OT-I/Bcl2l11^−/−) (24), and Blimp-1^gfp/+ (25) mice were maintained under specific pathogen-free conditions in the Walter and Eliza Hall Institute (WEHI) animal facilities and were used at 6–10 wk of age. All experiments were performed under the approval of the WEHI Animal Ethics Committee.

Spleens and lymph nodes were homogenized through a 70 μM cell strainer to generate single-cell suspensions. CD8⁺ T cells were isolated by negative selection using EasySep Mouse CD8⁺ T cell Isolation Kits according to the manufacturer’s protocol (STEMCELL Technologies).

Purified CD8⁺ T cells were resuspended in sterile PBS containing 0.1% BSA (PBS 0.1% BSA) and labeled with either 5, 2.5, or 0 μM CFSE (Invitrogen) at a density of ≤2 × 10⁷ cells/ml at 37°C for 10 min and washed twice with 10 ml of ice-cold RPMI 1640 containing 10% FCS.

Cells were resuspended in PBS 0.1% BSA, and those that were labeled with 5 μM CFSE were further labeled with either 5, 2.5, or 0 μM CellTrace Violet (CTV) (Invitrogen). Cells labeled with 2.5 μM CFSE were labeled with 5 μM CTV. Cells labeled with 0 μM CFSE were labeled with either 5 or 0 μM CTV. All labeling was performed at a density of ≤2 × 10⁷ cells/ml at 37°C for 20 min, and all cells were washed twice with 10 ml of ice-cold RPMI 1640 10% FCS.

Cells were resuspended in PBS 0.1% BSA and labeled with either 5 or 0 μM Cell Proliferation Dye eFluor 670 (CPD) (eBioscience) at a density of ≤2 × 10⁷ cells/ml at 37°C for 10 min and washed once with 10 ml of ice-cold RPMI 1640 10% FCS and once with tissue culture medium.

Purified CD8⁺ T cells were resuspended in PBS 0.1% BSA and labeled with either 10 or 0 μM CellTrace Yellow (CTY) (Invitrogen) at a density of ≤2 × 10⁷ cells/ml at 37°C for 20 min and washed twice with 10 ml of ice-cold RPMI 1640 10% FCS.

Cells were resuspended in PBS 0.1% BSA and labeled with either 5 or 0 μM CTV at a density of ≤2 × 10⁷ cells/ml at 37°C for 20 min and washed twice with 10 ml of ice-cold RPMI 1640 10% FCS.

Cells were resuspended in PBS 0.1% BSA and labeled with either 5 or 0 μM CPD at a density of ≤2 × 10⁷ cells/ml at 37°C for 10 min and washed once with 10 ml of ice-cold RPMI 1640 10% FCS and once with tissue culture medium.

Tissue culture medium was RPMI 1640 with 10% FCS, 1 mM sodium-pyruvate, 2 mM GlutaMAX, 10 mM HEPES, 100 U/ml penicillin, 100 μg/ml streptomycin (all Invitrogen), and 50 μM 2-ME (Sigma-Aldrich). Purified C57BL/6 CD8⁺ T cells were stimulated with 10 μg/ml plate-bound anti-CD3 mAb in flat-bottom 96-well plates (WEHI facility, clone 145-2C11). For some experiments, OT-I/Bcl2l11^−/− CD8⁺ T cells were stimulated with 0.01 μg/ml SIINFEKL (N4) peptide (Auspep). The use of Bim-deficient T cells enhances survival in vitro but does not affect proliferative or phenotypic behaviors (24, 26). These cells were stimulated in round-bottom 96-well plates at a density of 20,000 cells per well. This protocol leads to the self-presentation of peptide by T cells and is used as a minimal culture system to enable the addition of further costimulatory signals (22, 24, 27, 28). Cells were cultured in 200 μl of tissue culture medium in the presence of 25 μg/ml anti-mouse IL-2 mAb (WEHI mAb facility, clone S4B6) that inhibits the activity of mouse IL-2 but does not act on recombinant human IL-2 (rhIL-2) (29). rhIL-2 (PeproTech), anti-CD28 (WEHI mAb facility, clone 37.51), mouse IL-4 (purified from baculovirus-transfected Sf21 insect cells), and mouse IL-12 (130-096-707; Miltenyi Biotec) were added to cultures where indicated. Cells were incubated at 37°C in 5% CO₂.

Purified C57BL/6 and OT-I/Bcl2l11^−/− CD8⁺ T cells were sequentially labeled with CFSE, CTV, and CPD. The uniquely labeled cell populations were mixed (except for the unlabeled and CPD-only labeled controls) and stimulated under conditions indicated. After 22–26 h, prior to the first division, cells from across multiple wells stimulated under the same conditions were pooled and sorted according to their distinct fluorescence profiles into new wells such that each well contained a single cell from each unique labeling profile, with the exception of cells labeled with only 5 μM CPD or unlabeled cells. Wells contained the same conditions under which the cells were initially stimulated. Bulk population controls were also sorted into new wells, with 1000 cells from each labeling profile sorted into separate wells, in addition to 100 cells from each population sorted into the same well. This gave bulk populations of each fluorescence profile both separately and in combination. This included cells labeled with 5 μM CPD only and unlabeled cells.

Purified Blimp-1^gfp/+ CD8⁺ T cells were sequentially labeled with CTY, CTV, and CPD and stimulated with plate-bound anti-CD3 (10 μg/ml), rhIL-2 (31.6 U/ml), and mouse IL-12 (10 ng/ml) in the presence of S4B6 (25 μg/ml) and were subsequently cultured and sorted according to the same criteria as above. Sorting was performed on either a BD Biosciences FACSAria III or a BD Biosciences Influx cell sorter.

At time points indicated, cells were stained on ice with indicated Abs used at the following concentrations: 1:400 dilution anti–CD8-APC-Cy7 (clone 53-6.7; BD Biosciences), 1:1600 dilution anti–CD62L-PE (clone MEL-14; BD Biosciences), 1:1600 dilution anti–CD62L-APC-Cy7 (clone MEL-14; BD Biosciences), 1:400 dilution anti–CD25-PE-Cy7 (clone PC61; BD Biosciences), and 1:3200 dilution anti–CXCR3-PE-Cy7 (eBioscience). A total of 10⁴ beads (Rainbow Calibration Particles BD Biosciences) and 0.2 μg/ml propidium iodide (PI; Sigma-Aldrich) was also added to samples prior to analysis. An Ab-staining mixture containing all relevant Abs along with beads and PI was prepared for each experiment. Ab-staining mixture was added at staggered time points (∼2–3 min apart) to each sample in the 96-well culture plates and later transferred to 5-ml polystyrene tubes such that each sample was stained for as close to 30 min as possible prior to immediate acquisition of as much of the sample as possible (duration of acquisition lasted ∼2–3 min per sample). Analysis was performed on a BD Biosciences LSRFortessa X20. Gates were set using labeled bulk population controls, and these were then applied to clonal data. Live lymphocytes were identified using forward and side scatter and PI exclusion. Cells were separated into CPD⁺ and CPD⁻ populations, and division gates were identified for each labeling profile on CFSE versus CTV for C57BL/6 and OT-I/Bcl2l11^−/− or CTY versus CTV for Blimp-1^gfp/+. Clonal families were identified, and the division numbers and expression levels of surface markers or Blimp-1^gfp/+ of each cell was enumerated and exported for data visualization and further analysis.

We first give a general outline of the permutation testing procedure (30) and then provide specific detail for the data and hypotheses described in the main text. Given a data set of n-ordered observations D = (Z₁, Z_2,,…,Z_n), a permutation π of them is a reassignment of the labels of the individual datum, i → π(i), to create the reordered data set D_π = (Z_π(1), Z_π(2),…,Z_π(n)). For example, if π(i) = n + 1 − i for all i, then D_π = (Z_n, Z_n−1,…, Z₁) is the original data but in reverse order. The principle of permutation testing is to evaluate a statistic on the recorded data, T(D), in which the statistic depends on the data order. Under the null hypothesis, H₀, that a certain set of orders of the data were equally likely and characterized by a collection of permutations, π in a set Q, T(D) can be compared with the distribution of the statistic computed on the reordered data sets {T(D_π)}_π∈Q. In particular, denoting by |A| the number of elements in a set A, the proportion of permutations that lead to a statistic that is lower than that observed for the true data order is the lower p value,

whereas the proportion of permutations that lead to a statistic that is higher than that observed for the true data order is the upper p value

To realize a permutation test successfully, it is important that the collection of allowed permutations accurately describe the null hypothesis and that the test statistic tends to deviate from the true statistic if the null hypothesis is not true.

For many tests, the number of possible permutations |Q| is too large for $T(Dπ)$ to be computed for every permutation $π$⁠. For example, for data with $n$ interchangeable elements under a null hypothesis, there are $n$ factorial permutations, which grow faster than exponentially in $n$⁠. Thus, it is common to use Monte Carlo methods to estimate $pl$ and $pu$⁠. This is achieved by drawing a large number, $B$⁠, of samples from $Q$ uniformly at random and then making empirical estimates of the p values as

and

The data from the multiplex consist of a list of environments, an assignment of clones to those environments, and the generations and fluorescence levels of cells in those clones. Depending on the hypothesis to be tested, we compute a statistic on these data and compare it with the distribution of the same statistic for a collection of permutations. For Fig. 3C, we test the null hypothesis $H0$ that every cell’s fluorescence is equal in distribution irrespective of generation, clone, or environment. Our test statistic is the per-clone variance in fluorescence averaged across all clones. The set of allowable permutations is all possible reordering of cell labels, resulting in cells being reassigned among clones and environments. If in total there are $n$ cells, then there are $n$ factorial allowable permutations, and so Monte Carlo methods are needed to compute the one-sided test p values as described above.

In Fig. 3D, we test the null hypothesis that each clone’s expansion, recovery, and fluorescence levels are equal in distribution across environments. We define the statistic to be the per-environment variance in fluorescence averaged across environments. In this study, permutations are all possible relabelings of the environment label of clones, effectively swapping whole clones across environments.

In Supplemental Fig. 3A, we test the null hypothesis that, regardless of the environment in which they are found, each clone’s fluorescence levels are equal in distribution for clones at the same developmental stage (i.e., for clones that have the same number of cells in each generation). As in the previous test, the statistic is the average per-environment variance in fluorescence. What has changed is that not all permutations of clones are allowable. Instead, we first identify all families in which all cells were measured, which is achieved by testing if the sum across all cells in a clone divided by 2 to the power of the generation of each cell equals one. Of those, clones that have the same generation structure are interchangeable under the null hypothesis, and swapping these clones forms the basis of the permutations.

In Supplemental Fig. 3B, we test the null hypothesis that fluorescence levels are equal in distribution between cells from the same environment and generation, irrespective of their clone membership. The test statistic is the average per-clone variance as in Fig. 3C, but, again, not all permutations of cell labels are allowed. Instead, cells are only permuted with other cells of the same generation.

Finally, in Supplemental Fig. 3C, we test the null hypothesis that clonal expansion and recovery are equal in distribution across different environments. The statistic is the average per-environment variance in clone size. The collection of permutations is the swapping of clones across environments.

In all cases in the present paper, we report the lower p value, $pl$⁠, approximated via Monte Carlo with $B=250,000$permutations. This tests whether the data has lower average variance, and thus greater within-group (i.e., clone or environment) relatedness, than one would expect under the null hypothesis.

The premise of the method is to label the cells under consideration with distinct combinations and concentrations of division tracking dyes, generating multiple unique fluorescence profiles (Fig. 1A). After labeling, cells are sorted according to their fluorescence profile and placed in a system of interest, such as an in vitro or in vivo environment. In concert, bulk populations of labeled cells are used to identify generation-determining gates for each unique profile (Fig. 1B). On recovery at a later time, the lineage membership, generation number, and phenotypic state of each cell can be determined by flow cytometry (Fig. 1C). As a number of division diluting dyes with distinct fluorescent spectra are commercially available, the number of the combinatorially created distinguishable profiles generated can be optimized for the system of interest. Furthermore, the use of division tracking dyes to monitor clonal lineages is a significant feature of this approach, enabling simultaneous measurement of phenotypic changes and clonal division progression.

For illustration of the method, we analyzed the in vitro differentiation of stimulated, purified murine CD8⁺ T cells at the level of individual clones. Upon activation, CD8⁺ T cells generate substantial population-level heterogeneity, which is underpinned by a significant familial component (3, 4, 22, 31, 32). Purified murine CD8⁺ T cells were labeled with different combinations of three division tracking dyes, CFSE, CTV, and CPD, resulting in 10 distinct combinations and then stimulated with anti-CD3, anti-CD28, and rhIL-2. Anti-mouse IL-2 blocking Ab was also added to remove the effect of any endogenous production (24, 29). After 24 h, just prior to their first division (24), a single founder cell from each of the fluorescence profiles was sorted and mixed into each of 29 tissue culture wells, allowing analysis of up to 10 distinct, cocultured clonal families per well. In parallel cultures, cells from each fluorescence signature were sorted into new tissue culture plates. In all cases, the cells were maintained in the same stimulatory conditions as during the initial activation period. Sixty hours after initial stimulation, cells were harvested and analyzed for division progression (Fig. 1D) and expression of CD8, CD62L, and CD25 by flow cytometry (Fig. 1E). A known number of beads were added to each well to estimate sample recovery. Pooled across wells, 156 clonal families constituting a total of 865 cells spread over four generations were recovered. The resulting data are presented in Fig. 2 and permit the concurrent visualization of clonal lineage, marker expression level, and division progression. Additional independent experiments were performed using the same stimulation conditions described in Fig. 2 and analyzed at different time points as well as experiments using CD8⁺ T cells from distinct transgenic and reporter mice under different culture conditions. Data from these further experiments including the clonal expression of the transcription factor Blimp-1 and the chemokine receptor CXCR3 are provided in Supplemental Figs. 1 and 2. Upon visualization, it is clear that complementary to previously demonstrated division synchrony (22), clones display substantial familial homogeneity. For each marker, CD8, CD62L, and CD25, the overall distribution in expression level varies over one-to-two orders of magnitude across the CD8⁺ T cell population, whereas the intraclonal distribution of expression is far narrower. This suggests that, under these stimulation conditions, a key source of phenotypic heterogeneity across a population of CD8⁺ T cells early after in vitro activation is underpinned by intraclonal concordance and interclonal variation.

The data produced by the assay has an unusual structure that necessitates careful consideration for statistical hypothesis testing. The primary concern is that clones consist of a relatively small number of cells, so statistical tests based on asymptotic results may be inappropriate. A secondary concern is that the proportion of each clone recovered from a single captive environment (culture well, animal, etc.) can result in a systemic, rather than biological, statistical coupling between cohabiting clones that must be circumvented. Thus, to complement the experimental method, we propose a choice of simple-to-implement, nonasymptotic permutation tests to interrogate the data (Fig. 3A) for a range of null hypotheses. A natural exploratory statistic based on the label-permutated data can also be plotted, providing visual cues as to the likely outcome of such tests (Fig. 3B).

An all-encompassing hypothesis would posit that each cell’s expression level is independent of its clone, generation, and environment. To test that, it suffices to compute a clonal statistic, such as the average clonal variance of fluorescence, and compare it to the distribution of the same statistic when the data are permuted by reassignment of cells to clones (Fig. 3C). If one wished to challenge the null hypothesis that the expression levels of clones, rather than cells, are independent of their environment, one can evaluate a statistic such as the average per-environment variance of expression and compare it with the same statistic in which the clones to environment labels have been permuted (Fig. 3D). The resulting p value for these, and all other permutation tests, is the proportion of permuted assignments that result in a statistic that is at least as extreme as for the true assignment (Fig. 3E).

To challenge more nuanced hypotheses, a similar procedure can be used in conjunction with suitable restrictions on the class of allowed reassignments. For example, if one suspected that recovery of clones was environment-dependent but still wished to challenge if clonal expression was independent of the environment, one cannot arbitrarily reassign clones among environments, as the test described in Fig. 3D could fail because of correlations in the level of clone recovery rather than any inherent biological environmental dependence. Instead, the desired test can be achieved by restricting reassignments across environments only to clones that are fully recovered (i.e., for which every expected cell is measured) and have the same generation structure. This works as these clones are conditioned to not be subject to sampling bias. Similar approaches can be used to test several alternate and restricted hypotheses regarding other dependencies of clonal progression and cellular phenotype (Supplemental Fig. 3A–C).

For illustration, the tests described in Fig. 3C and 3D were applied to the data shown in Fig. 2. Fig. 3F plots the CD62L expression levels of all 865 cells pooled as well as fractionated per well (i.e., per environment) and per clone, in which the latter two are rank-ordered from highest mean geometric fluorescence to lowest. Fig. 3G plots the difference between the rank-ordered geometric mean fluorescence of the data and the geometric mean fluorescence of the label-reassigned data, averaged over reassignments, as well as 95% confidence intervals under the null hypothesis. For these data, the per-well statistic consistently lies within the confidence intervals, whereas the per-clone statistic falls far outside. The statistical significance of these observations is confirmed by the hypothesis tests (Fig. 3H), demonstrating strong evidence that CD62L expression depends on clone (p < 10⁻⁵), but no evidence of per-environment dependence (p = 0.72), for this system. Equivalent analysis for CD8 and CD25 is presented in Supplemental Fig. 3D.

As this method enables identification of clonal progression and phenotypic expression, it allows for the direct measurement of asymmetric expression among sibling cells after the first division following stimulation (Fig. 4). ACD is a key driver of cellular diversity during development (33) that has been implicated in mature stem cell systems (34) as well as the adaptive immune response (35–38). To determine if ACD has occurred, it is necessary to identify cells that are siblings and to measure properties of each. This is typically challenging, as generating statistically meaningful numbers of sibling cell pairs is highly time consuming by conventional methods, such as fixed-image microscopy or live filming, but is made much more attainable with this multiplex assay. On plotting the expression levels of siblings, one anticipates distinct patterns (Fig. 4, upper panels) dependent on the following: whether the underlying biology was ACD with identifiable sibling polarity, achievable by specific ligand-receptor labeling (39) or asymmetrically segregating endocytosed fluorescent beads (40), or ACD with undetermined sibling polarity; if there were no inheritance; or if there were symmetric inheritance. For illustration, we repeated the experiment setup described in Fig. 2, using plate-bound anti-CD3 in the presence of anti-CD28 and rhIL-2, but harvested cells 42 h poststimulation to observe more clonal families after only one division event. One hundred seventy-eight clonal families with two or more members were recovered, totaling 427 cells (see Supplemental Fig. 4). Of these, 96 clones consisted of two sibling cells in generation one, allowing us to examine their expression relationships. Plots for each of CD8, CD25, and CD62L are provided in Fig. 4 lower panels and are redolent of Fig. 4 upper right panel, indicating highly symmetric divisions for this system under these stimulation conditions.

The clonal basis of T cell activation and subsequent emergence of phenotypic heterogeneity is an important focus in furthering our understanding of lymphocyte biology (22, 41, 42). Using example data sets, we have demonstrated the utility of combining multiplexed division tracking dyes with single-cell sorting and conventional flow cytometry–based phenotyping to analyze the clonal lineage properties of CD8⁺ T cells with simplicity and high throughput.

Using this method, we observed a striking and significant concordance in marker expression among the progeny of single T cell clones after standard in vitro culture. These data imply that activated founder CD8⁺ T cells have the potential to pass on a heritable, phenotype-determining program to their progeny through multiple rounds of cell division. The nature of this program, and how it is preserved to such a precise degree through numerous repetitions of the cell cycle, is unknown. The relative contribution of shared heritable fate determinants, as seen in this study in vitro, and the imposition of lineage branching points by, for example, ACD or a chance encounter with a cytokine will require further experiments tracing cells during ongoing immune responses in vivo.

A key advantage of the method is the ability to undertake direct measurement of sibling phenotype generated after the first division following stimulation. As the system can identify siblings in the presence of other accessory cells, it will be possible to systematically investigate how manipulation of the activation conditions affects the fate of each sibling in a pair. For example, it has been suggested that the synapse that forms between a dendritic cell and a T cell provides polarity cues for an asymmetric division and that this cue is further enhanced by the affinity of interaction (22, 32, 41). The method is well suited to systematically measure how such culture and stimulation variables affect concordance and fate in first-generation siblings and later generation relatives.

Existing lineage-tracing technologies have contributed significantly to the understanding of the clonal basis of many biological processes. Those methods, however, have a number of caveats that leave important aspects of biological systems unmeasured. Measuring division progression, as enabled by the approach described in this article, alleviates some of these shortcomings and allows the development of the customized statistical methodology presented in this article alongside the clonal data. These tools provide prospective users with a robust means of assessing the relative impact of clonal lineage and environmental influence on cell fate selection.

Perhaps the most significant advantage of this method is its ease of implementation. By making use of affordable, commercially available reagents and widely accessible technology, any researcher with access to flow cytometry services can easily apply this method to study clonal dynamics in their system of interest. Therefore, although we have illustrated the method in this article for in vitro T cell systems, we believe it will find wide use including application to in vivo cell-tracing systems; however, this is not yet validated. This method does not require genetic manipulation, cell infection, or breeding of fluorescent or congenic reporter systems. It can identify lineages of adherent cells in vivo or in vitro within complex cultures that include additional cell types, provided they are labeled and/or identified using compatible cell-specific markers. Consequently, it is broadly applicable and well suited to address questions of expansion and differentiation at the level of clones.

The authors have no financial conflicts of interest.