Single-Cell RNA-Sequencing Reveals a Continuous Spectrum of Differentiation in Hematopoietic Cells

1Hematopoietic stem cells (HSCs) have the ability to self-renew and produce cells that give rise to all different blood cell types. Our understanding of the functional properties of these various hematopoietic cell types has been advanced mainly by population level analysis. Current methods of purifying hematopoietic cells to relative homogeneity are based on the expression of specific combinations of cell surface markers. However, a homogeneous population of cells, as determined by a well-defined set of cell surface markers, may include many functionally distinct populations. This was nicely illustrated in studies showing that within the HSC compartment, individual HSCs may have different reconstitution patterns (e.g., balanced production of myeloid and lymphoid cells or deficiency in lymphoid potential). More recently, it was demonstrated that common myeloid progenitors (CMP) are a mixed population of cells with distinct lineage potentials. The lack of CMPs as a separate cell entity with broad myeloid potential brings into question the traditional model of hematopoietic lineage development and further underscores the importance of revising the current view of lineage development in hematopoiesis. Therefore, there is a need to address the exact composition of the stem and progenitor populations in vivo, as well as the relationships between them. Single cell transcriptome analysis might provide answers to these outstanding questions.

2Among vertebrate models, the zebrafish provides a unique combination of advantages for the study of blood development at the single cell level. Zebrafish blood contains cells of all hematopoietic lineages and orthologs of most transcription factors involved in mammalian hematopoiesis. Importantly, transcriptional mechanisms and signaling pathways in hematopoiesis are well conserved between zebrafish and mammals, making them a clinically relevant model system.

3Over the past few years, a number of transgenic zebrafish lines were generated in which hematopoietic cell specific promoters drive expression of fluorescent molecules. These reporter lines provide a valuable resource of labeled cells ranging from HSCs to a wide range of mature blood cell types. As in mammals, adult hematopoiesis in zebrafish is both continuous and asynchronous. Thus, a single sample of kidney marrow (the analogous tissue to mammalian bone marrow) contains the full spectrum of hematopoietic cell types at various stages of differentiation at any one time. As this is the single site of hematopoiesis in zebrafish, and is easily accessible, the cells are minimally perturbed when sorted ex vivo, making this an ideal system to study basic principles of regulation of differentiation, both at the molecular and cellular levels.

4Here we used high-throughput single-cell RNA sequencing combined with fluorescence-activated cell sorting index sorting analysis of adult zebrafish marrow-derived hematopoietic cells. We ordered cells by their progression through differentiation based on gene expression profiles using no prior knowledge of which cell population they belong to, as defined by surface markers. Our analysis revealed the continuous nature of thrombocyte development and the coordinated transcriptional programs that govern differentiation progression. Interestingly, thrombocytes in zebrafish remain transcriptionally active even after leaving the kidney marrow and entering the circulation.

5Here, we used single-cell RNA-sequencing (RNA-seq) of zebrafish kidney cells to resolve the cellular hierarchy of lineage development in the myeloid branch of hematopoiesis. To focus on this lineage, we used expression of CD41 as a marker of HSCs and the megakaryocyte equivalent in fish (“thrombocytes”). CD41 in human is highly regulated during hematopoietic development, and in zebrafish, the Tg(cd41:EGFP) reporter line labels two distinct populations of cells that express the cd41-EGFP transgene. The weakly fluorescent (EGFPlow) subset marks HSCs and progenitor cells, and the brightly fluorescent (EGFPhigh) subset includes mature and differentiated thrombocytes.

6Using flow cytometry, we identified EGFPlow and EGFPhigh cells and sorted 188 cells from each population from a single kidney from a Tg(cd41:EGFP) reporter fish (Figure 1A; Figures S1A–S1I). Each EGFP+ cell was collected in a single well of a 96-well plate, and for each cell, its size (FSC), granularity (SSC), and EGFP fluorescence level were recorded. Single-cell mRNA-seq libraries were constructed and sequenced to a depth of around 2.5 million reads per library. Of 376 cells, 13 cells failed our quality control (QC) and were removed from further analysis (Experimental Procedures; Figures S2A and S2B). For the remaining 363 cells, we accurately quantified between 1,000 and 6,000 genes per cell.

7To identify groups of cells and order them in terms of their developmental progression, we used a multi-step approach. First, we used independent component analysis (ICA) to identify distinct factors that describe the variability of EGFP cells. ICA revealed four latent factors (hidden variables) that explain (1) a progression among EGFPlow cells (“within_small_component”), (2) a switch from EGFPlow cells toward EGFPhigh cells (“difference_component”), and (3) progression among the EGFPhigh cells (“within_large_component”). Finally, the fourth factor identified three outlier cells (“outlier_component”) (Figure S3A).

8To facilitate data depiction, we used non-linear dimensionality reduction (t-distributed stochastic neighbor embedding [t-SNE]; ) to represent the four latent factors in two dimensions (Figure 1B). ICA revealed a clear distinction between EGFPlow and EGFPhigh cells, implying sharp divergence at the transcriptional level (Figure S3A; Figure 1B).

9In addition, EGFPlow cells are a more heterogeneous group compared to EGFPhigh cells. To explore this further, we used hierarchical clustering to partition EGFP cells based on their independent components (Figure S3B). Interestingly, whereas EGFPlow cells were split into four distinct clusters (here named 1a, 1b, 2, and 3), EGFPhigh cells were all grouped into a single cluster (here named 4), confirming the substantial heterogeneity of the EGFPlow population of cells (Figure 1C).

10Differentiation of hematopoietic cells involves the acquisition of specific phenotypes that depend on the repression of genes characteristic of a multipotent cell state and expression of lineage-restricted genes. Thus, the whole process can be conceptualized as a temporal ordering of a highly coordinated transcriptional program through which each cell progresses. To examine the transcriptional transitions undergone by cd41-EGFP cells during differentiation, we ordered cells based on the cluster they belonged to, the latent factor that explains the variability of the cells within the cluster, and the level of EGFP fluorescence (details provided in the Experimental Procedures). Our model assumes gradual changes in gene expression during developmental progression of thrombocytes along a one-dimensional (i.e., non-branching) path. (We could not detect any apparent branch point in the data.) This ranking of cells through the entire process was treated as “pseudotime.”

11To ensure our pseudotime ordering was stable, we also ordered the cells using an alternative method, a Bayesian Gaussian process latent variable model (; see Experimental Procedures). Comparing the paths these orderings take when regressed into the t-SNE depiction, one can appreciate the similarity between them (Figure 2A). The two pseudotime orderings agreed very strongly (Spearman correlation 0.97; Figure 2B), giving us confidence in our method.

12When presented in pseudotime, the expression of endogenous cd41 (also known as itga2b) and EGFP, as well as EGFP fluorescence, recorded during sorting, were highly correlated and showed an expected increase through pseudotime (Spearman rho 0.85, 0.80, and 0.82, respectively) (Figure 2C). This supports our pseudotime ordering of the cells from the HSC to the differentiated thrombocyte extracted from a single kidney.

13To define the identity of cell types within the five clusters, we evaluated the expression of orthologs of transcription factors and other genes known to be relevant in mammalian hematopoiesis, including the expression of early (cd61, also known as itgb3a/b) and late (cd42b, also known as gp1bb) markers of megakaryocyte differentiation (Figure 3). The panel of genes analyzed was representative of HSCs (Tal1, Lmo2, Lyl1, Gata2, Runx1, Meis1, C-myb, and Erg;,,,, ), megakaryocyte/erythroid (Fli1, Gfi1b, Gata1, Cd61, Cd42b, Vwf, and Selp;,,, ), and myeloid- (Gfi1, Pu.1 also known as spi1a/b, and Cebp1;, ) lineage-affiliated genes.

14For each gene, we assessed the level of its expression in pseudotime, as well as the fraction of cells that expressed the gene of interest in each of the clusters (Figure 3). For example, c-myb was highly expressed in cluster 1a, as well as in clusters 1b, 2, and 3, but was downregulated in cluster 4. This is in line with previous reports that C-myb is expressed in immature hematopoietic cells and is downregulated during differentiation. Cells in cluster 1a had relatively high expression of lmo2, tal1, and meis1. These genes, together with fli1, showed a similar distribution of expression across pseudotime, whereas gata2 was more restricted to cluster 1a. The mammalian HSC genes runx1 and erg were expressed at a relatively low level overall, and in a small fraction of cells within all clusters. Overall, most of the mammalian HSC marker genes examined are expressed in cluster 1a, and to a lesser degree in 1b, 2, and 3.

15In contrast, Gata1 and Gfi1b are known to be expressed at high levels in the erythroid and megakaryocyte lineages but not in HSCs. In our dataset, gata1a and gfi1b were expressed in all clusters except cluster 1a. Furthermore, expression of both early (itgb3a/b) and late (gp1bb) markers of megakaryocyte differentiation started very early and peaked late in pseudotime, confirming that more mature thrombocytes are largely confined to cluster 4.

16We also assessed the expression of two well-known platelet genes, vWf (von Willebrand factor) and selp (P-selectin), through pseudotime (Figure 3). Our analysis revealed that, contrary to previous reports, thrombocytes in zebrafish do not express von Willebrand factor and P-selectin. This was confirmed by qPCR analysis of cd41 EGFPhigh thrombocytes from zebrafish kidney. We found, however, that vWf was expressed in the whole kidney sample and in fli1:GFP positive cells sorted from Tg(fli:EGFP) fish, suggesting that the vWf expression pattern differs somewhat in zebrafish compared to mammals.

17Surprisingly, myeloid lineage-affiliated genes (e.g., spi1, gfi1, and cebp1) were largely absent across all cells (Figure 3). This suggests that there is no common myeloid progenitor population in this dataset, which charts a continuous HSC to thrombocyte pathway. Altogether, our data are consistent with the notion that cells from cluster 1a belong to HSCs that transition directly to erythroid-thrombocyte progenitor cells, possibly circumventing the CMP step. Although this is surprising, there are other reports of direct, unconventional, HSC to megakaryocyte-erythroid progenitor transitions, such as a recent report in mouse.

18Identification of these progenitor and differentiated cell types prompted us to carry out additional analyses of the sets of genes that strongly correlate with the individual cell types. We used a machine learning method, random forest feature importance, to find genes whose expression “marks” distinct clusters of cells. The unique sets of genes expressed in each of the cell types provide an opportunity to reveal novel markers of the identified cell types, and at the same time, provide more insight into their biological function.

19Among the numerous newly identified cell-type markers (Table S1), we found several of particular interest (Figure 4A). Ccr9a is a member of the beta chemokine receptor family and is known to be expressed in HSCs; our data show that ccr9a expression is highly correlated with cluster 1a (Figure 4B). Transcription elongation factor A (SII), tcea3, was specifically expressed in cluster 1b (Figure 4B). Cells from cluster 1b can also be sorted by combining expression of plasminogen receptor gene (plgrkt) and ascc1 (Figure 4B). Good marker genes for cluster 2 included e2f8, which encodes a protein involved in progression through the cell cycle and top2a, a DNA topoisomerase involved in processes such as chromosome condensation and chromatid separation (Figure 4B). Interestingly, the overrepresented gene ontology (GO) enrichment terms for cluster 2 included cell division and cell cycle (Figure 4A), suggesting that an expansion phase precedes lineage commitment and terminal differentiation of thrombocytes.

20To experimentally validate the prediction of greater proliferation in this progenitor population, we sorted cells from clusters 1a/1b/2 versus 3 and 4, by distinguishing these three populations based on EGFP fluorescence, and SSC and FSC (Figures S4A–S4G). We compared the cell cycle distributions of the sorted populations using propidium iodide (PI) staining. The combined cells from clusters 1a/1b/2 had a significantly higher proportion of cells in S and G2/M phase compared to clusters 3 and 4 (Figure 4C), validating our finding that these cells proliferate faster.

21These results show that expression of EGFP together with SSC and FSC values could be used to efficiently separate cells from clusters 3 and 4 from the early progenitor populations (1a/1b/2) in the cd41 reporter line (Figures S4 and S5). Additional markers for cluster 3 included combined high expression of fzd8b and no expression of mibp (Figure 4B). For cluster 4, a high level of cd41 uniquely marks this population.

22Finally, we also assessed a unique set of genes expressed by the three outlier cells. GO enrichment analysis of their marker genes yielded only three statistically significant GO terms, all linked with immunity (Figure 4A). One plausible explanation is that these outlier cells represent macrophages that have engulfed or are attached to thrombocytes and hence retained a high level of EGFP fluorescence. Indeed, the outlier cells expressed an array of macrophage/monocyte affiliated genes such as mpeg (macrophage expressed gene 1), csf1r (colony-stimulating factor 1 receptor), csf3r (colony-stimulating factor 3 receptor) etc. Furthermore, compared to all other cells, the outlier cells had remarkably high FSC and SSC values, characteristic of macrophages (Figure 4D).

23Importantly, we validated many of our findings in a second set of single cell transcriptomics experiments on kidney cells, as well as circulating cells, from another fish. We sorted an additional 92 cells from cluster 1a/1b/2 (named here EarlyEnriched), 46 EGFPlow cells and 46 EFPhigh cells from the kidney of another Tg(cd41:EGFP) fish. We also sorted 24 EGFPlow and 68 EGFPhigh circulating cells from the same fish (Figure S6A). Our analysis confirmed that the pattern of ICA follows the same structure as observed in the previous experiment (Figure S6B). This means that the cell populations and their relative relationships are conserved in this biological replicate. Similarly, the pseudotime ordering of EarlyEnriched, EGFPlow, and EGFPhigh cells in the kidney recapitulated patterns we identified in the initial experiment (Figures 5A–5F).

24In addition, we discovered that EGFPhigh cells in circulation are transcriptionally identical to EGFPhigh cells in the kidney, with no significant change in the number of expressed genes (Figure 5B), RNA content (Figure 5C), or any gene’s expression pattern (likelihood ratio test, corrected for multiple testing with Holm-Sidak). We concluded, therefore, that the thrombocytes exit the kidney in a fully mature state and are maintained in a transcriptionally active state in circulation.

25In both datasets, the total number of genes and total mRNA content expressed per cell were correlated with its differentiation state (Figure S7). This was not due to a difference in the sequencing depth or cell size (Figure S7). Instead, it represents a biological difference between cells during development. This supports the idea that more differentiated, post-mitotic cells (clusters 3 and 4) have a specialized transcriptional program with expression of a small, focused set of genes (Figure S7).

26To find genes with similar trends in expression across pseudotime, we used a mixtures of hierarchical Gaussian processes model to cluster the pseudotime series. We identified 130 genes that are dynamically expressed through pseudotime. Clustering of these genes revealed three distinct patterns of their progression during differentiation (Figure 6A; Table S2). Genes upregulated early in pseudotime and then downregulated later (group I) were significantly enriched with the GO term “nucleic acid binding” and “chromosome maintenance” (Figures 6B and 6C; Table S2), possibly reflecting the increased proliferation of cells earlier in pseudotime. Genes gradually downregulated through pseudotime (group II) were highly enriched with the GO terms “eukaryotic translation elongation,” “ribosomes” etc. (Figures 6B and 6C; Table S2). Expression of these genes was highly correlated with the general trend of decreased RNA content over pseudotime (Spearman rho = 0.85), therefore suggesting a regulatory loop between the total RNA content in the cell and expression of genes that encode proteins relevant for ribosome synthesis. Finally, genes upregulated early and then maintained at a high level (group III) were highly enriched with the GO terms “ECM-receptor interaction,” “platelet aggregation,” and “hemostasis,” pointing to the genes important for thrombocyte function (Figures 6B and 6C; Table S2). Taken together, our analysis suggests that differentiation of thrombocytes is governed by coordinated transcriptional programs that limit the proliferation of cells and their translational capacity while simultaneously promoting genes relevant for thrombocyte function.

27Gene duplication is a common event in eukaryotic genomes and due to the teleost-specific genome duplication around 26% (i.e., 3,440) of zebrafish genes are duplicated. Gene duplicates that originate from genome duplication are called ohnologs. To assess the use of duplicated genes during thrombopoiesis in zebrafish, we examined the expression of ohnologs in each of the 363 cd41:EGFP cells. Of ∼8,000 ohnolog categories, we looked at 3,034 ohnolog categories that have only been duplicated once (ohnolog gene pairs). Of these 3,034 ohnologs, 2,107 were not expressed in our dataset. However the remaining 927 pairs can be divided into the following three major groups: (1) expression of ohnologs is mutually exclusive in individual cells (n = 177) (Figures 7A and 7B; Table S3). In this group, the expression of any one ohnolog appeared to be an independent event with an equal probability of happening. This suggests selective activation or silencing of these ohnologs in individual cells; (2) only one of the ohnologs is expressed in all cells (n = 430), Figures 7A and 7B; Table S3), and (3) both ohnologs are equally expressed in all cells (n = 218), (Figures 7A and 7B; Table S3). No patterns of ohnolog use over pseudotime were observed.

28Here we show the power of single cell transcriptome analysis to decipher the kinetics of hematopoietic lineage development. We ordered cd41 cells by their progression through differentiation based on gene expression profiles. Our analysis illustrates the continual nature of this process, where cells progressively transit through five transcriptional states that result in the generation of mature thrombocytes.

29Interestingly, myeloid lineage-affiliated genes were largely absent across all cells, suggesting direct HSC to thrombocyte-erythroid progenitor transition. The model of hematopoiesis generated recently, using single cells from over ten hematopoietic populations in mouse, implies that the megakaryocyte-erythroid lineage is closely linked to long-term repopulating HSCs and separates early from the lympho-myeloid lineage. The identification of platelet-primed stem cells within vWf-expressing long-term HSCs further confirmed that commitment to the megakaryocyte lineage starts in the most primitive stem cell compartment. Although in our dataset vWf was not expressed in any of the identified cell populations, the low expression of some of the thrombocyte lineage-affiliated genes in cluster 1a suggests that using our sorting strategy we are possibly capturing thrombocyte-primed stem cells. Therefore, HSCs in cluster 1a may represent a biased subpopulation within the wider pool of hematopoietic stem/progenitor cells present in the zebrafish kidney. Nevertheless, the gradual transition of cells during thrombocyte lineage development that we see in our dataset (e.g., gradual changes in the total number of genes as well as the total mRNA content) suggest that we do capture a continuous spectra of cells and that the common myeloid stage is not an obligatory step during thrombopoiesis.

30We also show that although each of the identified transcriptional states was characterized by substantial heterogeneity in the expression of the key lineage regulators, the underlying transcriptional program was highly coordinated. It included the simultaneous increase in the expression of genes important for thrombocyte function and suppression of genes relevant in cell proliferation and ribosomal biogenesis. Interestingly, although the maturation of thrombocytes was completed in the kidney, they maintained a transcriptionally active state in circulation. We did not, however, detect any qualitative or quantitative difference in the gene expression between circulating and kidney-based EGFPhigh thrombocytes. Surprisingly, unlike mammalian platelets, which have abundant expression of vWf, thrombocytes in zebrafish do not express vWf. Instead, our analysis suggests that other cells within the kidney marrow, such as endothelial cells (fli1:GFP positive cells), express vWf in zebrafish. Finally, we assessed use of duplicated genes during thrombopoiesis in zebrafish and identified patterns of their expression that would not be possible using a bulk transcriptomics approach.

31We used single-cell RNA-seq of zebrafish kidney cells to resolve the cellular hierarchy of lineage development in the myeloid branch of hematopoiesis and propose a refined model of developmental progression of hematopoietic cells.

32Our study addresses some of the basic questions of regulation of differentiation, both at the molecular and cellular levels. In this study, we focused on zebrafish thrombocyte development; however, a similar approach could be used in other systems and cell types.

33The maintenance of wild-type (Tubingen Long Fin) and transgenic zebrafish Tg(cd41:GFP) lines were performed in accordance with EU regulations on laboratory animals, as previously described.

34A single kidney from heterozygote Tg(cd41:EGFP) or wild-type fish was dissected and carefully passed through a strainer using the plunger of a 1 ml syringe. In the follow-up experiment, circulating GFP-positive cells were collected from the dissected heart of the same fish. Cells were collected in cold 1× PBS/5% fetal bovine serum. The kidney of a non-transgenic line was used to set up the gating and exclude autofluorescent cells. Dead cells were excluded based on PI staining. Individual cells were sorted using a Becton Dickinson Influx sorter with 488- and 561-nm lasers and collected in a single well of a 96-well plate containing 2.3 μl of 0.2% Triton X-100 supplemented with 1 U/μl SUPERase In RNase inhibitor (Ambion). At the same time, information about cell size and granularity and the level of the fluorescence were recorded. Whole transcriptome amplification and library preparation was performed using the Smart-seq2 protocol, with ERCC spike-in controls added at the same time as the oligo-dT and dNTP mixture. Twenty-five PCR cycles were performed during the amplification.

35GFP-positive cells from Tg(cd41:EGFP) kidney suspension were sorted using a Mo-Flo XDP (Beckman Coulter) with 488-, 561-, and 640-nm lasers. Cells were centrifuged at 1,200 rpm for 10 min at 4°C, resuspended in 100 μl 1× PBS and fixed by adding 300 μl ethanol. Cells were fixed overnight at 4°C, washed twice in 1× PBS, and re-suspended in 500 μl PI solution (25 μg/ml PI, 0.1% Triton X-100, 0.1% sodium citrate). Cells were incubated for 3 hr with RNase A (Sigma) and analyzed by BD LSR Fortessa (Becton Dickinson). Data were analyzed using FlowJo software.

36Sorted EGFP-positive cells were concentrated by cytocentrifugation at 350 rpm for 5 min onto SuperFrostPlus slides using a Shandon Cytospin 3 cytocentrifuge. Slides were fixed for 3 min in methanol and stained with May-Grünwald Giemsa (Sigma) as described elsewhere. Images were captured as described elsewhere.

37GFP-positive cells from Tg(cd41:EGFP) and Tg(fli1:EGFP) kidney suspensions were sorted using a Mo-Flo XDP (Beckman Coulter), along with an equal number of viable cells from the whole kidney, into 75 μl RLT buffer (QIAGEN) containing 1% β-mercaptoethanol. mRNA was extracted using Oligo (dT)25 Dynabeads (Ambion) and cDNA was prepared using SuperScript VILO (Invitrogen), according to the manufacturers’ instructions. qPCR reactions were performed using the 7900HT Real Time system (Life Technologies) with primers for vWf (F: CGGCAGCACATACACACATT and R: CGTTCCATCCACAGAGAGGT) and two housekeeping genes (eif1a F: GAGAAGTTCGAGAAGGAAGC and R: CGTAGTATTTGCTGGTCTCG, and b-actin F: CGAGCAGGAGATGGGAACC and R: CAACGGAAACGCTCATTGC). The ΔΔCt method was used for data analysis.

38Reads from RNA-seq were aligned to the zebrafish genome (Zv9.77) combined with sequences for eGFP and ERCC spike-ins as artificial chromosomes, using STAR (version 2.3;. The Ensembl Genes annotation track from UCSC was used with the read_distribution.py tool from the RSeQC tool suite to generate quality control information. Gene expression was quantified using the Salmon reads mode of Sailfish (; parameter -l IU) using Zv9 cDNA sequences from Ensembl version 77 as transcript sequences, together with ERCC spike-in and eGFP sequences as artificial transcripts. Based on comparison with empty control wells, samples with less than 50,000 paired reads and 1,000 expressed genes were considered unfit and were excluded from further analysis (Figure S2).

39For the follow-up experiment, expression was quantified the same way. We used a different stock and concentration of ERCC spike-ins, which changed the scales of the QC values. For these samples, we excluded cells with less than 200,000 paired reads and less than 150 expressed genes (Figure S6).

40Downstream analysis was performed using Transcripts per million (TPM) values reported by Salmon. The TPM unit is a measure of relative abundance of a gene, which is stable across samples. Before analysis expression for endogenous spike-ins were filtered out for each cell, and the TPM for each cell was rescaled to sum to a million. This gives us the interpretation that TPM of a gene will correspond to the concentration of mRNAs from a gene in a given cell.

41Unless stated otherwise, for all analyses, we filtered out genes expressed at a level higher than 1 TPM in only less than three cells, which leaves 20,556 genes.

42We used ICA to identify four latent factors (hidden variables modeling the data), as implemented in scikit-learn (with parameter random_state = 3,984 for the sake of reproducibility). The choice of four components was based on testing between one and ten components, and seeing diminishing returns on the Frobenius norm reconstruction error past four components. One latent factor explains a progression among EGFPlow cells; another factor explains a switch from EGFPlow cells toward the population of EGFPhigh cells. A third factor explains progression among EGFPhigh cells. The fourth factor identifies three outlier cells. We used the fluorescence levels of GFP to flip the orientation of the latent factors so that a higher factor value always corresponded to a higher GFP value. Because these factors are orthogonal, they are statistically independent. In other words, there are three distinct processes happening sequentially. We performed hierarchical Ward clustering of the cells in the four-dimensional ICA space, and assigned the cells to six clusters. (For exact commands, see Notebook 1 in Data S2.) Based on which cluster the cells belonged to, and which factor explains the variability of the cells of that cluster, we ordered cells along this three-stage progression. This ranking of cells through the entire process was treated as pseudotime. (For exact commands, see Notebook 3 in Data S2.)

43As an alternative way to estimate a pseudotime, we applied a Bayesian Gaussian process latent variable model with a one-dimensional latent variable. Briefly, the Bayesian GPLVM will infer a nonlinear function from an unobserved latent space to an observed high-dimensional space, using inducing inputs that are variationally inferred, which helps smooth the data and speed up computation. In our case, the latent space is the one-dimensional pseudotime, and the non-linear function will be a mapping from pseudotime to gene expression values. We used the BayesianGPLVM implementation in the GPy package using a Radial Basis Function (RBF) kernel on the log-transformed TPM values, all other parameters default. Without any information about the EGFP expression, the BayesianGPLVM recovers our original ordering, up to orientation (Spearman correlation 0.97; Figure 2B) (Notebook 7 in Data S2).

44To depict the structure of the data in a friendly way, we performed t-distributed stochastic neighbor embedding (t-SNE) of the four latent factors into two dimensions. The goal of the t-SNE algorithm is to attempt to preserve both global and local structures of higher dimensional data in two dimensions. It additionally tries to not crowd areas with too many points, making them hard to see. We set the perplexity parameter to 75 and used a fixed random seed to make sure the t-SNE plot would be reproducible (parameter random_state = 254 in the scikit-learn implementation of t-SNE).

45We can depict the inferred pseudotime by regressing it into the two-dimensional tSNE space (Figure 2A) and can see how well the two methods of constructing pseudotime agrees.

46To discover marker genes for the clusters of cells, we trained a random forest model for each cluster versus the rest of the cells. We used the Gini feature importance scores for each gene to order genes by how well they can distinguish a cluster from the rest of the cells. We used the ExtraTreesClassifier implementation in the Python machine learning package scikit-learn, with the parameter n_estimators = 100,000. (For the exact commands, see Notebook 2 in Data S2.)

47We treated the pseudotime progression order of the cells as a time series, and for each gene trained two Gaussian processes (GPs): one with a radial basis function (RBF) kernel (which can model change over time) and one with a constant kernel (which assumes that the expression of a gene does not change over time). After optimizing parameters for both models, we filtered the genes by the ratios of the likelihoods of the models. If the RBF kernel GP has a higher likelihood than the constant kernel GP, we can conclude that the gene in question has expression that is dynamic in time. Once we had identified genes that were pseudotime-dependent, we applied the mixtures of hierarchical Gaussian processes model to identify groupings of genes with similar pseudotime expression patterns. All functional enrichment analysis was performed with the gProfiler web service with the standard gene list as background (see Notebook 4 in Data S2 for exact commands).

48We obtained the list of duplicated genes arising from the teleost-specific genome duplication event from. We filtered the list to only retain pairs of genes whose IDs were present in version 77 of Ensembl. For these genes, we binarized the expression to “expressed” or “not expressed” in each cell based on whether the TPM was greater than 1. Using these binary values, for each Ohnolog pair we counted cells expressing either member of the pair, both members of the pair, or none of the members in the pair. Ohnolog pairs in which none of the members were expressed in more than 300 cells and were annotated as “Not expressed.” We defined a value “both_min_diff” as the difference between the smallest number of cells expressing only one of the members in a pair, and the number of cells expressing both members of the pair. Ohnolog pairs with a “both_min_diff”-value larger than 15 were annotated as “XOR Ohnologs.” To identify Ohnolog pairs in which only one member was used, we looked at the difference between the largest number of cells using one member compared to the largest number of cells using the other member. If this difference was larger than 60 cells, the Ohnolog pair was considered a “Single Ohnolog.” The remaining cells were dubbed “Mixed Ohnologs,” meaning cells with a mixture of both members of a pair. (See Notebook 5 in Data S2 for exact commands.)

49All analysis scripts are provided as IPython notebooks in the supplemental information (Data S1, Sample Information) together with a table of detailed information of each sample (Data S2, Analysis Files).

50I.C.M, C.L., and L.F. performed experiments. V.S. carried out the analysis. I.C.M, C.L., V.S., F.H., S.T., and A.C. contributed to the discussion of the results. T.V. oversaw implementation of the scRNA-seq pipeline. I.C.M, C.L., V.S., and A.C. wrote the manuscript. S.H. edited the manuscript. A.C. conceived the study. All authors approved the final version of the manuscript.