scRNA-Seq reveals distinct stem cell populations that drive hair cell regeneration after loss of Fgf and Notch signaling

1Non-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. 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. 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 (Figure 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. For example, mutations in genes causing deafness in humans also disrupt hair cells in the zebrafish lateral line. 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.

2In zebrafish regeneration occurs via support cell proliferation and differentiation, and in chicken and amphibians hair cells regenerate from proliferating and transdifferentiating support cells. 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.

3Using 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). A recent publication has confirmed these lineages. 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 four 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.

4To 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. 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 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.

5We 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); Figure 1A) and performed scRNA-Seq analyses using the 10X Chromium System (Supplementary file 1). The lateral line also possesses neuromasts with an epithelial planar cell polarity and corresponding gene expression pattern that is offset by 90o depending on which primordium they originated from (primI or primII,). The scRNA-Seq analysis contains cells from all of these neuromasts. For clarity we only discuss and illustrate primI-derived neuromasts.

6Unsupervised clustering partitioned 1521 neuromast cells into 14 different clusters. We combined some of the less well-defined clusters and identified seven major neuromast cell types (Figure 1C–D). For each population we identified genes specifically expressed or highly enriched (Figure 1C,E; Supplementary file 2). Dissociating tissues has the caveat that it likely triggers gene expression changes due to loss of adhesion molecules or to a global injury response. We controlled for a global gene expression response by only analyzing genes that are variable between clusters, however cluster-specific gene up- or down regulation can only be controlled for by performing in situ hybridization in intact organs.

7The t-distributed Stochastic Neighbor Embedding (t-SNE) plots for genes listed in Figure 1E are shown in Supplementary file 3. Based on marker gene expression, clusters 1, 2 and 4 encompass the hair cell lineage, with cluster one being differentiated hair cells, cluster two being young hair cells, and four representing proliferating hair cell progenitors (Figure 1D–E). The other 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 cell lineages are described in detail below. Cells in all other clusters represent different support cell populations (all data can be queried at https://piotrowskilab.shinyapps.io/neuromast_homeostasis_scrnaseq_2018/; see Materials and methods).

8To determine if the distinct cell clusters defined by scRNA-Seq can be detected in neuromasts, we performed in situ hybridization experiments with cluster marker genes (Figure 1E–Q). Mature, differentiated hair cells are centrally and apically located in a neuromast (tekt3, Figure 1F,R,S; dark green). Immediately above the mature hair cells are young hair cells that form a ring and express atoh1b (cluster 2, Figure 1G,R,S). Figure 1H 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 (Figure 1I,J,S,U; blue). lfng is also expressed in support cells that are situated underneath hair cells in the mouse cochlea. The central cell population in neuromasts expresses gata2a/b and slc1a3a/glasta, which are markers for hematopoietic and neural stem cells, respectively (Figure 3I; Supplementary file 6,; ). Within the central cell cluster, a subset of cells expresses other stem cell-associated genes, such as isl1 and fabp7a (clusters 7, 9; Figure 1K;;;; ). In addition, members of the retinoic acid pathway, such as crabp2a, dhrs3a and rdh10a are restricted to clusters 7 and 9 (Figure 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.

9Cells 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; ). As these cells self-renew 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 (Figure 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 (Figure 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).

10The most distinct support cell population are the mantle cells represented in clusters 5 and 6 and marked by ovgp1 and sfrp1a (Figure 1M–N). Mantle cells are the outermost cells in a neuromast and sit immediately adjacent to amplifying support cells (Figures 1B,M,N,R–U). Lineage tracing of mantle cells in medaka revealed that mantle cells give rise to support and hair cells and constitute long term stem cells. In addition, they give rise to postembryonic neuromasts during development and restore neuromasts on regenerating tail tips. In addition to representing stem cells, mantle cells may provide the amplifying support cells with niche factors.

11We also identified a number of genes that are expressed in a ring-like pattern, such as fndc7rs4, tfap2a, tcf7l2, hopx, cmah and alpl (Figure 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 (Supplementary 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. 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.

12As 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 (Figures 1D, 2A and D; Supplementary file 4; ). Genes that regulate early versus later phases of the cell cycle are expressed in complementary subsets of the pcna +cells (Figure 2B,C, Supplementary file 5). Genes with the associated GO terms DNA replication and DNA repair, such as mcm4 are expressed in cells closer to clusters 11 and 14, whereas genes involved in chromosome segregation and mitosis are expressed in cells closer to the younger hair cells in the t-SNE plot (cdk1; Figure 2B,C and D).

13Differential gene expression analysis between the proliferating cells in clusters 3, 4 and 10 revealed that only cluster 4 cells express genes characteristic for the hair cell lineage such as atoh1b and dld (Figure 1G–H; ). As cluster 10 and 3 cells are only defined by the presence of cell cycle genes, we wondered which non-cycling 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. 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 (Figure 2—figure supplement 1). Likewise, some of the cluster four hair cell progenitor cells now localize within the central support cells (Figure 2—figure supplement 1). Thus, cluster 3 cells likely belong to the cluster of amplifying support cells adjacent to mantle cells that give rise to two undifferentiated daughter cells (Figure 2E–G; support cells, red), whereas cluster 4 cells are more central support cells that give rise to two daughters that differentiate into hair cells (Figure 2E–G; green cells, ).

14Long term stem cells are often relatively quiescent in the absence of a severe or prolonged injury. We labeled 5dpf homeostatic neuromasts for 24 hrs 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, Figure 2G; reanalyzed data from ). 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. 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 (Figure 2H), however, 6.5% of the D/V cells proliferate, whereas only 0–1% proliferate in the A/P poles (Figure 2I). Thus, cells in the A/P poles and central cells beneath the hair cells are relatively quiescent during homeostasis. The expression of the zona pellucida-like domain-containing protein one 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 (Figure 2J–K). Therefore, A/P genes could play a role in regulating quiescence.

15Having established different support and hair cell populations and their proliferation status allows us to interrogate the expression pattern of any gene (see web app: https://piotrowskilab.shinyapps.io/neuromast_homeostasis_scrnaseq_2018/ ). For example, a heatmap of human deafness genes demonstrates which genes are expressed in hair cells and which ones are expressed in support cells, crucial information for elucidating their function (Figure 2—figure supplement 2; Supplementary file 6).

16Cell cycle analyses suggested the existence of two lineages of cycling cells that possibly represent the amplifying and differentiating support cells identified in time-lapse analyses. To identify possible lineage progression and relationships between non-cycling, cycling and differentiating cells, we generated a dendrogram of the 14 cell clusters using hierarchical clustering (Figure 3A and Figure 1D, see Materials and methods) with nodes and terminal branches represented by a number. We then produced heatmaps of genes enriched in each of the nodes and branches (Supplementary file 7, 8). Node 15 distinguishes the transcriptome of differentiated hair cells from support cells. In other species, SoxB1 genes characterize the prosensory domain from which hair cells and support cells arise. Node 17 shows that the SoxB1 family member sox2, its target sox21a and sox3 are expressed in support and mantle cells but are downregulated in differentiated hair cells (Figure 3D; Supplementary file 8, Supplementary file 3). Thus, mantle and all support cells constitute the prosensory domain in a lateral line neuromast.

17Genes that span multiple support cell clusters identified two lineages emerging within the prosensory domain (Figure 3B–C, Supplementary file 8). Clusters 7, 9, 8, 14 and 4 comprise the differentiating hair cell lineage, whereas clusters 5, 6, 12, 11, 10 and 3 encompass the amplifying lineage (Figure 3B–C). Indeed, lineage tracing experiments determined that central support cells give rise to hair cells, whereas mantle cells give rise to support and hair cells if traced for several months and constitute long term stem cells.

18Heatmaps 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. 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 (Figure 3B’ and Figure 3—figure supplement 1). 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 (Figure 3C’, Figure 3—figure supplement 1). The low ribosome synthesis levels suggest that central support and mantle cells resemble quiescent stem or progenitor cells.

19Mantle 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 a lineage relationship (amplifying lineage; Figure 3B,E; Supplementary file 8, nodes 18, 23). In addition, clusters 6, 12, 11, 10 and 3 express the D/V polar genes sost, wnt2, adcyap1b and fsta (Figure 3B,F; Supplementary file 8, node 24). These genes label the D/V compartments of neuromasts in which amplifying support cells reside next to mantle cells (Figure 2E–F; ). However, when sost+ cells are displaced towards the center of the neuromasts, they downregulate sost and differentiate into hair cells (lineage from clusters 11/14 to 4, green arrow in Figure 3B and G). The amplifying lineage is also supported by genes such as the pluripotency marker hopx, expressed in mantle and A/P cells as well as proliferating, non-differentiating cells in clusters 10 and 3 (Figures 1E,O and 3B,G; Supplementary file 8, nodes 18, 23, 3; ). We conclude that the cells adjacent to mantle cells (and possibly mantle cells themselves) constitute the amplifying lineage.

20The differentiating lineage is marked by atoh1a which specifies hair cells together with the downstream delta ligands (Figures 3H and 4A). atoh1a is expressed in non-cycling (clusters 7, 8, 9 and 14) and cycling hair cell progenitors (cluster 4), as well as young hair cells (cluster 2). These cells belong to a subset of central support cells that are marked by progenitor markers slc1a3a/glasta, lfng, ebf3a, gata2b and prox1a (Figure 1I–J and Figure 3I–K, Supplementary file 8, nodes 26, 22 and 7). Indeed, time lapse analyses of regenerating neuromasts in a prox1a reporter line show that central cells downregulate prox1a as they divide, while turning on the hair cell marker pou4f3:gfp (Figure 3—figure supplement 2, Figure 3-video 1), as described for the mouse cochlea. We conclude that central cells (clusters 7, 9, 8,14 and 4) contain or represent hair cell progenitors.

21After establishing that hair cell progenitors reside within central cells and that cluster one represents differentiated hair cells, we ordered these cells along a developmental trajectory (henceforth referred to as pseudotime). To define this trajectory, we generated a graph highlighting the relationships between each cluster using partition-based graph abstraction (PAGA,; ). Each node represents a cluster and the weight of the lines represents the statistical measure of connectivity between the nodes. The PAGA model is built on a neighborhood graph of single cells and represents a simplified structure from which a path of differentiation can be inferred. We selected clusters 14, 4, 2,1 to represent a putative path of hair cell differentiation. Cluster 14 cells represent the ‘root’ cell population and cluster one the terminal population, with clusters 4 and 2 representing intermediate populations. Cells in each of these four populations are ordered according to their distance to cluster 14, which is the basis for cell positions in our pseudotime heatmap (Figure 4I). It is important to note that amplifying support cells in clusters 3, 10, and 11 also share a relationship to differentiating support cells. This suggests that amplifying support cells may differentiate if displaced to the center of the neuromast (Figure 3M, arrow).

22Within the hair cell lineage, gene expression changes progressively from the non-cycling progenitors (cluster 14) to differentiated hair cells (cluster 1), reflecting developmental time (Figure 4A–H, Supplementary file 9 and 10). As progenitors are exiting the cell cycle they turn on differentiation genes, many of which are shared between clusters 2 and 1 (Supplementary file 3 and 8, node 16). However, hair cells are subdivided into younger hair cells (cluster 2) and mature hair cells based on differential gene expression (cluster 1; Supplementary file 3, nodes 1 and 2), but the younger hair cells in cluster two possess cilia and can be killed with neomycin (Figure 1R–S; Figure 4—figure supplement 1).

23To 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 (Figure 4I). This heatmap reveals clusters of genes that possess similar expression dynamics and likely form a regulatory network within each cell cluster (Figure 4I). This map of progressive gene activation serves as a blueprint for hair cell specification and differentiation.

24Hair cell specification and differentiation also depends on the downregulation of genes. Node 17 (Supplementary 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; Figure 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 (Figure 4K; ). Likewise, a Tg(NFκB:EGFP) reporter line shows that the NFκB pathway that regulates proliferation and self-renewal in other systems is expressed in support cells but that hair cells are GFP-negative (Figure 4L–N;;; ). A GO term analysis reveals that genes associated with regulation of transcription, translation, protein folding, cell cycle and Wnt signaling are downregulated in clusters 1 and 2 (Figure 4—figure supplement 2). Another group of genes is downregulated as hair cells develop from young to mature hair cells (Figure 4O, Supplementary file 8, node 2 and 24). These downregulated genes are associated with the GO terms translation, regulation of cell cycle and Notch signaling (dla, dlb, dlc, dld). Notch signaling plays an essential role in specifying hair cells versus support cells and a detailed expression analysis of Notch receptors, ligands and downstream targets is shown in Figure 4—figure supplement 3. In sum, the pseudotime heatmap provides a blueprint for the succession of gene activation and repression that occurs in support cells as they differentiate into hair cells, thus providing a framework for experimentally inducing hair cell differentiation in mammals.

25Our previous bulk RNA-Seq analysis of regenerating neuromasts revealed that Fgf pathway genes are downregulated 1 hr after hair cell death ((; Figure 5A–H’) suggesting that the downregulation of Fgf signaling could be involved in triggering regeneration, as was shown for Notch signaling. Indeed, hair cell regeneration is enhanced in fgf3 mutant larvae and even during homeostasis the total cell number per neuromast is increased (Figure 5I–L).

26Because 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 (Figure 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 (Figure 5A’). We generated a fgf3 knock-in line that recapitulates fgf3 expression (Figure 5N–O, Figure 5-video 1) and a cross with a prox1a reporter line shows that fgf3 and prox1a are co-expressed, as predicted by the scRNA-Seq data (Figure 5P–Q). To understand the cellular basis of the increased regeneration response, we performed BrdU analyses of homeostatic and regenerating fgf3-deficient neuromasts. During both homeostasis and regeneration, proliferation is significantly increased in fgf3-/- neuromasts (Figure 6A–F). Downregulation of Fgf signaling by expression of dn:fgfr1a during homeostasis also increases proliferation and neuromasts are significantly bigger in fgfr1a/fgfr2 double mutants, similarly to fgf3-/- (Figure 6—figure supplement 1, Figure 6G–L, Figure 5L). Therefore, fgf3 likely signals via fgfr1a and fgfr2 receptors. The BrdU plots and the rose diagrams show that in fgf3-/- homeostatic and regenerating neuromasts amplifying divisions are not restricted to the D/V poles (Figure 6B–B’ and D–D’). This pattern of proliferating cells looks almost identical to the ones observed after upregulation of Wnt signaling.

27To identify genes/pathways underlying the increased proliferation in fgf3-/-, we first performed bulk RNA-Seq analysis with five 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 1459 fgf3 mutant and 1932 sibling cells (Supplementary 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 (Figure 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-/- (Figures 7A–B and 1D).

28However, the scRNA-Seq analysis revealed fgf3 targets that are down- or upregulated in the mutants (Figure 7C, Figure 7—figure supplement 1). 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 (Figure 7C, Figure 7—figure supplement 1A). We validated the downregulation of Wnt inhibitor sost by in situ hybridization in fgf3-/- and dn:fgfr1 larvae (Figure 7D–E’). Interestingly, sost is also downregulated 1 hr after hair cell death in wild type larvae, suggesting that the downregulation of Fgf signaling after neomycin could also be responsible for the downregulation of the Wnt inhibitor sost (Figure 7F–F’). Also, the Wnt target gene wnt10a is upregulated, illustrating that Wnt signaling is increased in fgf3 mutants (Figure 7G–G’, Figure 7—figure supplement 1A).

29To interrogate if the increase in proliferation in Fgf pathway mutants is due to the upregulation of Wnt signaling, we abrogated Wnt signaling in fgf3-/- by crossing them with hs:dkk1 fish. hs:dkk1 significantly inhibits proliferation in siblings and it also reduces proliferation and hair cell numbers in fgf3-/- (Figure 8A–F). Therefore, in homeostatic neuromasts fgf3 inhibits proliferation by inhibiting Wnt signaling, possibly via sost. As Notch signaling also inhibits proliferation via inhibiting Wnt signaling, we wondered if Notch and Fgf act in the same pathway. Fgf signaling does not act upstream of Notch signaling, as Notch pathway members are not affected by in situ hybridization in fgf3-/- (Figure 8G–K’, Figure 9). fgf3 on the other hand, is slightly downregulated after Notch signaling inhibition with a γ-secretase inhibitor (Figure 8M’) while fgf10a and fgfr1a are not affected (Figure 8O–P’). It is therefore possible that Notch signaling plays some minor role in inducing/maintaining fgf3 expression. However, Notch signaling does not inhibit Wnt signaling via the upregulation of Fgf3 signaling, as shown by experiments in which we constitutively activated Notch signaling in fgf3-/- by expressing the Notch intracellular domain NICD. Activation of NICD in siblings during regeneration strongly inhibits proliferation and hair cell production (Figure 8Q,R,U,V,; ). Activation of NICD in fgf3-/- also strongly inhibits proliferation and hair cell numbers demonstrating that Notch signaling does not require Fgf to inhibit Wnt signaling (Figure 8S–V).

30We 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 (Figure 9). Thus Notch and Fgf need both to be independently and transiently downregulated for efficient hair cell regeneration.

31The scRNA-Seq analyses reported here have characterized in unprecedented detail the different support cell populations present in a homeostatic neuromast. As such, our findings have enabled us to detect subtle gene expression changes in fgf3-/-. Importantly, as dying hair cells are continuously replaced during homeostasis, our analyses also characterized support cells as they differentiated into hair cells. Even though lineage relationships cannot be inferred from scRNA-Seq data alone, the results of our pseudotime and cell classification delineate lineages that have been experimentally confirmed by time lapse and lineage tracing analyses (;;; Figure 2E–F and Figure 3A–C and L–M). The majority of hair cells originate from central support cells without bias to any of the poles (Figure 3M, clusters 7, 9, 8, 14 and 4), whereas amplifying divisions are strongly biased towards the D/V poles and occur in the periphery adjacent to mantle cells (Figure 3M, clusters 10 and 3).

32Interestingly, the PAGA graph (Figure 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. 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.

33Based on cell behavior and gene expression studies we postulate that lateral line neuromasts possess active and quiescent stem cells. Because amplifying support cells in the D/V poles self-renew continuously and only differentiate if displaced into the center, we postulate that amplifying cells represent active stem cells. In addition, neuromasts possibly possess two populations of quiescent stem cells. The first population are cells in the A/P compartment that are relatively quiescent during homeostasis and regeneration (Figure 2H,. These cells start to proliferate after manipulations of the Fgf, Notch or Wnt pathways (Figure 6D, ). Long term lineage tracing experiments have to be performed to assess if they indeed present stem cells. The second quiescent, long term stem cell population are the mantle cells that give rise to hair cells in long term lineage analyses but that do not respond to acute neomycin-induced hair cell death. Support cells also give rise to mantle cells when mantle cells are ablated arguing that active and quiescent stem cells can convert back and forth as shown in other systems. For example, in the intestine slow-cycling, label-retaining stem cells in the +4 position can give rise to stem cells at the bottom of the crypt, which conversely can also give rise to +4 stem cells.

34Quiescent stem cells in hematopoiesis, the brain, skeletal muscle, hair follicle and Drosophila germline stem cells are also distinguished by their low level of protein synthesis and transcription. A characteristic of stem cell activation is consequently the upregulation of protein synthesis and associated factors, such as ribosomal proteins. A heatmap of rpl genes and rps genes shows that in neuromasts these genes are lowly expressed in central support cells (clusters 7, 8, 9, 14), mantle cells (cluster 5, 6) and A/P cells (cluster 13) but are upregulated in dividing cells (clusters 3 and 4) and in D/V cells (clusters 11, 12; Figure 3B’, C’, Figure 3—figure supplement 1). Synthesis of these ribosomal proteins is completely shut down in mature hair cells (cluster 1; Figure 3C’). Based on their low ribosomal protein synthesis central support cells and mantle cells resemble quiescent stem cells, whereas D/V cells and dividing cells resemble activated stem cells or progenitor cells.

35The notion that some support cells constitute stem cells is also supported by the finding that their gene expression profiles share a number of genes with stem cells in the CNS, heart, intestine or hair follicles. For example, neural stem cells in the subventricular zone of the CNS and radial glial cells are characterized by glutamate aspartate transporter (GLAST) and Prominin-1 (also known as CD133) expression. In neuromasts, slc1a3a/glasta and prominin1a are expressed strongly in central and A/P support cells and prominin1b is specifically enriched in A/P cells. Likewise, genes associated with intestinal, hair follicle and neural stem cells (NSCs), such hopx and alpl are strongly expressed in mantle, A/P, and amplifying cells adjacent to mantle cells, but are downregulated in central support cells, forming a ring of expression similarly to fndc7rs4 (Figure 1O and Figure 3G,; ). Genes specifically expressed in central support cells beneath the hair cells likewise express genes that label stem cells in other organs, such as fabp7/Blbp, which labels glioma stem cells, radial glia cells, NSCs, and isl1 that is expressed in quiescent intestinal stem cells and stem cells in the heart. However, in lateral line neuromasts central cells only give rise to hair cells and we therefore do not consider them to be stem cells.

36A 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.

37Wnt signaling induces and is required for support cell proliferation in neuromasts. In fgf3-/- Wnt signaling is upregulated causing cells to proliferate (Figure 8). We found that of the neuromast-expressed secreted Wnt inhibitors sfrp1a, dkk2, sost and sostdc1b, only sost is downregulated in fgf3-/- (Figure 7D’, data not shown,; ). Since sost is also downregulated after hair cell death (Figure 7F’) it may have an important role in inhibiting hair cell regeneration. As Wnt signaling also induces proliferation in mammalian hair cell progenitors a role for Sost paralogs in this system should be tested. According to several data bases, Sost is not, or only lowly expressed in the mammalian ear sensory epithelia, whereas Sostdc1 is robustly expressed.

38Fgf signaling could also inhibit proliferation by regulating other receptor tyrosine kinase activities, for example via inhibition of EGFR by Sprouty2. sprouty2 is also down regulated in fgf3-/- after hair cell death and could perform a similar inhibitory function in homeostatic neuromasts (Supplementary file 14). Fgf signaling also inhibits proliferation in the regenerating utricle and the basilar papilla of chicken. Likewise, addition of Fgf2 or Fgf20 to auditory or vestibular cultures inhibits support cell proliferation, suggesting that the inhibitory effect of Fgf signaling on progenitor proliferation is evolutionary conserved in species that can regenerate their sensory hair cells.

39Experiments 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. 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.

40While 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. Likewise, in the mouse cochlea loss of Fgfr3 results in increased hair cells at the expense of support cells, 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. 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.

41Just 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. 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-/- (Figure 8G–K’). Also, while pharmacological inhibition of Notch activity during homeostasis decreases fgf3 expression it has no effect on Fgf target genes (Figure 8L–P’). Additionally, activating Notch signaling by over-expression of a Notch intracellular domain inhibits proliferation and differentiation during regeneration in both wild type and fgf3-/- (Figure 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. 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 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.

42The scRNA-Seq analysis revealed previously unidentified support cell populations and combined with in situ validation of these cell clusters 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. Such a comparison will become even more powerful once adult mouse single cell transcriptomes of support cells are available.

43ET(krt4:EGFP)SqGw57a, Tg(pou4f3:GAP-GFP)s356t, fgf3t26212, Tg(hsp70l:dkk11b-GFP)w32tg, Tg(hsp70l:MYC-notch1a,cryaa:Cerulean)fb12Tg, Tg(hsp70l:dnfgfr1a-EGFP)pd1tg, TgBAC(prox1a:KALTA4,4xUAS-ADV.E1b:TagRFP)nim5Tg, Tg(EPV.Tp1-Mmu.Hbb:EGFP)um14, Et(krt4:EGFP)sqet4ET and Et(krt4:EGFP)sqet20ET. Tg(6xNFkB:EGFP)nc1. fgfr1asa38715 and fgfr2sa30975 are from the Sanger Institute Zebrafish Mutation project. All experiments were performed according to guidelines established by the Stowers Institute IACUC review board.

44H2B-mturquoise2 was placed near the 5’ region of fgf3 using non-homologous repair with CRISPR/Cas9. A CRISPR recognition site (GGCCATGGAAACTAAATCTGCGG) was chosen 584 bp in front of the fgf3 transcription start site using CRISPRscan. The same recognition site was cloned by PCR onto both ends of a construct containing a 56 bp β-actin minimal promoter driving human-H2B fused at the c-terminus with mturquoise2 followed by the SV40 polyA from the Tol2 kit. This plasmid was mixed with the crRNA (GGCCATGGAAACTAAATCTG, IDT), tracrRNA (IDT) and Cas9 protein (PNA Bio) and complexed for 10 min at room temperature then placed on ice. The complex was then injected into the cell of a one cell stage zebrafish embryos from a cross of fgf3t26212 to wild type. 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. The founder generated also has the fgf3t26212 allele and is therefore heterozygous for the mutation. These fish are viable and fertile and show no obivious phenotypes.

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

46BrdU (Sigma-Aldrich) was added at 10 mM with 1% DMSO in E2. Larvae were treated for 24 hrs then fixed in 4% paraformaldehyde overnight at 4°C. Mouse anti-BrdU (Sigma-Aldrich) and rabbit anti-GFP (Invitrogen) immunohistochemistry with DAPI (Invitrogen) counterstain was performed as described. EdU (Carbosynth) was added at 3.3 mM with 1% DMSO in E2. Larvae were treated for 24 hrs then fixed in 4% paraformaldehyde overnight at 4°C. For staining, larvae were washed 3 times 10 min each in PBS/0.8% Trition-X (PBSTX), blocked for 1 hr in 3%BSA/PBSTX, washed 3 times 5 min each in PBS then put in fresh staining solution for 30 min. Staining solution contains 1xTris buffered saline, 0.1% Tween-20, 2 mM CuSO4, 2.5 mM Alexa-594-Azide (Invitrogen) and 50 mM ascorbic acid. After staining, larvae were washed extensively in PBSTX. Larvae were then processed for anti-GFP immunohistochemistry and DAPI staining as described above. Stained larvae were imaged on a Zeiss LSM780 confocal microscope at 40X. Three posterior lateral line neuromasts (L1, L2 and L3) were imaged per fish. Cell numbers were manually counted in Imaris software. Spatial positioning was performed as described.

47In situ hybridization was performed as described with modifications and http://research.stowers.org/piotrowskilab/in-situ-protocol.html. Incubation time in proteinase K (Roche) depended on the batch, and varied from 2 to 5 min at room temperature. Pre-hybridization was performed for at least 2 hrs at 65°C. Probes used were fgfr1a, fgf3, fgf10a, and etv4, fgfr2, dld, wnt2, wnt10a, wnt11r, gfp, atoh1a and notch3 and sfrp1a Additional probes were generated by PCR from zebrafish cDNA and cloned into the topo-pCRII vector (Invitrogen). See Data file 16 for list or primers used. Alkaline phosphatase staining was performed as described.

48Heat shock induction of transgene expression varied depending on the transgenic line used. Initially 5dpf larvae were heat shocked for 1 hr at 37°C then put back at 28°C for 1 hr. Larvae were then heat shocked for 1 hr at a higher temperature (39°C for notch1a-intraceullar 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 1 hr 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 min 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 2 hrs recovery at 28°C. The 1 hr heat shock followed by 2 hrs recovery was repeated for a total of 24 hrs post neomycin treatment. Larvae were then fixed in 4% paraformaldehyde overnight at 4°C.

49TgBAC(prox1a:KALTA4,4xUAS-ADV.E1b:TagRFP) were crossed to Tg(pou4f3:GAP-GFP)s356t and larvae were raised to 5dpf. Time-lapse imaging after neomycin treatment was carried out on a Zeiss LSM780 confocal microscope as described. Tg(fgf3:H2B-mturquoise2)psi60Tg/Tg(pou4f3:GAP-GFP)s356t or Tg(fgf3:H2B-mturquoise2)psi60Tg/TgBAC(prox1a:KALTA4,4xUAS-ADV.E1b:TagRFP) double transgenic larvae were imaged live as above. Three-dimensional rendering and image analyses were done using Imaris (Bitplane).

505 dpf Tg(pou4f3:GAP-GF)/Tg(GW57a:GFP) zebrafish embryos were dissociated by adding 1.5 ml (per 100 embryos) of 0.25% Trypsin-EDTA (Thermo Fisher Scientific, Waltham, MA. USA) and triturated with 1 ml pipette tip for 3 min on ice. To collect dissociated cells, cells were filtered with 70 μm Cup Filcons (BD Biosciences, San Jose, CA. USA) and washed with ice-cold DPBS (centrifugation at 2000 rpm for 6 mins at 4°C). Cells were stained with Hoechst 33342 (final concentration: 0.005 μg/μl) or Draq5 (1:2000) (biostatus, UK) and 7-AAD (final concentration: 0.5 μg/ml) or DAPI (5 μg per ml) to gate out dead cell populations. FACS was performed at the Cytometry Core facility (Stowers Institute for Medical Research) using BD Influx Cell Sorter (BD Biosciences, San Jose, CA. USA).

51scRNA-seq was carried out with 10X Chromium single cell platform (10X Genomics, Pleasanton, CA. USA). FAC-sorted live or MeOH-fixed cells were used as the input source for the scRNA-seq. MeOH-fixed cells were rehydrated with rehydration buffer (0.5% BSA and 0.5 U/μl RNase-inhibitor in ice-cold DPBS) following manufacturer’s instructions (10X Genomics). The maximum recommended volume of single cell suspension (34 μl) was loaded on a Chromium Single Cell Controller (10x Genomics) targeting ~1500–2000 cells per sample. Chromium Single Cell 3' Library and Gel Bead Kit v2 (10X Genomics) was used for libraries preparation according to manufacturer’s instructions. Resulting short fragment libraries were checked for quality and quantity using an Agilent 2100 Bioanalyzer and Invitrogen Qubit Fluorometer. Libraries were sequenced to a depth of ~160–330M reads each on an Illumina HiSeq 2500 instrument using Rapid SBS v2 chemistry with the following paired read lengths: 26 bp Read1, 8 bp I7 Index and 98 bp Read2.

52Raw 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 wild type and version 2.1.1 for fgf3-/- data sets). 1,666 cell barcodes were obtained for wild type embryos, 1932 for fgf3 siblings and 1459 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).

53Subsequent analyses on UMI count matrix for all three data sets were performed using the R package Seurat (version 2.3.4, ) 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 five 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 six principal components (PCs) for wild type 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 two-dimensional visualization, we performed a second round of dimensional reduction using t-SNE. The dendrogram in Figure 3A was generated using the Seurat function BuildClusterTree. This function performs hierarchical clustering based on a Euclidian distance matrix, which is created from the average expression of the variable genes (the same set used for dimensional reduction) in each of the 14 clusters. 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 non-specific cell types that contained markers for skin and blood cells which were removed from the final analysis.

54All data contained within our processed Seurat object for the wild type data set was converted to the AnnaData format for pseudotime analysis in Scanpy (version 1.2.2, ), 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 one 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).

55The GO term analysis was performed in DAVID (Database for Annotation, Visualization and Integrated Discovery,.

56The web app (https://piotrowskilab.shinyapps.io/neuromast_homeostasis_scrnaseq_2018/) for our homeostatic data set was generated using the Shiny framework in R. Shiny allows one to create a web-based GUI (graphical user interface) for data exploration using R code (the source code for our app can be found at https://github.com/Piotrowski-Lab/Shiny-Apps;; copy archived at https://github.com/elifesciences-publications/Shiny-Apps). The user can choose a zebrafish gene(s) of interest, provided it is expressed in our dataset, and view expression values on a heatmap, violin, or t-SNE plot. Instructions are provided in the welcome page of the web 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).

57Total 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.

58Q-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 min), quantitation (40 cycles of 95°C for 15 s, 60°C for 60 s) and melting curve analysis (95°C for 15 s, 60°C for 60 s, 95°C for 15 s). The experiment was performed in triplicate. All the signals were normalised to the ef1α expression level. All primer sequences are provided in Data file 15.