Single cell RNA-Seq reveals distinct stem cell populations that drive sensory hair cell regeneration in response to loss of Fgf and Notch signaling

Loss of sensory hair cells leads to deafness and balance deficiencies. In contrast to mammalian hair cells, zebrafish ear and lateral line hair cells regenerate from poorly characterized, proliferating support cells. Equally ill-defined is the gene regulatory network underlying the progression of support cells to cycling hair cell progenitors and differentiated hair cells. We used single cell RNA-Sequencing (scRNA-Seq) of lateral line sensory organs and uncovered five different support cell types, including quiescent and activated stem cells. In silico ordering of support cells along a developmental trajectory identified cells that self-renew and new groups of genes required for hair cell differentiation. scRNA-Seq analyses of fgf3 mutants, in which hair cell regeneration is increased, demonstrates that Fgf and Notch signaling inhibit proliferation of support cells in parallel by inhibiting Wnt signaling. Our scRNA-Seq analyses set the foundation for mechanistic studies of sensory organ regeneration and is crucial for identifying factors to trigger hair cell production in mammals. As a resource, we implemented a shiny application that allows the community to interrogate cell type specific expression of genes of interest.


Introduction
Non-mammalian vertebrates readily regenerate sensory hair cells during homeostasis and after injury, whereas in mammals hair cell loss leads to permanent hearing and vestibular loss (Bermingham-McDonogh and Rubel, 2003;Brignull et al., 2009;Corwin and Cotanche, 1988;Cruz et al., 2015;Ryals and Rubel, 1988). The molecular basis for the inability of mammals to trigger proliferation and a regenerative response is still unknown. Understanding hair cell production in regenerating species is essential for elucidating how regeneration is blocked in mammals. We and others showed that the zebrafish lateral line system is a powerful in vivo model to study the cellular and molecular basis of hair cell regeneration (Kniss et al., 2016;Lush and Piotrowski, 2014b;Ma and Raible, 2009;Romero-Carvajal et al., 2015;Viader-Llargues et al., 2018). The lateral line is a sensory system that allows aquatic vertebrates to orient themselves by detecting water motion. The lateral line organs (neuromasts), distributed on the head and along the body contain approximately 60 cells, composed of central sensory hair cells surrounded by support cells and an outer ring of mantle cells (Fig. 1A, B). The lateral line system is one of the few sensory organs where stem cell behaviors can be observed at the single cell level in vivo because the organs are located in the skin of the animal, are experimentally accessible and easy to image. These properties make it an ideal system to study hair cell regeneration. Despite the unusual location of the hair cells on the trunk, lateral line and ear hair cells develop by similar developmental mechanisms. Importantly, the morphology and function of sensory hair cells are evolutionarily conserved from fish to mammals (Duncan and Fritzsch, 2012;Nicolson, 2005;Whitfield, 2002). For example, mutations in genes causing deafness in humans also disrupt hair cells in the zebrafish lateral line (Nicolson, 2005, Fig. S2, Data file 6). We therefore hypothesize that the basic gene regulatory network required for hair cell regeneration could be very similar in zebrafish and mammals. In support of this hypothesis, our findings that Notch signaling needs to be downregulated to activate Wnt-induced proliferation during hair cell regeneration is also true in the mouse cochlea (Li et al., 2015b;Romero-Carvajal et al., 2015).
In zebrafish regeneration occurs via support cell proliferation and differentiation, and in chicken and amphibians hair cells regenerate from proliferating and transdifferentiating support cells (Balak et al., 1990;Bermingham-McDonogh and Rubel, 2003;Harris et al., 2003;Jones and Corwin, 1996;Lush and Piotrowski, 2014b;Ma et al., 2008;Wibowo et al., 2011b;Williams and Holder, 2000). Yet, even in regenerating species, support cells are not well characterized due to a dearth of molecular markers and the lack of distinct cytological characteristics.
Using time-lapse analyses and tracking of all dividing cells in regenerating neuromasts, coupled with cell fate analyses, we previously identified two major support cell lineages: 1) support cells that divide symmetrically to form two progenitor cells (amplifying divisions); and 2) support cells that divide to form two hair cells (differentiating divisions) (Romero-Carvajal et al., 2015). A recent publication has confirmed these lineages (Viader-Llargues et al., 2018). The two cell behaviors display a striking spatial compartmentalization. Amplifying divisions occur in the dorsal-ventral (D/V) poles and differentiating divisions occur in the center. Mantle cells surrounding the support cells only divide after severe injury to the neuromasts but rarely divide if only hair cells are killed. In addition, we observed quiescent support cells. Thus, there are at least 4 support cell types in a neuromast that likely play different functions in balancing progenitor maintenance with differentiation to ensure the life-long ability to regenerate.
To identify the gene regulatory network that triggers support cell proliferation and hair cell differentiation, we previously performed bulk RNA-Seq of support cells at different time points during regeneration (Jiang et al., 2014). These studies revealed the dynamic changes in signaling pathway activations over time. However, the unexpected diversity and mosaicism in support cells that we discovered during fate analyses (Romero-Carvajal et al., 2015) was masked in the RNA-Seq analysis of pooled support cells. To determine how many support cell populations exist in a neuromast and to characterize their transcriptomes, we performed scRNA-Seq analysis on 1521 purified, homeostatic neuromast cells from a transgenic line. As the lateral line neuromasts are GFP-positive and only consist of about 60 cells we were able to purify a rich cohort of lateral line cells (25x coverage). Our analysis identified seven major neuromast cell populations, revealing genes that are specifically expressed in these cells and characterized the transcriptional dynamics of the process of differentiation and of progenitor maintenance. These results led to the hypothesis that some support cell populations are involved in signaling to trigger regeneration, which we tested by scRNA-Seq analyses of fgf3 mutants that strikingly show increased proliferation and hair cell regeneration. Our scRNA-Seq analysis identified fgf3 targets that we could not identify in bulk RNA-Seq analyses. Importantly, we show that Notch and Fgf signaling act in parallel and that both need to be downregulated together to induce efficient regeneration. Knowing the temporal dynamics and identity of genes required for proliferation and hair cell differentiation are essential for devising strategies to induce hair cell regeneration in mammals.

Single cell RNA-Seq reveals support cell heterogeneity
We reasoned that transcriptional profiling of homeostatic neuromast cells would identify known and previously uncharacterized support cell populations. In addition, as hair cells are continuously replaced, we aimed to identify amplifying and differentiating support cells at different stages of differentiation. We isolated neuromast cells by fluorescence activated cell sorting (FACS) from 5 day post-fertilization (dpf) dissociated transgenic zebrafish in which hair cells, as well as support cells are GFP-positive (Et(Gw57a);Tg(pou4f3); Fig. 1A) and performed scRNA-Seq analyses using the 10X Chromium System (Data file 1). The lateral line also possesses neuromasts with an epithelial planar cell polarity that is offset by 90 o depending on which primordium they originated from (primI or primII, (López-Schier et al., 2004)) and the scRNA-Seq analysis contains cells from all of these neuromasts. For clarity we only discuss gene expression patterns in primI-derived neuromasts.
Unsupervised clustering partitioned 1,521 neuromast cells into 14 different clusters (Butler et al., 2018). We combined some of the less well-defined clusters and identified 7 major neuromast cell types (Fig. 1C, D). For each population we identified genes specifically expressed or highly enriched (Fig. 1C, E; Data file 2). The t-distributed Stochastic Neighbor Embedding (t-SNE) plots for Table 1E are shown in Data file 3.
Based on marker gene expression, groups 1, 2 and 4 encompass the hair cell lineage, with cluster 1 being differentiated hair cells, cluster 2 being young hair cells, and 4 representing proliferating hair cell progenitors (Fig. 1D). The other groups of cycling cells belong to clusters 3 and 10, and because they fail to express hair cell lineage genes, they likely represent amplifying support cells. The proposed hair cell and amplifying support shows that delta ligands are only expressed in a subset of the young hair cells (light green). Lfng and ebf3a mark the most basal, central support cells (Fig. 1I, J, S, U; blue).
Lfng is also expressed in support cells that are situated underneath hair cells in the mouse cochlea (Maass et al., 2016). The central cell population in neuromasts expresses gata2a/b and slc1a3a/glasta, which are markers for hematopoietic and neural stem cells, respectively ( Fig.3 I; Data file 6, (Hewitt et al., 2016;Llorens-Bobadilla et al., 2015)).
Within the central cell cluster, a subset of cells expresses other stem cell-associated genes, such as isl1 and fabp7a (clusters 7, 9; Fig. 1K; (Kim et al., 2016;Makarev and Gorivodsky, 2014;Morihiro et al., 2013;Shin et al., 2007)). In addition, members of the retinoic acid pathway, such as crabp2a, dhrs3a and rdh10a are restricted to clusters 7 and 9 (Fig. 1E). Even though central cells express genes characteristic for stem cells in other systems, our lineage tracing experiments showed that they only give rise to hair cells and do not self-renew (Romero-Carvajal et al., 2015).
Cells in the D/V poles of neuromasts that express wnt2 are located immediately adjacent to the mantle cells, and proliferate to generate more support cells that do not differentiate into hair cells (see below; (Romero-Carvajal et al., 2015)). As these cells selfrenew and possibly represent a stem cell population, we were particularly interested in characterizing new markers for these cells and tested the expression of sost, fsta, srrt/ars, six2b and adcyap1b (Fig. 1E, L, T; orange cells). However, all of these polar genes are expressed in more cells than just the ones immediately adjacent to mantle cells, precluding us from obtaining a specific marker for the amplifying cells ( Fig. 1T; red cells).
Moreover, D/V polar cells do not form their own cluster but are distributed throughout several clusters including some mantle cells (clusters 5,6), central cells (cluster 11) and amplifying, non-differentiating cells (cluster 3).
The most distinct support cell population are the mantle cells represented in clusters 5 and 6 and marked by ovgp1 and sfrp1a (Fig. 1M, N). Mantle cells are the outermost cells in a neuromast and sit immediately adjacent to amplifying support cells give rise to support and hair cells and constitute long term stem cells (Seleit et al., 2017).
In addition, they give rise to postembryonic neuromasts during development and restore neuromasts on regenerating tail tips (Dufourcq et al., 2006;Ghysen and Dambly-Chaudière, 2007;Jones and Corwin, 1993;Ledent, 2002;Wada et al., 2010). In addition to representing stem cells, mantle cells may provide the amplifying support cells with niche factors.
We also identified a number of genes that are expressed in a ring-like pattern, such as fndc7rs4, tfap2a, tcf7l2, hopx, cmah and alpl (Fig. 1O, data not shown). These genes are not restricted to any cluster but are expressed in mantle cells, anterior-posterior (A/P) cells and polar cells. Expression is relatively low in central cells and absent in the hair cell lineage (Data file 3). Interestingly, hopx, cmah and alpl are stem cell markers in different systems raising the possibility that they also mark stem cells in neuromasts (Fuchs, 2009;Takeda et al., 2011). In summary we identified and validated the presence of previously unknown support cell populations, some of which are signaling to trigger regeneration, as shown below.

Cycling cells characterize the amplifying and differentiating lineages
As proliferation is the basis for zebrafish hair cell regeneration we were particularly interested in identifying cycling support cells. Cells in clusters 10, 3 and 4 express pcna (proliferating cell nuclear antigen), required for DNA replication and repair, as well as the mitotic spindle regulator stmn1a (Fig. 1D, 2A, D; Data file 4; (Rubin and Atweh, 2004).
Genes that regulate early versus later phases of the cell cycle are expressed in complementary subsets of the pcna+ cells ( Fig (Cai and Groves, 2015)). As cluster 10 and 3 cells are only defined by the presence of cell cycle genes, we wondered which noncycling support cell type they might be most closely related to. To mitigate the effect of cell cycle genes, we regressed out S and G2/M phase genes using Seurat's cell cycle scoring function (Butler et al., 2018). Using the original cluster classification on the newly generated t-SNE plot, we observed that several cluster 3 (red) cells are now intermingled with D/V support cells in clusters 10, 11, 12 (Fig. S1). Likewise, some of the cluster 4 hair cell progenitor cells now localize within the central support cells (Fig. S1). Thus, cluster 3 cells likely belong to the group of amplifying support cells adjacent to mantle cells that give rise to two undifferentiated daughter cells (Fig. 2E, F, G; support cells, red), whereas cluster 4 cells are more central support cells that give rise to two daughters that differentiate into hair cells (Fig. 2E, F, G; green cells, (Romero-Carvajal et al., 2015)).
Long term stem cells are often relatively quiescent in the absence of a severe or prolonged injury. We labeled 5dpf homeostatic neuromasts for 24 hr with BrdU and subsequently scored non-proliferating support cells (grey squares), and progeny of the dividing cells that differentiated into GFP+ hair cells (green squares) or remained support cells (red squares, Fig. 2G; reanalyzed data from Romero-Carvajal et al., 2015). To visualize and compute the ratio of quiescent cells, we plotted the location of each progeny and calculated the BrdU index and spatial distribution of the different cell types (Romero-Carvajal et al., 2015;Venero Galanternik et al., 2016). Amplifying divisions are restricted to the D/V poles, whereas differentiating divisions are more centrally located but do not show a bias toward any quadrant. D/V poles possess more cells than the A/P poles ( Fig.  2H), however, 6.5% of the D/V cells proliferate, whereas only 0-1% proliferate in the A/P poles (Fig. 2I). Thus, cells in the A/P poles and central cells beneath the hair cells are relatively quiescent during homeostasis (Cruz et al., 2015;Romero-Carvajal et al., 2015).
The expression of the zona pellucida-like domain-containing protein 1 gene si:ch73-261i21.5 in the A/P poles has a striking complementary expression pattern to the D/V maker wnt2 and is expressed in quiescent cells ( Fig. 2J and K). Therefore, A/P genes could play a role in regulating quiescence.
Having established different support and hair cell populations and their proliferation status allows us to interrogate the expression pattern of any gene (see shiny app).

Heatmaps and the ordering of cells along pseudotime reveals how gene networks change in different lineages
Cell cycle analyses suggested the existence of two lineages of cycling cells. To long term stem cells (Seleit et al., 2017).
Heatmaps of factors involved in ribosome and protein synthesis also provide lineage information. Quiescent neural stem cells show low levels of ribosomal subunits and protein synthesis (Llorens-Bobadilla et al., 2015). Likewise, we observed mantle (clusters 5,6) and A/P cells (cluster 13) expressing low levels of rpl (ribosomal protein-L) and rps (ribosomal protein small subunit) genes but these levels significantly increase in clusters 12 and 11 and the dividing cells in clusters 10 and 3 ( Fig. 3B', Fig. S3). Also, the hair cell lineage and the central cells (clusters 7, 8, 9 and 14) show little ribosome synthesis, which drastically increases in dividing hair cell progenitors ( Fig. 3C', Fig. S3).
The low ribosome synthesis levels suggest that central support and mantle cells resemble quiescent stem or progenitor cells.
Mantle cell genes show fairly specific gene expression, such as tspan1 but also share genes with clusters 12, 13, and amplifying support cells in clusters 10, 3, suggesting To visualize the expression dynamics of all detected hair cell lineage genes we generated a heatmap in which cells are ordered along pseudotime on the x-axis (Fig. 4I).
The heatmap reveals groups of genes that possess similar expression dynamics and likely form a regulatory network within each cell cluster (Fig. 4I). This map of progressive gene activation serves as a blueprint for hair cell specification and differentiation.
Hair cell specification and differentiation also depends on the downregulation of genes (Matern et al., 2018). Node 17 (Data file 8) shows such genes that are enriched in support cell types and are downregulated in the hair cell lineage (clusters 4,2 and 1; Fig.   4J). For example, in situ hybridization with fndc7a shows that it labels support cells as in the mouse and that the more apically located young and mature hair cells are unlabeled

scRNA-seq reveals that loss of fgf3 in central support cells leads to increased proliferation and regeneration
Our previous bulk RNA-Seq analysis of regenerating neuromasts revealed that Fgf regeneration, as we showed for Notch signaling (Ma et al., 2008;Romero-Carvajal et al., 2015). Indeed, hair cell regeneration is enhanced in fgf3 mutant larvae and even during homeostasis the total cell number per neuromast is increased ( Fig. 5I-L).
Because fgf3 disappears as hair cells die, we wondered if fgf3 is expressed in hair cells. The scRNA-Seq analysis of homeostatic neuromasts shows that ligand and receptor expression is complex, and that Fgf signaling is not active in young or mature hair cells ( Fig. 5M; clusters 2,1). fgf3 is expressed exclusively in central support cells (clusters 7,8,9) and is downregulated in response to death of the overlying hair cells ( To identify genes/pathways underlying the increased proliferation in fgf3 -/-, we first performed bulk RNA-Seq analysis with 5 dpf homeostatic fgf3 -/and siblings. However, the differences in gene expression between siblings and mutants was too low. We therefore performed scRNA-Seq analyses on 1,459 fgf3 mutant and 1,932 sibling cells (Data file 13). A t-SNE plot in which we plotted mutant and sibling datasets together demonstrates that the variance between the two datasets is small as the two datasets intermingle (Fig. 7A, sibling blue, mutant red). The plot also shows that mutant cells contribute to each cluster and that therefore no major cell type is missing in fgf3 -/- (Fig.  7A, B, 1D).
However, the scRNA-Seq analysis revealed fgf3 targets that are down-or upregulated in the mutants (Fig. 7C, Fig. S9). We were particularly interested in genes that regulate the Wnt pathway and/or proliferation and identified that some of the D/V polar genes, such as sost and adcyap1b are downregulated in the mutants (Fig. 7C, Fig.   S9C). We validated the downregulation of Wnt inhibitor sost by in situ hybridization in fgf3 -/and dn:fgfr1 larvae ( Fig. 7D-E'). Interestingly, sost is also downregulated 1 hour after hair cell death in wildtype larvae, suggesting that the downregulation of Fgf signaling after neomycin could also be responsible for the downregulation of the Wnt inhibitor sost ( We conclude from these data that Notch and Fgf signaling largely act in parallel to inhibit Wnt signaling, with a small amount of proliferation being inhibited by Notch via the upregulation of Fgf signaling (Fig. S10). Thus Notch and Fgf need both to be independently and transiently downregulated for efficient hair cell regeneration.

Discussion
The

Cells can transition from an amplifying to differentiating lineage
Interestingly, the PAGA graph (Fig. 3M) shows a connection between cluster 3 (the amplifying support cells) and cluster 4 (the proliferating differentiating support cells) suggesting that the amplifying support cells can turn on differentiation genes as they are displaced toward the center of the neuromast and also become hair cells. This hypothesis is supported by time lapse movies of amplifying support cells during homeostasis and regeneration (Romero-Carvajal et al., 2015). Our finding shows that the amplifying and hair cell lineages are not predetermined, but rather that the location of cells within the neuromast and the signals they are exposed to determine their fate.

Support cells share gene expression profiles with stem cells in other organs
The notion that some support cells constitute stem cells is also supported by the finding that their gene expression profiles share many genes with stem cells in the CNS, heart, intestine or hair follicles. other organs, such as fabp7/Blbp, which labels glioma stem cells, radial glia cells, NSCs, (Kim et al., 2016;Morihiro et al., 2013) and isl1 that is expressed in quiescent intestinal stem cells and stem cells in the heart (Makarev and Gorivodsky, 2014;Shin et al., 2015).
However, in lateral line neuromasts central cells only give rise to hair cells and we therefore do not consider them to be stem cells.
A comparison of the transcriptional profiles of support cells in regenerating species, such as zebrafish and chicken with mouse support cells will be highly informative. Do mammalian support cells also still express many of the above-mentioned stem cell genes, or do mammalian support cells represent a more differentiated cell population? The results from such analyses will help us determine if mammalian support cells need to be reprogrammed for efficient induction of regeneration.

Signals that control multipotency versus differentiation
It is still unknown what signals prevent amplifying support cells from differentiating in the D/V poles, as Notch signaling does not seem to be active in these regions based on Notch reporter expression (Tg(Tp1bglob:eGFP, (Romero-Carvajal et al., 2015)).
Likewise, loss of Fgf3 signaling does not lead to the loss of amplifying divisions in the D/V poles demonstrating that Fgf3 signaling is not involved in preventing amplifying support cells from differentiating. D/V polar gene expression, such as fsta, wnt2 and sost, correlates with the occurrence of amplifying support cells and these genes are downregulated as cells differentiate into hair cells. However, their expression extends further centrally beyond the row of cells immediately next to the mantle cells suggesting that they are not specifically regulating amplifying cells. We therefore hypothesize that, as immediate neighbors of amplifying support cells, mantle cells are involved in sending anti-differentiation signals. These signals still remain to be identified.

Fgf3 inhibits Wnt signaling and proliferation possibly via Sost
Wnt signaling induces and is required for support cell proliferation in neuromasts sprouty2 is also down regulated in fgf3 -/after hair cell death and could perform a similar inhibitory function in homeostatic neuromasts (Data file 14). Fgf signaling also inhibits proliferation in the regenerating utricle and the basilar papilla of chicken (Jiang et al., 2018;Ku et al., 2014). Likewise, addition of Fgf2 or Fgf20 to auditory or vestibular cultures inhibits support cell proliferation (Ku et al., 2014;Oesterle et al., 2000), suggesting that the inhibitory effect of Fgf signaling on progenitor proliferation is evolutionary conserved in species that can regenerate their sensory hair cells.
Experiments utilizing chemical inhibition of Fgfr have been less clear and in some studies had no effect on proliferation or even led to the inhibition of proliferation (Jacques et al., 2012, Ku et al., 2014, Lee et al., 2005. The differences in the effect of chemical inhibition can likely be attributed to differences in culture conditions, timing or doses of drug treatment, underscoring the importance of utilizing multiple methods of signaling pathway inhibition, particularly gene mutations.

A role for Fgf is specification?
While our data suggest that the main role of fgf3 in mature neuromasts is in regulating proliferation, it is also possible that fgf3 acts to maintain hair cell progenitors in an undifferentiated, non-sensory state as has been observed in the developing zebrafish ear (Maier and Whitfield, 2014). Likewise, in the mouse cochlea loss of Fgfr3 results in increased hair cells at the expense of support cells (Hayashi et al., 2007;Puligilla et al., 2007), while activating mutations in Fgfr3 or loss of function mutation in the Fgfr inhibitor Spry2 lead to transformation of one support cell type into another (Mansour et al., 2013;Mansour et al., 2009;Shim et al., 2005). However, a possible inhibitory effect of fgf3 on hair cell fate has to be rather subtle, as we only observe a limited increase in hair cell numbers, and our scRNA-Seq analyses did not reveal any obvious candidates that might regulate cell fate in a Fgf-dependent fashion.

Interactions between the Fgf and Notch pathways
Just like Fgf, Notch signaling is immediately downregulated after hair cell death causing cell proliferation. During homeostasis both pathways inhibit proliferation through negative regulation of Wnt activity (Ma et al., 2008;Romero-Carvajal et al., 2015). To test if Fgf and Notch act in the same pathway we performed epistasis experiments. In fgf3 -/mutants expression of a Notch reporter and Notch target genes are not affected suggesting that Notch signaling is largely intact in fgf3 -/- (Fig. 8G-K'). Also, while pharmacological inhibition of Notch activity during homeostasis decreases fgf3 expression it has no effect on Fgf target genes ( Fig. 8L-P'). Additionally, activating Notch signaling by over-expression of a Notch intracellular domain inhibits proliferation and differentiation during regeneration in both wildtype and fgf3 -/- (Fig. 8R, T, V). This finding shows that Notch inhibits Wnt signaling even in the absence of Fgf3. These data argue that Fgf and Notch signaling are functioning largely in parallel. However, hs:nicd by itself is more effective in inhibiting proliferation than hs:nicd in fgf3 -/-, arguing that Notch signaling likely inhibits a small fraction of proliferation via Fgf (Fig. 8V, Fig. S10). We therefore postulate that the transient down regulation of both of these pathways immediately after hair cell death is required to maximally activate Wnt signaling and induce proliferation. Nevertheless, the Fgf and Notch pathways are reactivated before hair cell regeneration is complete (Jiang et al., 2014) and play additional roles later in regeneration, such as ensuring that not too many support cells differentiate. As such, in mammals a short inhibition of one or both of these pathways is more likely to result in functional regeneration than prolonged treatments.

Conclusion
The scRNA-Seq analysis revealed previously unidentified support cell populations and combined with in situ validation of these cell groups identified lineages that either lead to stem cell self-renewal or hair cell differentiation. Importantly, we have identified the cascade of gene activation and repression that leads to hair cell differentiation. Our analyses led to the hypothesis that some of the support cell populations are involved in signaling to trigger regeneration, which we tested by scRNA-Seq analyses of fgf3 mutants that strikingly show increased proliferation and hair cell regeneration. These experiments identified fgf3 targets that we could not identify in bulk RNA-Seq analyses. Having characterized the support cell transcriptome of a regenerating species allows us to identify commonalities and differences with mouse support cells that do not trigger a meaningful regenerative response (Burns et al., 2015;Maass et al., 2016). Such a comparison will become even more powerful once adult mouse single cell transcriptomes of support cells are available.
(B) Heatmap of genes upregulated in fgf3 mutant lateral line cells.

Supplementary Material
Data file 1. Related to Figure 1: excel file of genes that are expressed in at least three cells.
Data file 2. Related to Figure 1E: excel file of cluster marker genes.
Data file 3. Related to Figure 1E: t-SNE plots of all cluster marker genes.
Data file 4. Related to Figure 2D: excel file of cell cycle genes.
Data file 5. Related to Figure 2D: t-SNE plots of cell cycle genes.
Data file 6. Related to Figure
Generation of Tg(fgf3:H2B-mturquoise2) psi60Tg H2B-mturquoise2 was placed near the 5' region of fgf3 using non-homologous repair with Crisper/Cas9 (Auer et al., 2014;Kimura et al., 2014). A Crisper recognition site (GGCCATGGAAACTAAATCTGCGG) was chosen 584bp in front of the fgf3 transcription start site. The same recognition site was cloned by PCR onto both ends of a construct containing a 56bp -actin minimal promoter driving human-H2B fused at the c-terminus with mturquoise2 followed by the SV40 polyA from the Tol2 kit (Kwan et al., 2007). This plasmid was mixed with the crRNA (GGCCATGGAAACTAAATCTG, IDT), tracrRna (IDT) and Cas9 protein (PNA Bio) and complexed for 10 minutes at room temperature then placed on ice. The complex was then injected into the cell of a one cell stage zebrafish embryo. Integrated DNA will just contain the minimal -actin promoter, H2B-mturquoise2 and the polyA sequence without any plasmid DNA. In a few larvae, fluorescence could be seen by 24 hrs and onward. Fluorescent embryos were sorted and raised to identify founders that showed H2B-mturquoise expression similar to fgf3 expression. All experiments were performed according to guidelines established by the Stowers Institute IACUC review board.

Sensory hair cell ablation
For hair cell ablation experiments 5 days post fertilization (dpf) fish were exposed to 300mM Neomycin (Sigma-Aldrich) for 30 minutes at room temperature. Neomycin was then washed out and larvae were allowed to recover at 28°C for as long as the experiment lasted.

Heat shock paradigm
Heat shock induction of transgene expression varied depending on the transgenic line used. Initially 5dpf larvae were heat shocked for 1hr at 37°C then put back at 28°C for 1hr. Larvae were then heat shocked for 1hr at a higher temperature (39°C for notch1aintraceullar and 40°C for dkk1a and dnfgfr1a), put at 28°C for 1 hr, followed by another higher temperature heat shock (39 °C or 40 °C) for 1 hr. Larvae were then allowed to recover for 1hr at 28°C then fixed in 4% paraformaldehyde for in situ hybridization or continued with neomycin and proliferation analysis. For neomycin experiments, larvae were treated for 30 minutes in 300µM, then extensively washed and immediately placed in EdU as described above. Larvae were then heat shock at the higher temperature followed by 2hrs recovery at 28°C. The 1hr heat shock followed by 2hrs recovery was repeated for a total of 24hrs post neomycin treatment. Larvae were then fixed in 4% paraformaldehyde overnight at 4°C.

scRNA-Seq read alignment and quantification
Raw reads were demultiplexed and aligned to version 10 of the zebrafish genome (GRCz10) using the Cell Ranger pipeline from 10X Genomics (version 1.3.1 for wildtype and version 2.1.1 for fgf3 -/data sets). 1,666 cell barcodes were obtained for wildtype embryos, 1,932 for fgf3 siblings and 1,459 for fgf3 mutants. These quantities were estimated using Cell Ranger's barcode ranking algorithm, which estimates cell counts by obtaining barcodes that vary within one order of magnitude of the top 1 percent of barcodes by top UMI counts. The resulting barcodes (henceforth referred to as cells) were used to generate a UMI count matrix for downstream analyses. Data deposition: the BAM files and count matrices produced by Cell Ranger have been deposited in the Gene Expression Omnibus (GEO) database, www.ncbi.nlm.nih.gov/geo (accession no. GSE123241).

Quality control, dimensional reduction, and cell classification
Subsequent analyses on UMI count matrix for all three data sets were performed using the R package Seurat (version 2.3.4, (Butler et al., 2018)) following the standard workflow outlined in the pbmc3k example on the Satija lab webpage (https://satijalab.org).
Both fgf3 sibling and mutant data sets were analyzed independently using the same parameters and arguments, and then each data set was merged using the Seurat function MergeSeurat(). Initial gene quality control was performed by filtering out genes expressed in less than 3 cells for the WT data set, and 5 for the fgf3 sibling/mutant data sets. The remaining UMI counts were then log normalized. Gene selection for dimensional reduction was accomplished using the Seurat function FindVariableGenes() with the following arguments for the WT data set: x.low.cutoff = 0.001, x.high.cutoff = 3.0, and y.cutoff = 0.5; and for the fgf3 sibling/mutant data sets: x.low.cutoff = 0.20 and y.cutoff = -0.20. Following gene selection, all log-normalized expression values were then scaled and centered using ScaleData(). For dimensional reduction, we chose to use the first 6 principal components (PCs) for wildtype and the first 19 for the fgf3 sibling/mutant data sets. PCs were chosen according to the PCA elbow plot, which orders PCs from highest to lowest based on the percentage of variance explained by each PC. Thus, each set of PCs chosen showed the highest percentage of variance explained on the elbow plot.
Next, we performed clustering on each set of principal components, and for twodimensional visualization, we performed a second round of dimensional reduction using t-SNE. Cells for all data sets were classified according to their marker gene expression.
Markers for each cell type were identified based on their differential gene expression using the Seurat function FindAllMarkers(). Genes with an adjusted p-value less than 0.001 were retained. For the pairwise comparison between fgf3 mutant and sibling data sets, we used the function FindMarkers() and retained genes with a fold change greater than 0.10, or less than -0.10. Genes differentially expressed between dendrogram nodes were calculated using the function FindAllMarkersNode(), and we kept the top 100 genes with the highest p-values for each node comparison. All three data sets contained nonspecific cell types that contained markers for skin and blood cells which were removed from the final analysis.

Pseudotime Analysis
All data contained within our processed Seurat object for the wildtype data set was converted to the AnnaData format for pseudotime analysis in Scanpy (version 1.2.2, (Wolf et al., 2018)), using the Seurat function convert. We recalculated k-nearest neighbors at k = 15 and chose cluster 14 as our putative "stem cell" population. Pseudotime was calculated using Scanpy's partitioned-based graph abstraction function, paga. To visualize gene expression in pseudotime, cells from clusters 14, 4, 2, and 1 were subsetted into a separate matrix and ordered according to their pseudotime values from least to greatest. Next, we chose genes differentially expressed between clusters, and differentially expressed between selected nodes in our cluster dendrogram. The resulting genes were ordered according to our understanding of neuromast biology and pseudotime expression patterns. The final count matrix was then log normalized and rendered as a heatmap using the python package Seaborn (version 0.8.1).

GO term analysis
The GO term analysis was performed in DAVID (Database for Annotation, Visualization and Integrated Discovery, Huang da et al., 2009).

Shiny Apps and Data repository
The Shiny app web interface generates t-SNE plots, violin plots, co-expression t-SNE plots and heatmaps from the homeostasis scRNA-Seq data. The user can choose the genes to plot and the results can be saved as pdfs. Instructions are provided in the welcome page of the Shiny app. We are also in the process of uploading the data into gene Expression Analysis Resource (gEAR), a website for visualization and comparative analysis of multi-omic data, with an emphasis on hearing research (https://umgear.org).

RNA extraction and cDNA synthesis
Total RNA was extracted from FAC-sorted cells using Trizol (Thermo Fisher Scientific, Waltham, MA. USA), chloroform and isopropanol. During isopropanol precipitation, to enhance the pellet visibility, GlycoBlue coprecipitant (Thermo Fisher Scientific, Waltham, MA. USA) was used. Following that, total RNA was washed with 80% ice-cold Ethanol and resuspended in RNase-free water. The first-strand cDNA synthesis and cDNA amplification were done using SMART-Seq v4 Ultra Low Input RNA kit (Takara Bio USA, Mountain View, CA. USA) following manufacturer's instructions.

Quantitative PCR
Q-PCR was carried out using ABI SYBR Green master mix (Applied Biosystems, Foster City, USA) in the QuantStudio 7 Real-Time PCR System with a 384-Well Block (Applied Biosystems, Foster City, USA). The reaction program consisted of four steps: UDG treatment (50°C for 10 minutes), quantitation (40 cycles of 95°C for 15s, 60°C for 60s) and melting curve analysis (95°C for 15s, 60°C for 60s, 95°C for 15s). The experiment was performed in triplicate. All the signals were normalised to the ef1a expression level. All primer sequences are provided in Data file 15.