Feng Chen,Zizhang Li,Xiaoyu Zhang,Peng Wu,Wenjing Yang,Junnan Yang,Xiaoshu Chen,Jian-Rong Yang
Phylogenetic Comparative Analysis of Single-Cell Transcriptomes Reveals Constrained Accumulation of Gene Expression Heterogeneity during Clonal Expansion
phylogenetic comparative analyses,cell lineage tree,transcriptomic heterogeneity,functional homeostasis,robustness
Abstract In the same way that a phylogeny summarizes the evolutionary history of species, a cell lineage tree describes the process of clonal expansion, in which gene expression differences between cells naturally accrue as a result of stochastic partitioning and imperfect expression control. How is functional homeostasis, a key factor in the biological function of any population of cells, maintained in the face of such continuous accumulation of transcriptomic heterogeneity remains largely unresolved. To answer this question, we experimentally determined the single-cell transcriptomes and lineage relationships of up to 50% cells in single-HEK293–seeded colonies. Phylogenetic comparative analyses of the single-cell transcriptomes on the cell lineage tree revealed three lines of evidence for the constrained accumulation of transcriptome heterogeneity among cells, including rapid saturation of transcriptomic heterogeneity upon four cell divisions, reduced expression differences within subtrees closer to expression boundaries, and cofluctuations among genes. Our analyses showcased the applicability of phylogenetic comparative methods in cell lineage trees, demonstrated the constrained accumulation of transcriptomic heterogeneity, and provided novel insight into the functional homeostasis of cell populations.
Phylogenetic Comparative Analysis of Single-Cell Transcriptomes Reveals Constrained Accumulation of Gene Expression Heterogeneity during Clonal Expansion

In the same way that a phylogeny summarizes the evolutionary history of species, a cell lineage tree describes the process of clonal expansion, in which gene expression differences between cells naturally accrue as a result of stochastic partitioning and imperfect expression control. How is functional homeostasis, a key factor in the biological function of any population of cells, maintained in the face of such continuous accumulation of transcriptomic heterogeneity remains largely unresolved. To answer this question, we experimentally determined the single-cell transcriptomes and lineage relationships of up to 50% cells in single-HEK293–seeded colonies. Phylogenetic comparative analyses of the single-cell transcriptomes on the cell lineage tree revealed three lines of evidence for the constrained accumulation of transcriptome heterogeneity among cells, including rapid saturation of transcriptomic heterogeneity upon four cell divisions, reduced expression differences within subtrees closer to expression boundaries, and cofluctuations among genes. Our analyses showcased the applicability of phylogenetic comparative methods in cell lineage trees, demonstrated the constrained accumulation of transcriptomic heterogeneity, and provided novel insight into the functional homeostasis of cell populations.

In the same way that all species on earth are connected by a single phylogenetic tree, all cells within a clonal population of cells (e.g., a multicellular organism) can be described using a cell lineage tree. To understand the biology behind phenotypic differences between species, phylogenetic comparative analyses based on species trees are commonly employed and have been proved to be successful (Cornwell and Nakagawa 2017). They are, however, rarely used to investigate the intercellular phenotypic differences based on cell lineage trees. In this context, the ability of any population of normal cells to control intercellular phenotypic differences and maintain functional homeostasis is essential for their biological function (Mo et al. 2021). Specifically, it is well known that significant expression differences can be seen even among cells in a single clone, which is a major manifestation of expression noise (Raser and O'Shea 2004; Raj and van Oudenaarden 2008). In spite of the fact that such stochasticity of gene expression has been demonstrated to be genetically regulated (Munsky et al. 2012; Chen and Zhang 2016) and nonneutral to natural selection (Metzger et al. 2015; Schmiedel et al. 2019; Sun and Zhang 2020), it is unclear to what extent transcriptomic heterogeneity can be tolerated without affecting the functional homeostasis of the population. Understanding this biological constraints on transcriptomic heterogeneity may also provide valuable insight into the developmental process, since expression noise beyond the constraint triggers cellular differentiation (Chang et al. 2008; Eldar and Elowitz 2010). One way to approach this problem is to consider the simple yet fundamental scenario in which one cell continuously divides and grows into a clonal population of cells although maintaining functional homeostasis. During this process of clonal expansion, transcriptomic heterogeneity naturally accrued due to imperfect expression control (Golding et al. 2005; Raj et al. 2006; Munsky et al. 2012; Baudrimont et al. 2019) and stochastic partitioning in cell division (Huh and Paulsson 2011; Kysela et al. 2013). Meanwhile, the transcriptomic heterogeneity should also be constrained by gene expression dynamics (Golding et al. 2005; Raj et al. 2006; Munsky et al. 2012; Baudrimont et al. 2019), shared ancestry among cells (Kimmerling et al. 2016), and timely divisions (Bertaux et al. 2018). However, to the best of our knowledge, many questions regarding this simple scenario remain unanswered. For example, how much transcriptomic heterogeneity can be tolerated without affecting functional homeostasis, and when will that level be reached during the clonal expansion? Is heterogeneity regulated independently for each gene and cell, or is there some kind of coordination involved? By addressing these questions in the simple scenario of clonal expansion, one will not only reveal key parameters in a general framework for the growth of all types of cell population but also help understand specific clonal expansion processes, such as the developmental process by which single T-lymphocytes generate a variety of progeny that are both acutely responsive to a specific antigen and protective in the event of a future exposure (Chang et al. 2007; Arsenio et al. 2014, 2015; Flossdorf et al. 2015). To answer these questions, one potential strategy would be to measure the transcriptomic heterogeneity at different stages of clonal expansion, as has been done in developmental studies (Petropoulos et al. 2016). However, since transcriptome measurements routinely require RNA extraction and therefore killing the cells, such experiments would need to be conducted independently for clones grown at different stages. Consequently, the comparability of data from different expansion stages may have been compromised due to differences in progenitor cells, environmental heterogeneity across experiments, and the potential errors introduced by merging data across time points (Tritschler et al. 2019). In this regard, recent technological advances that have allowed for the reconstruction of cell lineage trees (e.g., scGESTALT (Raj, Wagner, et al. 2018) and CARLIN (Bowling, et al. 2020)) have provided an alternative approach. Specifically, the different stages of clonal expansion can be represented by groups of cells with varying kinship distances (extractable from the lineage tree), and these groups of cells have progenitors of minimal differences and grow in the same environment. Moreover, cell lineage trees can be subject to phylogenetic comparative analyses, which have been proven useful in analyzing species trees. With these considerations, we reconstructed cell lineage trees of colonies seeded by individual HEK293 cells. The lineage positions and single-cell transcriptomes of up to 50% of the cells in the colonies were determined. This data set facilitated phylogenetic comparative analyses of single-cell transcriptomes during clonal expansion with an unprecedented resolution and revealed noise saturation after four divisions, decelerated accumulation of transcriptomic heterogeneity in subtrees rooted with cells inferred to have expression levels near an upper limit, and cofluctuation among gene pairs whose dosage balance is functionally essential. Our results showcased the applicability of phylogenetic comparative methods in cell lineage trees and demonstrated the constrained dynamics of transcriptomic heterogeneity during clonal expansion, a previously understudied aspect of functional homeostasis despite stochasticity in gene expression.

Transcriptomic heterogeneity studied here is, by definition, about the differences between individual cells. It is therefore natural to anchor the analysis to the time scale of cell divisions. In a full cell lineage tree containing all cells that appear during clonal expansion, all division events are represented as bifurcations on the tree. In terms of retrieving the cell division events, the experimental reconstruction of single-cell lineage trees using Cas9-based genome editing and 10× Chromium is highly informative, but not ideal. A major problem lies in the fact that the vast majority of cells were lost during clonal expansion, sample digestion, and library preparation, and the loss of a cell results in one missed cell and one missed cell division on the lineage tree. Moreover, this issue is difficult to solve by modeling the dynamics of Cas9 editing, which is highly variable and uncoupled from the division of cells (Feng et al. 2021). We reasoned that keeping as many cells in the colony as possible on the tree should alleviate (but not resolve) this problem and therefore attempted to reconstruct for the clonal expansion process a lineage tree with a high cell sampling rate by adapting the scGESTALT protocol (Raj, Gagnon, et al. 2018; Raj, Wagner, et al. 2018) with some modifications. These include 1) the introduction of up to three mismatches between sgRNAs and the target site to reduce editing efficiency and hence the probability of multisite long deletions in the lineage barcode (fig. 1A and supplementary table S1, Supplementary Material online. See also Materials and Methods); 2) the use of an inducible Cas9 to facilitate synchronized initiation of colony expansion and editing of the lineage barcode (fig. 1A and supplementary fig. S1AandB, Supplementary Material online); 3) the close monitoring of the expression levels of transgenes (supplementary fig. S1C–F and table S2, Supplementary Material online) and the morphology of the colony (fig. 1B and supplementary fig. S2, Supplementary Material online) to allow selection of colonies that are most likely to exhibit high proportions of surviving cells with diversely edited lineage barcodes (supplementary figs. S3 and S4, Supplementary Material online); 4) a four-step digestion protocol to minimize the number of cells killed during digestion (fig. 1C and supplementary fig. S5, Supplementary Material online. See also Materials and Methods); and 5) the use of PacBio HiFi sequencing for more accurate calling of the lineage barcode allele (fig. 1C and see Materials and Methods). We obtained two HEK-293T colonies (hereinafter “sample 1” and “sample 2”) in Cumate-containing medium, each initiated by a single cell carrying the lineage tracing system. The two colonies grew to ∼5,000 cells and ∼35,000 cells, respectively, maintaining their identity as HEK-293T cells (fig. 1B and supplementary fig. S6, Supplementary Material online) before they were digested. Single-cell cDNA libraries were then generated for each colony using a GemCode Single Cell Instrument (10× Genomics) (see Materials and Methods). Each library was split into two halves: one half was subjected to NovaSeq-based RNA-seq for single-cell transcriptomes, and the other half was subjected to amplification of the lineage barcode, followed by PacBio Sequel-based HiFi sequencing of the lineage barcode (fig. 1C and also see supplementary tables S3 and S4, Supplementary Material online). A lineage tree was then constructed for each sample based on the recovered lineage barcode alleles using the maximum likelihood (ML) method (see Materials and Methods). For sample 1(fig. 1D), the lineage tree contained 2,667 single cells and exhibited a low frequency of long deletions (8.8% of all unique alleles), a high percentage of single-cell tips (91.7%), and a high depth (the maximum number of ancestor nodes for a tip is 18) and sampled the majority (2,667/4,900 ≈ 54.4%) of cells present in the colony at digestion. The lineage tree of sample 2 (supplementary fig. S7, Supplementary Material online) was of similar quality, with the exception that only 6.9% of the cells in the colony were sampled (supplementary tables S4 and S5, Supplementary material online. See Materials and Methods). We conducted a couple of tests for the lineage tree to ensure that the hierarchical relationships between cells were largely correct. On the one hand, our bootstrap analysis (see Materials and Methods) showed that the internal nodes in the tree of sample 1 received median bootstrap support of 82%, which is comparable to previous reports (Liu et al. 2021). On the other hand, we assumed that ancestral barcode alleles could be recovered in descendent cells as yet-to-decay transcripts, so that shared ancestral alleles could be used as independent evidence for close kinship between cells. Using PacBio reads from each sampled cell, we identified ancestral alleles by extracting reads corresponding to a subset (but not all) of mutations found in the lineage barcode allele of the focal cell (see fig. 1E for an example). As expected, the ancestral barcode alleles were particularly abundant in cells with high expression of lineage barcodes (fig. 1F). Most importantly, we found shared ancestral alleles in ∼75% of the closest cell pairs (i.e., their immediate ancestor is the same internal node) in sample 1 (fig. 1G). As the barcode transcript decays, the probability of sharing an ancestral allele should decrease as cell pairs become more distant, which is indeed observed (fig. 1G). In sample 2, the internal nodes were supported by approximately 85% of bootstraps (supplementary fig. S7A, Supplementary Material online), and shared ancestral alleles were detected in 25% of the closest cell pairs (supplementary figs. S7BandC, Supplementary Material online), as expected based on the lower cell sampling rate. Overall, we concluded that the cell lineage tree adequately captures the hierarchical relationships/kinship among the cells. Nevertheless, we emphasize that the reconstructed lineage trees are still too coarse for detailed developmental analyses, such as those conducted on Caenorhabditis elegans (Sulston et al. 1983; Du et al. 2014, 2015).

The gradual accumulation of heterogeneity in expression during clonal expansion should result in greater average transcriptomic differences in larger subtrees/subclones, as they are at a later stage of clonal expansion. To test this, we compared the heterogeneity of expression within a subtree with the height of that subtree (i.e., the number of nodes on the longest branch in the subtree; see fig. 2A and Materials and Methods). Note here that the height of a subtree can also be viewed as a proxy for the divergence time of any cell pairs whose latest common ancestor is the root of the subtree. On the level of individual genes, we evaluated the expression heterogeneity by the coefficient of variation (CV) of a gene's expression level among single cells at the tip of a subtree. Such estimated expression heterogeneities for the top 100 highly expressed genes turned out to be gradually increasing in typical subtrees in sample 1 with heights of 1, 2, and 3 (figs. 2B and C). Similarly on the transcriptome level, a principal component analysis involving all genes indicated that subtrees with larger heights showed greater variation in gene expression (fig. 2D). To obtain a more comprehensive view, the average transcriptomic heterogeneity was calculated for every subtree by the Euclidean distance of expression levels between all pairs of single-cell tips that cannot be paired up in a smaller subtree (because if they do, they should be assigned to the smaller subtree and therefore an earlier stage of clonal expansion). Such calculated transcriptomic heterogeneity has a positive correlation with the height of the subtree (sample 1: fig. 2E; Spearman's rank correlation coefficient ρ = 0.32, P < 10−8. Sample 2: fig. 2F; ρ = 0.12, P = 0.02). This result demonstrated the gradual accumulation of transcriptomic heterogeneity during clonal expansion. As the transcriptomic heterogeneity continues to accrue, it should eventually disrupt the population's homeostasis, which contrasts with common observations for the HEK-293 cell line and suggests that noise accumulation will eventually reach saturation. In order to verify this notion, we fitted the growth of transcriptomic heterogeneity among single-cell tips within a subtree to the height of the subtree (i.e., the divergence time between the cells) using a segmented linear model (Muggeo 2008) (see Materials and Methods; hereafter named as “saturation” model; red line in fig. 2E and F). This model allows a saturated stage beyond a certain divergence time and is often used to identify saturated phenotypic divergence in phylogenetic comparative analyses (Muschick et al. 2014). We found that this saturation model significantly outperformed other linear models assuming either unlimited noise accumulation along clonal expansion (“unlimited” model) or complete independence between noise and lineage distance (“independent” model) (P < 0.005, pairwise likelihood ratio test; see Materials and Methods). According to the saturation model, transcriptome heterogeneity was saturated at the height of 3.35 and 2.02 for sample 1 and 2, respectively. The lower saturation height of sample 2 is expected by its lower cell sampling rate compared with sample 1 (see Materials and Methods, and fig. 2G). Using data from sample 1 and assuming that the height of the lineage tree corresponds with the rounds of cell divisions, we further estimated that noise saturation requires four rounds of cell divisions/generations in HEK293 colonies (fig. 2G and see Materials and Methods). As far as we are aware, this is the first experiment-based estimation of the time required for saturation of transcriptomic heterogeneity in clonal expansion. To further resolve this pattern of noise saturation, we analyzed the accumulation of expression noise for each gene using the saturation model. An individual gene is considered to exhibit noise saturation if, in the segmented linear (i.e., saturation) model fitted to the experimental data, 1) the slope of the first stage of heterogeneity accumulation is significantly positive and 2) there is a breakpoint (i.e., the intercept of the second stage is significant). We found that 4,254 genes in sample 1 and 3,694 genes in sample 2 displayed patterns of noise saturation, with a significant overlap of 933 genes (fig. 3A, hypergeometric P < 10−10. Overlapping genes listed in supplementary table S6, Supplementary Material online). Furthermore, a correlation was observed between sample 1 and sample 2 in terms of the divergence time at which noise reaches saturation for individual genes (fig. 3B, Pearson's R = 0.091, P = 0.005), as well as the saturated level of expression heterogeneity in terms of CV (fig. 3C, R = 0.74, P < 10−159). Consistent with theoretical expectations, we found that genes with rapid noise saturation encode mRNAs with shorter half-lives (Schueler et al. 2014) (supplementary fig. S8A, Supplementary Material online) and exhibit a lower static level of expression noise (supplementary fig. S8B, Supplementary Material online). There is a significant variation in saturation time across genes (supplementary table S6, Supplementary Material online). For example, some genes, such as the Ribosomal Protein L7 (RPL7), showed noise saturation at height ∼ 2. Some genes, such as Zinc Finger Protein 91 (ZNF91), showed much later noise saturation at height ∼ 5. These results indicate that the noise saturation of individual genes during clonal expansion is also readily detectable with the lineage tree and that at least partially the saturation of transcriptomic heterogeneity can be attributed to constrained noise accumulation of individual genes during clonal expansion. To gain a deeper understanding of the functional characteristics of these 933 genes that exhibited noise saturation, we conducted a Gene Ontology (GO)–based functional enrichment analysis. It was found that these genes are significantly enriched with genes that are involved in ribosome biogenesis (fig. 3D) and genes that are components of the ribosome (fig. 3E). This result is not merely caused by the high expression and therefore greater statistical power of the ribosome-related genes, because a similar functional enrichment is observed for genes best described by the saturation model (compared with the unlimited and the independent model) in both sample 1 and sample 2 (supplementary fig. S9, Supplementary Material online), despite that high expression should enhance statistical power similarly for all three models. The ribosome-related genes’ phenomenon of noise accumulation across cell divisions and saturation later can be understood from two different perspectives. The ribosome-related genes are usually highly expressed and are located in broad regions of open chromatin (Bajic and Poyatos 2012), which makes them relatively resistant to transcriptional bursting, a major source of expression noise (Dar et al. 2016). It is therefore likely that the expression heterogeneity of ribosome-related genes will be limited within a cell cycle, leaving the possibility for further accumulation of heterogeneity due to stochastic partitioning over successive divisions of the cell. On the other hand, as a major factor driving extrinsic noise (Swain et al. 2002; das Neves et al. 2010; Taniguchi et al. 2010; Johnston et al. 2012), noise saturation of ribosome-related genes is also compatible with the previous finding that natural selection favors the reduction of global expression noise (Metzger et al. 2015; Sun and Zhang 2020). In summary, these results indicate that noise accumulation is subject to biologically regulated constraints.

What is the exact nature of the biological constraint limiting the accumulation of transcriptomic heterogeneity during clonal expansion? Phenotypes were often found to have a bounded distribution in phylogenetic analyses (Koch et al. 2004; Kaiser et al. 2007; Donovan et al. 2011). Similarly, we hypothesized that expression boundaries for individual genes could be one mechanism of noise constraint, and we tested this hypothesis using a simple geometric model for gene expression dynamics (fig. 4A and see Materials and Methods). Briefly, the transcriptome of each cell can be viewed as a moving particle in an N-dimensional space (N equals the number of genes, in this case the 933 genes listed in supplementary table S6, Supplementary Material online). Under this model, inheritance dictates that the transcriptome of all tips within a subtree should be limited by an N-sphere approximately centered at the root of the subtree, that is, a limited range of expression variation (the circles in fig. 4A). Whenever an expression boundary is present, a subtree whose roots are closer to the boundary should exhibit a smaller range of expression variation than another subtree whose roots are farther away from the boundary (e.g., gene 1 in fig. 4A). It is important to note that, although there should be both lower and upper boundaries for expression, we have only focused on the upper boundary because the lower boundary is most likely undetectable given the overwhelming technical noise associated with scRNA-seq data at the lower end of the expression spectrum. For each gene, we inferred the upper boundary of expression based on its maximal observed expression within our scRNA-seq data. Comparison with other published scRNA-seq data (Raj, Wagner et al. 2018; Bowling et al. 2020) suggested that the inferred expression boundaries are unlikely physical limits, but instead biologically regulated limits (supplementary fig. S10, Supplementary Material online). The gene's expression at the root of a subtree were further inferred by hierarchically averaging its expression at the tips (see Materials and Methods). We then estimated a subtree root's distance from the boundary expression by calculating the relative difference between the expression of a gene at the root and at the boundary, followed by geometrically averaging these differences across all genes. We found that roots closer to the boundary expression tend to give rise to subtrees with transcriptomically less heterogeneous/noisy (fig. 4B). By limiting the analysis to highly expressed genes, whose expression measurements are likely more accurate, the correlation was strengthened (dots with error bars in fig. 4C). Moreover, we found that this pattern of less-noisy tips in subtrees with roots close to the boundary expression is dependent on the real lineage relationship between the cells, as it disappeared when we randomly rearranged the cells on the lineage tree (violin plots in fig. 4C). In addition, tips within subtrees whose roots are near the boundary expression tend to display weaker tendencies of expression divergence, as they exhibit smaller angular distances of RNA velocity (fig. 4D). To further support the existence of a constraint on expression heterogeneity, we fitted our data with the Ornstein–Uhlenbeck (OU) model (see Materials and Methods), which is commonly used in phylogenetic studies to infer mean-reverting tendencies of quantitative traits (Butler and King 2004). As a result, we noticed that the 933 genes with noise saturation were enriched with genes showing strong mean-reverting tendencies (50% or 20% genes with top ɑ value in OU models, respectively shown by solid and dashed lines in fig. 4E), a pattern that becomes increasingly dramatic when we considered only saturating genes with higher expression levels (x axis in fig. 4E). Collectively, our findings suggest that, at least for some genes, expression levels are individually constrained, resulting in limited transcriptomic heterogeneity.

Gene expression constraints shown above limit expression heterogeneity individually for each gene and are independent of other genes. However, another type of constraint, cofluctuation between genes (fig. 5A), is also possible. It should be noted that cofluctuation is different from coexpression, since cofluctuation refers to correlated gene expression among single cells, although coexpression refers to that among bulk samples. To detect cofluctuation among gene pairs, we calculated the Pearson's correlation coefficient for each pair of the 17,281 human genes with nonzero expression in our data set. With a threshold of P < 0.05 by Pearson's correlation test after Benjamini–Hochberg (BH) correction of multiple testing, we identified over 1.2 million gene pairs that were cofluctuating in both samples 1 and 2 (fig. 5B inset). These gene pairs represent a significant overlap between the two samples (hypergeometric P < 10−300), which becomes even stronger under more stringent criteria for the correlation (fig. 5B). Similarly to previous studies (Sun and Zhang 2019, 2020), our results demonstrated the existence of cofluctuations in gene expression levels using scRNA-seq data, which should have contributed to the constraint of transcriptome heterogeneity. In addition, the lineage tree may facilitate further explorations of the molecular mechanisms underlying the observed cofluctuations. More specifically, expression cofluctuation can be caused by two nonmutually exclusive mechanisms: transcript cosegregation during asymmetric division (fig. 5C) and transcriptional coregulation (fig. 5D). Transcript cosegregation, in this case, is the correlated biased distribution of transcripts during cell division resulting from molecular packaging/segregation (e.g., transcripts within one mitochondrion), which leads to two sister cells with very different transcript abundances. Accordingly, the tips of the two subtrees rooted by these pair of sister cells may exhibit expression levels clustered around their respective roots, thereby showing signals of cofluctuation between genes (fig. 5C). Alternatively, gene pairs with cofluctuation maintained by transcriptional coregulation are likely to exhibit strong correlations even during individual cell cycles (fig. 5D). For genes that are sensitive to changes in relative dosage ratios, such as those belonging to the same protein complex (Fraser et al. 2004), this feature is particularly beneficial, because too little or too much of one subunit of a protein complex can compromise the proper assembly of the complex and waste the energy invested in producing the other subunits (Fraser et al. 2004). We reasoned that with the knowledge of the lineage tree, transcriptional coregulation can be discriminated from cosegregation using phylogenetically independent contrasts (PICs) (Felsenstein 1985), because coregulation is analogous to correlated evolution of quantitative traits, whereas cosegregation is analogous to correlated traits due to shared ancestry. Using PIC-based correlations that is significant after BH correction for multiple testing, we found approximately 0.5 million pairs of genes with signal of transcriptional coregulation in both samples (fig. 5E, inset), which is a significant overlap (hypergeometric P < 10−300) that is further enlarged under more stringent criteria for the PIC-based correlation (fig. 5E). Compared with cofluctuations implied by significant correlation in the original expression levels, the additional requirement of significant PIC correlation can enhance the enrichment of gene pairs within the same protein complex (fig. 5F) and remove gene pairs known to encode cosegregating transcripts (fig. 5G). A closer look at the data further revealed some specific cases with apparent biological relevance. For example, OS9 Endoplasmic Reticulum Lectin (OS9) is known to interact with Heat Shock Protein 90 Beta Family Member 1 (HSP90B1, also named as GRP94) to form a chaperone–lectin complex and participated in the endoplasmic reticulum–associated degradation quality control pathway (Christianson et al. 2008; Määttänen et al. 2010; Seidler et al. 2014; Di et al. 2016). The functional relationship of these two genes is consistent with the strong cofluctuation inferred by correlation from both PIC and the original expression level (fig. 5H, left). On the contrary, the mRNA of OS9 and Vacuole Membrane Protein 1 (VMP1) are colocalized on the outer mitochondrial membrane (Fazal et al. 2019) and should be cosegregating during cell division. They appeared cofluctuating according to the correlation of original expression levels, but not PIC (fig. 5H, right). Finally, a significant PIC correlation (at adjusted P < 0.05) was observed for approximately 35.8% of gene pairs with significant OEL correlation (at adjusted P < 0.05), suggesting that transcriptional coregulation may account for a third of the observed cofluctuation among genes. As the biological signal is most likely obscured by technical noise, we conclude that transcriptional coregulation is likely the major mechanism driving cofluctuation among genes and, consequently, constraining the transcriptomic heterogeneity during clonal expansion.

To investigate the accumulation and constraint of transcriptomic heterogeneity during clonal expansion, we utilized a scGESTALT-based assay to determine the single-cell transcriptomes as well as the lineage trees of two single-cell-seeded colonies consisting of ∼5,000 and ∼35,000 HEK-293 cells, respectively. Our data set revealed rapid transcriptomic and gene-specific saturation of expression heterogeneity, biological boundaries of gene expression, and pervasive cofluctuation between different genes, suggesting a highly constrained regulatory landscape of gene expression noise. To the best of our knowledge, this is the first comprehensive study of the regulation of transcriptomic heterogeneity among cells that takes into consideration their kinship, an underappreciated factor in previous studies of expression noise using scRNA-seq. Our results thus significantly enhance our understanding of how gene expression stochasticity is constrained in the context of ongoing cell divisions. Our results have a few limitations worth discussing. First, we estimated the progress of clonal expansion of each subclone represented by a node within the lineage tree based on the height of the subtree. In this case, several types of errors can occur, including 1) underestimating the height due to inevitable cell losses during the growth of a single-cell clone or library preparation using 10× Chromium technology, 2) overestimating due to false editing events caused by sequencing errors in the lineage barcode, and 3) incorrect topology due to parallel editing/mutation or sequencing errors. To minimize the impact of these errors, we optimized the experimental protocol to retain as many cells as possible, used PacBio HiFi reads to ensure highly accurate readout of lineage barcodes, and chose an estimate of the subtree's height that was robust to underestimation due to cell losses (see Materials and Methods). More importantly, as the height of the subtree increased, the probability of sharing ancestral barcode alleles decreased and expression heterogeneity increased, indicating that the hierarchical relationship/kinship between the cells was generally accurate. Second, the single-cell transcriptome quantified with 10× Genomics Chromium technology exhibits a high level of technological noise. For instance, droplet-based scRNA-seq has shown to have a CV of at least 0.5 for the counts per million (CPM) values of the majority of genes (Zhang et al. 2019). Nevertheless, the patterns revealed in our analysis, as well as by other researchers (Sun and Zhang 2019), suggest that such data may still capture some major trends of expression noise, such as cofluctuations between genes (Sun and Zhang 2019). Rather, our results may have been undermined by the technological noise, such that there are even stronger biological signals at work behind our observations. Third, we did not distinguish between extrinsic noise, which is caused by a difference in cell state, and intrinsic noise, which is caused by stochasticity during gene expression in cells of the same state. We, however, suspected that the heterogeneity constraints were mostly for extrinsic noise, since both expression boundary and cofluctuation likely require extrinsic regulatory factors. It will be interesting to further disentangle intrinsic and extrinsic noise using methods such as allele-specific scRNA-seq data (Sun and Zhang 2020). Our study represents a technological effort in a novel direction of scGESTALT-like assays for the reconstruction of cell lineage trees, in which there is more emphasis on the sampling density in terms of single cells, instead of the breadth of sampling (Raj, Wagner et al. 2018; Bowling et al. 2020). Thus far, the value of the nonsampled branches/cells in the developmental cell lineage tree has largely been ignored. But just as the ghost lineages (ancient extinct species or unknown modern species) in the phylogenetic tree can strongly influence multiple phylogenetic interpretations (Tricou et al. 2022a, 2022b), retrieving the nonsampled cells may provide valuable new insights. For example, observing transcriptomic heterogeneity saturation at four cell divisions would have been impossible in lineage trees with cell sampling density < 6.25% (=1/24), since all cell pairs are expected to be at least four divisions apart in such trees. Similarly, the dense sampling of single cells should enhance the statistical power of analyzing cofluctuations between genes (Ackerly 2000). Furthermore, in terms of approximating divergence times between individual cells, high-density sampling for cells might be a useful alternative to statistical inference based on assumed editing efficiency of Cas9. The reason for this is that Cas9 editing efficiency is much more variable than the distribution of loss cells, especially when cells are sampled at sufficiently high density. It is possible to improve experimentally for higher sampling densities by using longer lineage barcodes for a larger recording capacity, cell cycle-coupled editing to improve division event resolution, and using other genomic editing systems instead of CRISPR/Cas9 (which causes DNA double-strand breaks and, thus, DNA damage–induced cell death) (Liu et al. 2021). Regardless, the present study demonstrated a unique advantage of dense lineage trees in analyzing transcriptomic heterogeneity among related cells. How the accumulation of transcriptomic heterogeneity is constrained during clonal expansion or cellular population growth is critical not only for HEK-293 cells but also for other cell types whose functions are dependent upon functional homeostasis of the cellular population. The technological limitations notwithstanding, our results provide a clear answer: transcriptomic heterogeneity among cells is biologically constrained at levels ranging from individual genes (expression boundary) to interactions among genes within the same cell (cofluctuation), to populations of cells (noise saturation) (fig. 6). Furthermore, these results indicate potential applications for assessing the differentiation of cells and thus defining cell types objectively (Vickaryous and Hall 2006; Arendt et al. 2016; Yuan et al. 2020) based on an unconstrained accumulation of expression heterogeneity among cells, as well as identifying candidate strategies for manipulating cell fate based on the expression boundaries of individual genes or cofluctuation between genes. In the future, it will be interesting to extend our strategy to assess the differentiating and nondifferentiating lineage trees of the same type of cells, or to compare the lineage tree underlying the expression geometry among normal cells with that among cancerous cells (Liang et al. 2015). “Nothing in biology makes sense except in the light of evolution” (Dobzhansky 1973). Our study demonstrated the applicability of phylogenetic comparative analyses to single-cell transcriptomes on cell lineage trees. With the continued development of relevant experimental methods and the accumulation of more cell lineage tree data sets, phylogenetic comparative methods will likely facilitate further insights into clonal expansion or other biological processes that can broadly be categorized as somatic evolution. In addition to being a great opportunity for evolutionary wisdom to find its next application, it is also a chance to re-examine many evolutionary theories in a novel type of phylogeny—the cell lineage tree.

We generated over 1 million randomized 20 bp sequences as sgRNA candidates. Following this, candidates with no more than three substitutions relative to any fragment of the human genome were removed for further analysis in order to avoid off-targeting the genome. We selected the spacer sequence 5′-TATTCGCGACGGTTCGTACG-3′ from the remaining sgRNA candidates, which we refer to as sgRNA1. Then, 13 protospacer sequences were designed based on sgRNA1 according to the following criteria: 1) each protospacer contained 2 to 3 bp of mismatches with sgRNA1, which reduces the probability of multisite long deletions (Quinn and Jones 2021); 2) there is no recurrence of any sequence of 9 bp or longer across protospacers; and 3) consecutive repeats of the same nucleotide for more than 2 bp were completely absent. The 13 protospacers along with protospacer-adjacent motif (PAM) sequences designed following the previous GESTALT strategy (McKenna, et al. 2016) were concatenated in an order of decreasing cutting frequency determination (CFD) scores (Doench et al. 2016) to obtain the full lineage barcode. To further ensure that the lineage recording capacity of the lineage barcode was fully utilized, we designed sgRNA2/3/4, which targeted the 9th/12th/13th protospacer sequences (rarely edited sequences in preliminary experiments) with 100% identity but showed small CFD scores (<0.55) with other protospacers. The designed lineage barcode and sgRNA sequences are listed in supplementary table S1, Supplementary Material online. To obtain a lineage tree with a high cell sampling density, it is imperative that the lineage barcode can be successfully reverse transcribed during the preparation of a single-cell library using 10× Chromium technology. Several approaches were used to enhance the capture rate of the lineage barcode by 10× Chromium technology, which includes 1) using the strong promoter EF1α, 2) concatenating the lineage barcode to the 3′ end of the EGFP gene so that green fluorescence could be used as a marker for single-cell colonies with appropriate expression of the lineage barcode, and 3) adding a 20 nt poly-dA (A20) sequence at the 3′ end of the lineage barcode to facilitate capture by the poly-dT RT primer on 10× gel beads. Along with other accessory elements, the final lineage barcode vector (LV-EF1A>EGFP:T2A:Bsd:V1. VectorBuilder, no. VB170906-1089zqd) was constructed by using Gateway technology (Hartley 2003; Chee and Foan 2015) on the backbone of LV-3d. After the transfection of the lineage barcode vector into HEK-293 cells (see below), we applied strong selection by blasticidin (15 µg/ml) for 12 days to ensure a high expression level of the lineage barcode. For the sgRNAs, two DNA cassettes, U6>sgRNA1-U6>sgRNA2 and U6>sgRNA3-U6>sgRNA4, were synthesized and subsequently combined by Golden Gate ligation in the sgRNA vector LV-U6>sgRNA1-U6>sgRNA2-U6>sgRNA3-U6>sgRNA4 (VectorBuilder, no. VB180515-1178wxn).

The two lentiviral particles of the barcode and sgRNA vectors were individually transfected into HEK-293FT cells (Cyagen, no. HEKFT-30001) with lentivirus packaging plasmids using Lipofectamine 2000 (ThermoFisher, no. 11668027). To obtain the viral particles, the supernatant was harvested three days posttransfection, in which cellular debris was removed by centrifugation at 3,000 rpm for 15 min and filtering with a 0.45 µm disposable needle filter (Millipore, no. SLHP033RB), followed by centrifugation at 50,000 × g for 90 min. To gauge the viral titration, the viral particle–containing supernatant was diluted to various concentrations, which were then used to transduce a second batch of HEK-293FT cells. The genomic DNA (gDNA) of these HEK-293FT cells was isolated using DNeasy Blood & Tissue Kits (Qiagen, no. 69504). The gDNA was further analyzed by qPCR targeting the barcode and the endogenous gene BMP2 using the QuantiNova SYBR Green PCR kit (Qiagen, no. 208057-500T) so that the viral titration could be determined. A HEK-293 cell line (abm, no. T3352) with an inducible Cas9 was used to reconstruct the lineage tree of the clonal expansion of a single-cell-seeded colony. Specifically, Cas9 was driven by an inducible CuO promoter (CuO-Cas9) and its suppressor CymR was inserted into the genome. Therefore, Cas9 could be activated by adding Cumate to the cell culture (Kaczmarczyk et al. 2013; Li et al. 2017), enabling clonal expansion and lineage barcode editing to start simultaneously. These cells were then transduced with lentiviral particles containing lineage barcodes at a low MOI (0.1), followed by selection with 15 µg/ml blasticidin (InvivoGen, no. ant-bl-1) for 12 days, giving rise to HEK-293 cells harboring a single copy of the lineage barcode. The sgRNA lentiviral particles were transduced at different MOIs (10, 20, and 30), and positive cells were selected with 400 µg/ml hygromycin (InvivoGen, no. ant-hg-1) for 5 days. PCR or RT-qPCR (supplementary fig. S1A–E, Supplementary Material online) was used to detect the expression levels of Cas9, the lineage barcode, and sgRNA1 using the primers listed in supplementary table S2, Supplementary Material online. The inducibility of Cas9 was confirmed by Western blotting (supplementary fig. S1F, Supplementary Material online) for cells cultured in the presence of different concentrations of Cumate (abm, no. CH065). The experiments hereinafter were all conducted in the constructed HEK-293 cells containing inducible Cas9 and the lineage barcode/sgRNA.

Approximately 106 HEK-293 cells were grown in 10% serum media containing 30 µg/ml Cumate in a 6-well dish. After three days of culture, cells were harvested for gDNA extraction using DNeasy Blood & Tissue Kits (Qiagen, no. 69504). Alleles of the lineage barcode were amplified from gDNA using Phanta Max Super-Fidelity DNA Polymerase (Vazyme, no. P505) with the gDNA-V1f and gDNA-V1r primers (supplementary table S2, Supplementary Material online) and cloned into the CE Entry Vector (Vazyme, no. C114). Recombinant plasmids were obtained by blue-white selection, and 20 randomly selected clones were subjected to Sanger sequencing to assess the editing efficiency. To assess the feasibility of clonal expansion starting from a single cell, we prepared 384 single-cell-seeded colonies with HEK-293 cells (LV-sgRNA MOI = 30). These single cells were then grown in three 96-well plates for 10 days, during which 10 µl of CloneR (STEMCELL, 05888) was added to the medium of each colony on the 0th and 1st days, although 5 and 30 µg/ml Cumate were added to the media on the 0th, 2nd, 4th, 6th, and 8th days. The colonies were visually inspected and categorized as follows: normal (type I, 43 colonies), disassociated/substantial death (type II, 162 colonies), or not-growing/excessive death (type III, 179 colonies) (supplementary fig. S2, Supplementary Material online). Despite a relatively low success rate, these results suggest that a HEK-293 colony seeded by a single cell is possible. To gauge the editing efficiency of the lineage barcodes in single-cell-seeded colonies, 288 single cells were first manually selected by micromanipulation from the log-phase population of HEK-293 cells constructed by using three different titrations of sgRNA lentiviral particles (LV-sgRNA, MOI = 10, 20, and 30). Similar to the cultivation procedure described above, we added 10 µl of CloneR and 30 µg/ml of Cumate to the medium of each colony at different time points (see above) and allowed the colonies to grow for 10 days. The colonies were visually inspected and categorized as follows: normal (type I), disassociated/substantial death (type II), or not-growing/excessive death (type III). Thirty type I colonies with different expression levels of sgRNA and Cas9 were chosen. Then, approximately 50 ng of gDNA was harvested from each colony using the QIAamp DNA Micro Kit (Qiagen, no. 56304). The gDNA was PCR amplified in a 50 µl reaction targeting the lineage barcode, whose editing status was further confirmed by Sanger sequencing (using 13 µl of the product), TA cloning (1 µl of the product), or Illumina HiSeq-PE250 sequencing (36 µl of the product) as follows. Based on the Sanger sequencing data, we found that ∼93% of the lineage barcode alleles from these single-cell-seeded colonies were edited (supplementary fig. S3A, Supplementary Material online). As determined by TA cloning and subsequent Sanger sequencing of individual lineage barcode alleles, colonies that had high copy numbers of sgRNA in the genome (transformation with MOI = 20 or 30) and were grown with a high concentration of Cumate (5 or 30 g/ml) displayed a reasonably high diversity of alleles (supplementary fig. S3B–G, Supplementary Material online). The sufficient editing efficiency and diversity of the lineage barcode was further corroborated by the sequencing data from Illumina HiSeq-PE250. Briefly, the clean data of each sample sequenced by HiSeq-PE250 were trimmed with fqtrim (https://ccb.jhu.edu/software/fqtrim/) using the default parameters. Reads 1 and 2 were spliced with flash (Magoč and Salzberg 2011) using 30 bp of overlapping sequence and 2% mismatches. Nonspecifically amplified sequences generated during PCR were removed by Bowtie2 (Langmead and Salzberg 2012) with alignment to the human reference genome using default parameters. The sequences that did not contain the two expected sub-sequences of "TCTATACGATCCGGACCTC" and "ATTCAGAGCAGTGCTCG" (two mismatches allowed) were removed as falsely amplified products. Finally, the sequences of the lineage barcode were aligned with the reference sequence by MUSCLE (Edgar 2004) using the default parameters, and the editing events in each sequence were called according to the previously proposed method in the GESTALT paper (McKenna et al. 2016) (supplementary fig. S4A–H, Supplementary Material online). For the sample (supplementary fig. S4H, Supplementary Material online) with the highest sgRNA copy number (MOI = 30) and strongest Cumate induction (30 µg/ml), ∼82.9% of the deletion events revealed by HiSeq were shorter than 20 bp (supplementary fig. S4I, Supplementary Material online), and majority of the lineage barcode alleles contained no intersite deletions (supplementary fig. S4J, Supplementary Material online). More importantly, the frequency of editing events suggested that editing started as early as the two-cell stage of the colony (supplementary fig. S4K, Supplementary Material online), and intersite deletion events were absent until the 16-cell stage (supplementary fig. S4L, Supplementary Material online). These results suggest that the recording capacity of the barcodes has been effectively utilized in our experiment.

We developed an optimized protocol for single-cell digestion to minimize cell loss, which is detailed as follows. First, the culture medium was removed by gentle suction, and morphologically normal colonies (type I colonies remaining after suction, according to visual inspection under an optimal microscope) were retained. Second, without washing the single-cell colonies with any solution, we directly added 100 µl of ACCUTASE (StemPro, no. A1110501) to digest the colonies at 37°C for 20 ∼ 30 min and continuously examined the status of the single-cell clones after 20 min. Third, the colonies were disassociated into single cells via gentle, repetitive pipetting, after which 900 µl of PBS was added to each colony to terminate the digestion. Finally, each 1 ml single-cell suspension was centrifuged at 500 × g for 10 min. The supernatant was removed by gentle suction and discarded, and the remaining single cells were resuspended in 20 µl of PBS containing 5% BSA. As a final batch of single-cell-seeded colonies to be used for scRNA-seq, 576 colonies each seeded by one single HEK-293 cell (LV-sgRNA MOI = 30) were cultured for 10 days (see the previous section). We discarded type II and type III colonies and colonies with weak GFP fluorescence (i.e., lineage barcode expression), resulting in 15 remaining colonies that were subjected to further digestion as described in the previous paragraph. Altogether, we obtained single-cell suspensions for eight samples (among 15 colonies at the beginning of digestion). We used three of these samples to determine the fraction of living single cells by trypan blue staining and cell counting with Countstar BioMed IM1200, which indicated 98.49% (among ∼13800 cells), 96.74% (∼5480 cells), and 83.08% (∼8530 cells) living cells, respectively (supplementary fig. S5, Supplementary Material online).

We chose five single-cell suspension samples and loaded them individually onto a GemCode Single-Cell Instrument (10× Genomics, Pleasanton, CA, USA) to obtain the GEM-RT product following the Single Cell 3′ Reagent Kit v2 User Guide. Half of the GEM-RT product was directly amplified without being fragmented, targeting the lineage barcode (Sc-barcode-f and Sc-barcode-r listed in supplementary table S2, Supplementary Material online), and the product was later subjected to TA cloning followed by Sanger sequencing of ∼20 lineage barcode alleles to obtain a rough estimation of the diversity of the lineage barcode. The two samples showing the highest diversity of lineage barcodes were chosen and sequenced with PacBio Sequel technology via standard library preparation and sequencing protocols. The remaining half of the GEM-RT product of each of these two chosen samples was cleaned, fragmented, amplified, and used to prepare a sequencing library with an average fragment size of 450 bp. Library sequencing was finally performed on the Illumina NovaSeq platform following a standard 10× scRNA-seq protocol.

For the HEK-293 samples assayed in this study, raw scRNA-seq reads were obtained from the Illumina Nova-seq platform. Mapping and quantification were performed with the reference GRCh38 using Cell Ranger (Zheng et al. 2017) according to the official 10× Genomics protocols, which resulted in a count table (of unique molecular identifiers or UMIs) for each gene in each droplet/cell. Cells with <200 detected genes were removed. As a result, 2,863 and 20,893 cells were retrieved from 10× scRNA-seq in sample 1 and sample 2, respectively (supplementary table S4, Supplementary Material online). Based on the ratio between the area of the whole clone and the area of a cell, sample 1 contained 4,900 cells (supplementary fig. S6, Supplementary Material online). In other words, expression profiles were obtained for 58.4% (2,863/4,900) of the cells in the colony. Assuming the same rate of single-cell transcriptome retrieval for the parallel colonies, we estimated there are ∼35,775 cells in sample 2, whose cell density is too high for area-based estimation. Expression data were normalized for sequencing depth using the NormalizeData function in the Seurat R package (Butler et al. 2018). For the analysis of the physical or biological boundaries of expression, raw expression matrixes of HEK-293 cells and multiple other scRNA-seq data sets (Raj, Wagner, et al. 2018; Bowling et al. 2020) generated with Cell Ranger were combined and normalized according to sequencing depth by Seurat. Under the geometric model of expression levels, the boundary of expression for each individual gene was estimated by the maximal observed expression level among all single cells. Then, for each cell, the distances to the expression boundary for all genes were combined according to their geometric average (x axis of fig. 4B).

To distinguish the heteroduplex DNA formed during PCR, all subreads from each zero-mode waveguide (ZMW) well on the PacBio Sequel platform were first categorized as belonging to the positive or negative strand according to the existence of 12 successive A or T bases, respectively, within the raw sequence (allowing three mismatches except in the leading three A/T bases). Consensus sequences were generated for both positive and negative strands by using the CCS program (Anzalone et al. 2020) from PacBio, with a minimum of three passes/subreads per strand. The consensus sequences were then aligned to the human genome by using Bowtie2 (Langmead and Salzberg 2012) with default parameters, and the sequences with significant hits were removed because they represented nonspecific amplification of the lineage barcode primers. We further checked the existence of four sequences ("TGCTGCCCTCTGGTTATGTG", "AGATCGGAAGAGCGTCGTGTAG", "TCTATACGATCCGGACCTC", and "ATTCAGAGCAGTGCTCG") within the consensus sequence, allowing two mismatches for each primer. Next, the cell barcode and UMI were extracted from the consensus sequence and were refined according to the cell barcode and UMI in the transcriptome data. That is, if the cell barcode and UMI from Sequel or NovaSeq differed by only one nucleotide, they were considered to be the same. Finally, the lineage barcode sequences were extracted from the remainder of consensus sequences between the fragments of "TCTATACGATCCGGACCTC" and "ATTCAGAGCAGTGCTCG". Detailed summary statistics of the read numbers throughout the analyses presented above are listed in supplementary table S3, Supplementary Material online. To further correct the sequencing error on the PacBio Sequel Platform, the lineage barcodes with the same cell barcode and UMI were merged by multiple sequence alignment via MUSCLE (Edgar 2004), followed by selecting the nucleotide with the highest frequency at each site. Then, the lineage barcodes were aligned with the reference sequence by MUSCLE using the default parameters, and the editing events of each sequence (i.e., the lineage barcode allele) were called according to a previous method employed by GESTALT (McKenna et al. 2016). Then, for each lineage barcode allele, the frequency was calculated as the total number of UMIs of the allele and its ancestral allele. Here, the ancestral allele of a specific allele was defined as any allele in which the observed editing events were a subset of the editing events in the focal allele; that is the focal allele could be generated by further editing of the ancestral allele (fig. 1E). Finally, the lineage barcode allele of a cell was defined as the allele with the highest frequency, prioritizing the alleles with more editing events if the frequencies were equal. We retrieved 2,667 (54.4% of all cells within the colony of sample 1) and 2,480 (6.9% of all cells within the colony of sample 2) cells from the two samples, which have both single-cell transcriptomes and lineage barcodes (supplementary table S4, Supplementary Material online). Finally, all recovered lineage barcodes of the two samples were reconstructed into two multifurcating lineage trees using ML implemented via the IQ-TREE LG model (Umkehrer et al. 2020). To process our data with IQ-TREE, we replaced the default amino acid substitution matrix in IQ-TREE with a matrix containing integers from 1 to 190, which represented the weights of editing events. Accordingly, all lineage alleles were encoded in a “sequence” of “amino acids,” with each “residue” in the sequence representing one editing event, and the ancestral and derived states of each “residue” were chosen according to the aforementioned substitution matrix, such that rare editing events received small weight and common editing events received large weight (McKenna et al. 2016). Finally, we obtained a lineage tree of which >89% of the leaves were single cell. The detailed statistics during the above data processing steps are listed in supplementary table S4, Supplementary Material online. Programming codes for reproducing our analyses can be found on GitHub. Additionally, using the same lineage allele encoding and mutation weighing scheme, we also tried reconstructing the cell lineage tree via maximum parsimony (MP) and neighbor-joining (NJ) methods. Due to the pervasive polytomies in all trees, we compared the ML trees with the MP and NJ trees using Nye similarity (Nye et al. 2006), which accounts for fractional similarity between “split” by internal nodes. As a result, we found that MP and NJ support the majority of topologies in the ML-based phylogenies. Specifically, Nye similarity between the ML tree and NJ tree is 0.69 for sample 1 and 0.73 for sample 2, whereas Nye similarity between the ML tree and the MP tree is 0.58 for sample 1 and 0.58 for sample 2.

Assuming that the subtrees were unbiased representations of clonal expansion in the same environment, a quantitative metric was necessary to put each subtree in the context of a specific phase of clonal expansion, allowing the accumulation of transcriptomic heterogeneity to be compared with clonal expansion progress. It would seem intuitive to use the number of tips (or cells in the tips) of the subtree as such a metric. However, this may be inappropriate due to the large number of noncaptured cells. Instead, we utilized the height of the subtree, defined as the number of nodes in the longest branch (fig. 2A), which has several advantages. First, the height of the subtree should be proportional to the number of rounds of cell divisions that a clone has undergone, which is more appropriate for our analysis within the context of cell divisions. Second, the height of the subtree was more robust against loss of cells than tree size (number of tips) in estimating the progress of clonal expansion. Due to the sheer number of small subtrees, the impact of errors in height estimates should be averaged out for the smaller subtrees, although for larger subtrees, the effects should be negligible, as the height should be accurate as long as there are only a minor number of intact branches (with few cell losses). Third, increasing the sampling density of cells, which is a significant feature of our experimental pipeline, may reduce errors for heights, but not for those based on the assumption of Cas9 editing efficiency of lineage barcodes. Regarding the estimation of divergence time between any two cells, for the same reasons discussed above, and to limit the number of metrics used, instead of alternative metrics such as topological distance (number of nodes along the path connecting two tips), we resorted to the height of the subtree rooted at the most recent common ancestor of the focal cells (fig. 2A). Nevertheless, we acknowledge that there are fundamental limitations in estimating divergence time based on Cas9 editing that is not coupled with or more frequent than cell division. For example, even with methods more sophisticated than tree height or topological distance, it is still impossible to resolve cell division during which no Cas9 editing occurred. In the future, it may be desirable to use a lineage tracing system with a higher temporal resolution, such as SMALT (Liu et al. 2021), to resolve the problem.

In order to determine whether noise accumulation will eventually reach saturation, a segmented linear model that allows a saturation stage above a certain height (the saturation model) (Muschick, et al. 2014) was used. Briefly, the average transcriptomic heterogeneity was calculated for every subtree by the Euclidean distance between all pairs of single-cell tips that cannot be paired up in a smaller subtree (because if they do, they should be assign to the smaller subtree and therefore an early stage of clonal expansion). These values from all subtrees of the same height were pooled together, and their median value was fitted to the height using either a simple linear model (using the R function “lm.fit”) allowing a nonzero slope to reflect unlimited noise accumulation or a linear model with a fixed slope = 0 to reflect noise levels independent of the height or the segmented linear model (using the R package “segmented”) to reflect noise accumulation with a saturation stage. Using the “AIC” and “lrtest” functions in R, the performance of these three models was compared. We noticed that when fitting noise data for individual genes, multiple fitting trials could sometimes result in different saturation heights (i.e., the “breakpoint”). As a means of ensuring the robustness of our conclusion, we performed 20 repeated fittings with each data set and excluded genes whose 20 fitted breakpoints did not round to the same integer. Approximately 1% of the genes were removed as a result of this procedure.

In order to demonstrate the impact of cell sampling density on the analysis of cell lineage trees, we simulated sparse cell sampling by down-sampling the cells in sample 1. More specifically, we randomly selected a predetermined proportion of tips from sample 1 and dropped the remaining tips and corresponding internal nodes although maintaining the tree's topology. Essentially, this simulates a lineage tree from sample 1 but with a lower cell sampling rate. Using the same computational pipeline as we analyzed sample 1, we estimate the height of the subtree at which transcriptomic heterogeneity reaches saturation. Based on the results of this procedure repeated for 1,000 down-sampled lineage trees, we found that a lower cell sampling density would result in smaller saturation heights (fig. 2F, cyan points). It is also noteworthy that the difference between saturation heights in samples 1 and 2 (3.35 and 2.02, respectively) is also in agreement with their respective sampling densities (54.4% and 6.9%, respectively). Most importantly, the results also demonstrate that noise saturation could not be observed at lower cell sampling densities. With a similar logic, we aimed to estimate the actual saturation height if the cell sampling rate was raised to 100%. For a rough simulation of this scenario, we randomly insert new tips as sisters of existing tips in sample 1 (i.e., the original sisters become cousins). The transcriptome of the new tips was either copied from the sister cells (30% of the new tips) or cousin cells (70% of the new tips) of the original tip. We randomly generated 1,000 “up-sampled” lineage trees, which showed noise saturation at a height of 4.62 ± 0.25 (fig. 2G, red point).

Functional annotation and GO term enrichment analyses of genes with signal of expression noise saturation were performed by the BiNGO (Maere et al. 2005) plugin of Cytoscape (Shannon et al. 2003), with the default parameters. Note that the GO-based functional enrichment was conducted over two gene sets, including 933 genes displaying a significantly positive slope in the first segment and significant intercept in the second segment (existence of a breakpoint) in the saturation (segmented linear) model in both samples (fig. 3D and E), and 360 genes best described by the saturation model (compared with the unlimited and independent model) in both samples (supplementary fig. S9, Supplementary Material online).

For each cell, the RNA velocity vectors for the expression level of each saturated gene was estimated by the “velocyto.R” packages (La Manno et al. 2018) in R. Here, the RNA velocity vector of one or more genes can be seen as the difference between their current and future expression levels, the latter of which was estimated by the amount of unspliced transcripts (La Manno et al. 2018). Consider the simple geometric model for gene expression dynamics in figure 4A, assuming that the RNA velocities of all single cells in a subtree represent the tip of a continued arrow of divergence starting from the center of the circle (the root of the subtree); because not all directions are possible in subtrees near the boundary (the yellow or blue circles in fig. 4A), the average angular distance between RNA velocities should be smaller therein compared with subtrees farther away from the expression boundary (the red full circles in fig. 4A). The angular distance of the RNA velocity was calculated for each cell pair using “angle” function of “matlib” package in R.

For each gene, its normalized expression levels in all individual single cell tips, in combination with the overall lineage tree, were processed by the “ouch” package (Butler and King 2004) in R to fit the OU model. To demonstrate the relationship between the OU model and the signal of noise saturation as detected by the segmented linear model, we simulated changes in gene expression levels along the cell lineage tree of sample 1 using either the Brownian motion or OU model implemented by the “fastBM” function in the R package “phytools”(Revell 2012), using different mean expression levels (supplementary fig. S11, Supplementary Material online, x axis). Without loss of generality, we set the expression at the root to the specified mean expression level, and the ɑ and σ parameters of the OU model to the average values of the genes in our data set. We found that under the BM model (no mean-reverting tendency), none of the genes are best described by the saturation model compared with the unlimited or independent model (the green bars in supplementary fig. S11, Supplementary Material online). Additionally, 9% to 13% of the genes exhibit noise saturation under the OU model. More importantly, the fraction of genes with noise saturation found in simulation is quite similar to the empirical observation (8.5% ∼ 18.2%) made with our data set. There results suggest that the simple mean-reverting behavior (i.e., the OU model) may be sufficient to explain our observation of noise saturation.

Take the enrichment of cofluctuating gene pairs in protein complex members as an example, we constructed a 2 × 2 contingency table: by categorizing each pair of genes (only those detected by scRNA-seq) into one of the four cells in the table on the basis of 1) whether the pair of genes showed correlated expression (by Pearson's correlation coefficient of either the original expression level or phylogenetic independent contrasts (PIC) of the expression level, see below) and 2) whether both genes in the pair were members of the same protein complex. The fold enrichment was calculated as (A/(A + C))/((A + B)/(A + B + C + D)). In a similar manner, fold enrichment analyses such as coregulation or cosegregation versus PIC-based correlation between genes, as well as genes with noise saturation versus genes with mean-reverting tendency (larger ɑ value in the fitted OU model), were conducted.

PIC were calculated in R using normalized expression levels of individual genes and the lineage tree. Note that the linear models fitted to find correlation between PICs of different genes are required as going through the origin. But the same requirement did not apply to correlations for the original expression levels. The list of human protein complexes and their members were obtained from the CORUM database (Giurgiu et al. 2019). To examine whether cofluctuations inferred by PIC can eliminate cosegregating gene pairs, we also extracted the list of genes localized to the outer mitochondrial membrane (OMM) from the previous study (Fazal et al. 2019), which was defined as mRNAs with log2(number of reads in OMM APEX-seq/number of reads in negative controls) > 1.5.

Click here for additional data file.