All zebrafish strains were bred and housed in the Danio rerio research facilities within the University of Melbourne and the Walter and Eliza Hall Institute of Medical (ethics approval ID 22235 and 10400), in accordance with local guidelines. All procedures were approved by the Faculty of Science Animal Ethics Committee at the University of Melbourne (approval IDs 1,814,542 and 10,232). Zebrafish embryos were raised at 28.5°C in E3 medium (5 mM NaCl, 0.17 mM KCl, 0.33 mM CaCl₂, 0.33 mM MgSO₄) at a maximum density of n = 50 per 50 mL petri dish. E3 medium was supplemented with 0.003% 1-phenyl-2-thiourea (PTU) from 24 h post-fertilization (hpf) to prevent pigmentation and to maintain transparency of zebrafish embryos, prior to sorting for fluorescence signal. Prior to sorting for experiments, transgenic Tg(gfap:GFP) larvae were temporarily anesthetized in 0.04% tricaine methanesulfonate (MS-222; Sigma-Aldrich, cat. Number E10521-50G) and could be sorted for GFP signal at any age over 24 hpf. Larvae for Tg(opn1lws2:nfsb-mCherry)^uom3, Tg(opn1sws2:nfsb-mCherry)^uom4, and Tg(xops:nfsb-mCherry)^uom5 strains were treated with the same anesthetic bath to allow for sorting of the nfsb-mCherry signal in larvae older than 72 hpf.

DNA constructs were injected into one-cell stage wildtype AB zebrafish embryos, using a FemtoJet microinjector (Eppendorf) and borosilicate glass capillary needle (1.0 mm O.D./0.78 mm I.D./100 mm long capillary). Tg(gfap:GFP)^m2001 zebrafish were generated by Bernardos and Raymond (2006). The photoreceptor promoters were chosen to specifically ablate the subtype expressing the distinct opsin protein, which is a hallmark of the different photoreceptor types. The lws2 promoter (1.77 kbp) and xops promoter (1.38 kbp, Xenopus rhodopsin) were amplified by PCR using specific primers. The primers for the lws2 promoter were 5′-GGCCAGATGGGCCCTGTTGTGCACCAGATCTGAGT-3′ and 5′-TGGTCCAGCCTGCTTTTTGGAAACCCTGAAGATCA-3′, while the xops promoter was amplified from pFIN-XOPS-tdTOMP using forward (5′-TATAGGGCGAATTGGGGCCGCAGATCTTTATACATTGC-3′) and reverse (5′-CCGGTGGATCCCAAACCCTCGAGATCCCTAGAAGCCTGTGAT-3′) primers. The products were subcloned into the pTol-uas:nfsB-mCherry plasmid using the ClonExpress MultiS One Step Cloning Kit to replace the UAS promoter by homologous recombination (Wang et al., 2020). The pTol-uas:nfsB-mCherry plasmid was a gift from Prof. Toshio Ohshima’s lab at Waseda University, Tokyo, Japan. The pTol-sws2:nfsb-mcherry plasmid was a gift from Prof. Rachel Wong at University of Washington, Seattle, Washington, United States (Yoshimatsu et al., 2016). The resulting plasmids, pTol-lws2:nfsB-mCherry, pTol-xops:nfsB-mCherry, and pTol-sws2:nfsb-mcherry, were co-injected with Tol2 transposase mRNA into AB embryos at the one-cell stage. F1 was identified by screening for the mCherry signal. The pFIN-XOPS-tdTOM was a gift from Susan Semple-Rowland (Addgene plasmid # 44359¹; RRID:Addgene_44359).

Plasmid linearization of Tol2 DNA was achieved using a Quiaquick Kit, while RNA synthesis was carried out using the Ambion message machine kit Sp6.

Zebrafish lines Tg(her4.1:dRFP/gfap:GFP) was used for control (non-ablated) tissue, while Tg(her4.1:dRFP/gfap:GFP/lws2:nfsb-mCherry) and Tg(her4.1:dRFP/gfap:GFP/sws2:nfsb-mCherry) were used to obtain long wavelength sensitive (Lws2) cone-ablated and short wavelength sensitive (Sws2) cone-ablated retinal tissue, respectively. Zebrafish larvae from all lines at 6 days post-fertilization (dpf) were exposed to a 10 mM solution of metronidazole for 48 h, either as non-injury control (no nfsb-mCherry) or to specifically eliminate the red or blue cone photoreceptor. At 8 dpf, fish were rinsed in fresh system water and housed in standard conditions. At 9 dpf (3 days post-injury; dpi), fish were humanely killed. Single-cell suspensions of 9 dpf zebrafish (age-matched control and post-injury larvae) were prepared using a specific protocol (Lopez-Ramirez et al., 2016). Retinae were dissected and digested in 350 μL papain solution at 37°C for 15 min. The papain solution was prepared as follows: 100 μL papain (Worthington, LS003126), 100 μL of 1% DNAse (Sigma, DN25), and 200 μL of 12 mg/mL L-cysteine (Sigma, C6852) were added to a 5 mL DMEM/F12 (Invitrogen, 11,330,032). During digestion, retinal tissue was mixed by repeated pipetting 4 to 10 times. Following digestion, 1,400 μL of washing buffer was added, containing 65 μL of 45% glucose (Invitrogen, 04196545 SB), 50 μL of 1 M HEPES (Sigma, H4034), and 500 μL FBS (Gibco, 10,270,106) in 9.385 mL of 1x DPBS (Invitrogen, 14,190–144). All solutions were filtered through a 0.22 μm filter (Millipore) and stored at 4°C prior to use.

To perform single-cell RNA sequencing (scRNA-seq), cells isolated through fluorescent activated cell sorting (FACS) were loaded onto a Chromium Single Cell Chip (10x Genomics, United States) according to the manufacturer’s protocol. The perdurance of the GFP in the Tg(gfap:GFP)^m2001 enables us to track all glia derived cells including those that have de-differentiated into progenitors as well as those that start differentiating into mature cells (glia and neurons). The scRNA-seq libraries were generated using the GemCode Single-Cell Instrument and Single Cell 3’ Library and Gel Bead kit v2 and v3 Chip kit (10x Genomics, 120,237). Library quantification and quality assessments were performed by Qubit fluorometric assay (Invitrogen) and dsDNA High Sensitivity Assay Kit (AATI, DNF-474-0500). Analysis of DNA fragments was performed using the High Sensitivity Large Fragment -50 kb Analysis Kit (AATI, DNF-464). The indexed library was tested for quality and sequenced using an Illumina NovaSeq 6,000 sequencer with the S2 flow cell using paired-end 150 base pair reads. Sequencing depth was 60 K reads per cell.

Filtered matrix raw data files were analyzed in R using Seurat. Prior to downstream filtering, the following describes cell numbers for datasets: 8325 (no ablation control), 10,520 (Lws2-ablation) and 7,851 (Sws2-ablation). Post-filtering and exclusion of non-Müller glia-derived cell types achieved the following cell numbers: 7884 (no ablation control), 9,789 (Lws2-ablation), and 7,542 (Sws2-ablation). Low quality cells or cells containing doublets were excluded from all datasets; reads between 200 and 3,500 genes per cell (no ablation control), reads between 200 and 4,500 genes per cell (Lws2-ablation), and reads between 200 and 4,000 genes per cell (Sws2-ablation) were obtained. Cells with a percentage of mitochondrial gene expression of greater than 20% (no ablation control), 35% (Lws2-ablation), and 30% (Sws2-ablation) were excluded from further analysis.

Samples were analyzed using Seurat::NormalizeData, while variable features for downstream analysis were identified using Seurat::FindVariableFeatures, and scaled using Seurat::ScaleData. An optimal number of principal components (PCs), generated through Seurat::RunPCA, for dimensional reduction were selected using the function Seurat::ElbowPlot. PCs containing the greatest variance were selected. For each sample, 20 PCs were specified, and we selected 13 (no ablation control), 15 (Lws2-ablation) and 15 (Sws2-ablation) PCs. Clustering was performed using the shared nearest neighbor (SNN), graph-based approach through Seurat::FindNeighbours. Uniform manifold approximation and projection (UMAP) plots containing 11 (no ablation control), 12 (Lws2-ablation) and 12 (Sws2-ablation) clusters were prepared using Seurat::FindClusters.

Samples were integrated by a standard integration protocol for R package Seurat using Seurat::FindIntegrationAnchors and Seurat::IntegrateData. Downstream analysis including data scaling, dimensional reduction and clustering was conducted as per above. A total of 13 PCs (no ablation control vs. Lws2-ablation) and 15 PCs (Lws2-ablation vs. Sws2-ablation) were selected for dimensional reduction.

Markers driving the characterization of each cluster were resolved using Seurat::FindMarkers, with only positive features expressed in at least 25% of cells in each group of cells. Lists of differentially expressed genes generated through this approach were used for gene set enrichment analysis (GSEA), conducted using gprofiler::gost, with statistical significance evaluated using a Benjamini–Hochberg FDR test set at a threshold of 0.05. Relevant terms were selected and plotted in accordance to their negative adjusted p value, along with gene ratio (number of genes attributed to the relevant GO term/number of genes in the cluster gene list). Genes mentioned from this cluster characterization analysis were within the top 50 markers.

Pseudotime analysis was conducted using the R package monocle3. Prior to analysis of trajectory and searching for differentially expressed genes, ribosomal protein genes or mitochondrial genes were identified from lists and excluded. The Seurat object was converted to a cell data set (cds), and a subset of cells from quiescent to activated Müller glia cell states was selected using Monocle3::choose_graph_segments and Monocle3::choose_cells. Next, Monocle3::order_cells was used to select the root node, which was in the quiescent cell cluster furthest from the activated cell population in the UMAP plot. In Monocle3::graph_test, we used the principal graph option with Moran’s test statistic to identify significantly differentially expressed genes along the specified trajectory. The q value was set to <0.05 and Moran’s I value set to >0.25. The top 40 differentially expressed genes are presented in a heatmap, hierarchically clustered based on gene modules generated from Monocle::graph_test.

RNAscope probes for efnb2a, fgf24, and rdh10a were generated by ACD company (Shanghai). The experiments were carried out as follows. Larval zebrafish were fixed in 4% PFA at 4°C overnight followed by cryoprotection in 30% sucrose and then were cryosectioned at a thickness of 12 μm. The slices were post-fixed in 4% PFA at room temperature for 15 min and washed with 1 × PBS at room temperature for 3 min. To block the activity of endogenous peroxidase, all slides were treated with 0.1% H₂O₂ at room temperature for 30 min. After being washed twice with 1 × PBS at room temperature for 3 min, slides were treated with 10 μg/mL proteinase K (Sigma) diluted in TE (10 mM Tris–HCl, pH 8.0, and 1 mM EDTA, pH 8.0) at 37°C for 8 min, then treated with 4% PFA at room temperature for 10 min. Subsequently, all slides were washed with 1 × PBS at RT for 3 min, followed by the incubation in 0.2 M HCl at RT for 10 min. After washing with 1 × PBS for 5 min, all slices were then incubated with 0.1 M triethanol amine-HCl (662.5 μL triethanolamine and 1.35 mL 1 M HCl; adding water to the final volume of 50 mL, pH 8.0) at room temperature for 1 min and in 0.1 M triethanol amine-HCl containing 0.25% acetic anhydrate at room temperature for 10 min with gentle shaking. Slides were then washed by 1 × PBS at room temperature for 5 min then dehydrated in a series of 60, 80, 95% ethanol baths, and finally twice in 100% ethanol at room temperature for 90 s, respectively. Slides were incubated in the hybridization buffer (50% formamide (Sigma), 10 mM Tris–HCl, pH 8.0, 200 μg/mL yeast tRNA (Invitrogen), 1 × Denhart buffer, SDS, EDTA and 10% dextran sulfate (Ambion) containing 1 μg/mL probes at 60°C overnight). On the second day, slides were washed sequentially in 5 × SSC at 65°C for 30 min, 2 × SSC with 50% formamide at 65°C for 30 min, TNE buffer (100 mL TNE consisting of 1 mL 1 M Tris–HCl, pH 7.5, 10 mL 5 M NaCl, and 0.2 mL 0.5 M EDTA) at 37°C for 10 min and then in TNE buffer with 20 μg/mL RNaseA at 37°C for 30 min. Slides were then incubated with 2 × SSC at 60°C for 20 min, 0.2 × SSC at 60°C for 20 min, and 0.1 × SSC at RT for 20 min. Next, slides were blocked by TN buffer at room temperature for 5 min (200 mL TN buffer consisting of 20 mL 1 M Tris–HCl, pH 7.5, 6 mL 5 M NaCl, and 174 mL water) followed by TNB buffer (TN buffer and 0.5% blocking reagent; Roche) at room temperature for 5 min. Finally, slides were incubated in TNB buffer with anti–DIG-POD (1:500; Roche) at 4°C overnight. On the third day, the signal was detected by the TSATM Plus Cyanine 3/Fluorescein System (PerkinElmer, NEL753001KT).

For Tg(lws2:nfsb-mCherry), Tg(sws2:nfsb-mCherry), or Tg(xops:nfsb-mCherry) zebrafish larvae processed for immunohistochemical analysis, metronidazole treatment was conducted at 4 dpf for 48 h, leading to ablation of either Lws2 cones, Sws2 cones, or rod photoreceptors, respectively. Zebrafish larvae were exposed to a solution of 10 mM metronidazole (Sigma-Aldrich, cat. Number M3761-100G) in standard fish water. To induce a lesser degree of injury in Tg(lws2:nfsb-mCherry) zebrafish, larvae were swum in a solution of 5 mM metronidazole in standard fish water for 2 h. Zebrafish larvae were placed into solution at maximum densities of n = 50 larvae per 50 mL petri dish and kept at 28.5°C for the duration of treatment. Control larvae were placed in standard fish water. Zebrafish were rinsed in standard fish water after treatment. For adult zebrafish, these were immersed in a solution of 10 mM metronidazole in standard fish water for 3 consecutive days, while control adult zebrafish were placed in standard fish water, both with these respective baths changed twice a day.

Zebrafish larvae were humanely killed in 4 g/L tricaine, fixed in 4% paraformaldehyde (PFA; Sigma-Aldrich, cat. Number 158127-500G) in PBS overnight at 4°C and subsequently rinsed in PBS. Larvae were then cryoprotected in a 30% sucrose solution in 1x PBS overnight at 4°C before being embedded in OCT (Tissue-Tek) in 10 mm/10 mm/5 mm Tissue-Tek Cryomolds (ProSciTech). Molds containing these tissues were frozen and stored at −20°C. Samples were sectioned at 12 μm thickness using a Microtome Blade (Arthur Bailey Surgico) on a Leica (CM1860) cryostat, with sections transferred onto room-temperature 25 mm/75 mm Menzel Superfrost Plus glass slides (Grale Scientific), allowed to dry for 1 h, and used immediately for immunohistochemistry or stored at −20°C.

Antibody staining was carried out at room temperature using standard protocols. Antigen retrieval to detect epitopes for antibodies was performed by incubating slides in 150 mM Tris–HCl (pH 9) at 70°C for 20 min, then rinsed for 30 min before overnight incubation in a primary antibody prepared with 5% fetal bovine serum (FBS) blocking solution. Primary antibodies used for this study were as follows: mouse anti-proliferating cell nuclear antigen (PCNA; Santa Cruz Biotechnology, sc-25,290), 1:500; rabbit anti-mCherry (Invitrogen, PA5-34974), 1:500, rabbit anti-Bmpr1b (GeneTex, GTX128200). The next day, slides were rinsed thrice with PBS and then incubated for 2 h in secondary antibody, diluted in the same blocking solution. Secondary antibodies used were as follows: Alexa Fluor 647 nm donkey anti-mouse, 1:500; Alexa Fluor 546 nm goat anti-rabbit, 1:500. Following antibody staining, slides were rinsed thrice with PBS and then stained with 4′,6-diamidino-2-phenylindole (DAPI; 1:10000, Sigma-Aldrich, D9542-10MG) in PBS. Finally, sections were mounted in Mowiol (Sigma-Aldrich, cat. Number 81381-250G) using 24 mm/60 mm glass coverslips (ProSciTech). Slides were stored at 4°C.

Terminal deoxynucleotidyl transferase dUTP nick end labeling was carried out on sections (Roche, 11,684,795,910). All steps were performed at room temperature unless otherwise indicated. Tissue sections were post-fixed with 4% PFA in PBS for 20 min. Slides were rinsed twice in PBS and incubated in 0.1% triton X-100/0.1% sodium citrate buffer for 20 min and rinsed twice in PBS. Following 5 min of air drying, slides were incubated in TUNEL reaction mixture in a dark humid chamber at 37°C for 1 h. Slides were rinsed twice with PBS and then stained with DAPI (1,10,000) in PBS and coverslipped in Mowiol.

Images of stained immobilized sections were captured using a Nikon A1 confocal microscope with a 40 x pan-fluor oil-immersion objective lens, with 1.3 numerical aperture. Z-stacks were acquired with a step size of 1 μm. For RNAscope samples, images were taken using an inverted confocal microscope system (FV1200, Olympus) using a 30 x (silicon oil, 1.05 NA) or 60 x (silicon oil, 1.3 NA) objective lens. For each fish, an image of a single retina was obtained. Quantification was conducted using FIJI/ImageJ. To define dorsal, central and ventral regions of the retina, three segments of 60° were partitioned for quantification studies therein. For TUNEL and mCherry comparison between uninjured and each of the injury paradigms, we used FIJI to draw ROIs in dorsal, central or ventral regions of the retinal sections. The loss of mCherry after injury was quantified as pixels/area.

The number of TUNEL labeled cells was manually quantified within each ROI.

Data are expressed as mean ± standard error. The number of fish in each condition are specified in the respective graphs. Statistical analyses were conducted using Prism 8 (GraphPad) using a two-way ANOVA followed by a Bonferroni post hoc test for multiple comparisons. To determine significance between mCherry-positive glia in the activated population compared to the total glia analyzed, a chi-square test of independence was conducted. To determine significance in TUNEL+ and mCherry comparison before and after injury, the unpaired student’s t-test was used (with Welch’s correction applied in case of unequal variances). Statistical significance are as follows: * = p ≤ 0.05, ** = p ≤ 0.01, *** = p ≤ 0.001, or **** = p ≤ 0.0001.

Figures were compiled using Adobe Illustrator. Microscope images were arranged in Adobe Photoshop, where brightness was adjusted evenly across channels and images, prior to placement into Illustrator.

Following neural ablation or cell death in the zebrafish retina, Müller glia that proliferate must first undergo reprogramming from a quiescent to an activated state. First, we characterized the gene expression networks which underlie the switch from quiescence to activation following two injury paradigms; specific ablation of long wavelength sensitive cone photoreceptors (Lws2) in Tg(lws2:nfsb-mCherry), or short wavelength sensitive cone photoreceptors (Sws2) in Tg(sws2:nfsb-mCherry) zebrafish. Expression of nitroreductase (nfsb) in these cell populations specifically leads to cytotoxicity upon metronidazole exposure (Curado et al., 2007, 2008). We combined our injury lines with Tg(gfap:GFP) zebrafish and performed single cell RNA sequencing (scRNA-seq) of fluorescent activated cell sorting (FACS) Gfap-positive cells, 3 days post injury (dpi) in larval zebrafish. By integrating scRNA-seq datasets of Müller glia from Lws2-ablated and Sws2-ablated retinas, we sought to determine commonalities and differences in the Müller glia response (Figure 1). Regardless of injury, quiescent and proliferative clusters were conserved across both datasets. However, the main differences were observed in clusters containing differentiating Müller glia-derived progenitor cells, expressing gria2b, pou3f3b, elav3 and elav4, all markers of neuronal differentiation (Figures 1A,B). Understanding the molecular pathways that guide the activation of stem cells is of great interest for studying regeneration. To investigate this, we isolated quiescent and activated Müller glia clusters and constructed a pseudotime trajectory (Figures 1C,D) to determine the activation of gene modules which accompany the changes between these two cellular states (Figures 1C–G).

Eight temporally expressed gene modules along pseudotime were observed (Figure 1E). The first and second modules showed expression of genes related to mature Müller glia (glula, glulb, and rlbp1a) and retinal homeostasis (slc1a2b, apoeb, and atp1a1b). Fos and Jun-family members, transcriptional regulators that dimerize to form AP-1 transcriptional regulators involved in regulation of cell proliferation and differentiation, were also highly expressed in these Müller glia. Additionally, cdo1 and cebpd associated with cell proliferation and differentiation were upregulated, as well as mlc1a, associated with Müller glia resting state across vertebrate species (Hoang et al., 2020). However, as Müller glia transitioned to an activated state, they repressed expression of these genes, and temporally activated a module including transient expression of genes including notch3 and her12, id1 (Inhibitor of DNA binding), the glial differentiation gene meteorin (metrn) and the tight junction protein-encoding gene tjp2b (Figure 1G, Supplementary Figure 1).

The third module of gene expression involved ribosomal genes faua and rack1, si:dkey-151 g10.6, and Eef-family genes (eef1g, eef1a1l1, eef1b2 and eef2b). These genes, which are induced as Müller glia dedifferentiate (lose glia marker expression) and enter the cell cycle, are likely important in driving this regenerative cascade (Figures 1E,G, Supplementary Figure 2). Upon entering the cell cycle, Müller glia upregulate ribonucleotide reductase encoding genes rrm1 and rrm2.1 (module four), as well as the S-phase markers pcna (proliferating cell nuclear antigen) and dut (deoxyuridine triphosphatase), in module five. Next, the following two modules six were enriched for genes relating to chromatin remodeling, expressing histone mobility group (HMG) genes hmgb2e, hmgn2, hmgb2b, hmga1a as well as other histone-associated genes (h2afvb and h3f3a), as well as the cytoskeleton, expressing alpha tubulin (tuba1a) and beta tubulin (tubb2b and tubb4b) genes. Finally, module eight reflects Müller glia progressing into G2/M phase, with expression of cdk1, cdc20, nusap1 and the kinesin family member genes kifc1 and kif11 (Figure 1E, Supplementary Figure 2). Thus, these data suggest functional competencies adopted by Müller glia that sequentially underpin gene expression as they exit their quiescent state to move toward a proliferative state.

Previous studies indicate that the recruitment of Müller glia to enter this neurogenic pathway is dependent on the proximity to, and size of neural injury (Wan et al., 2014; Powell et al., 2016; Ng Chi Kei et al., 2017). However, whether all Müller glia are capable of proliferating and subsequently contributing to regeneration has not been systematically tested. We generated a variety ablation models that target specific photoreceptor subtypes (Figures 2A–H), which differ in their relative abundances and distribution across the larval zebrafish retina (Noel et al., 2021). These were the previously mentioned cone photoreceptor ablation models Tg(lws2:nfsb-mCherry) and Tg(sws2:nfsb-mCherry), as well as the rod (Xops) photoreceptor ablation model Tg(xops-nfsb-mCherry). Additionally, we optimized a lesser injury model in Tg(lws2:nfsb-mCherry) fish through weaker chemical induction of photoreceptor damage that was similar in extent to the ablation observed in Sws2 and Xops ablation paradigms, allowing us to minimize differences driven by injury extent. The injury extent across different spatial regions of the retina (correlated to the relative abundance of these different photoreceptors in each of the spatial regions) was quantified by measuring loss of mCherry photoreceptors and increase in TUNEL labeled cells in the outer nuclear layer at 2 dpi (Supplementary Figure 3). As can be seen, most mCherry labeled photoreceptors have started to be cleared from the three larger injury by this time. In the small Lws2 injury, we see individual photoreceptors with abnormal morphology and increased mCherry staining, particularly in the central retina (Figure 2F, Supplementary Figure 3). This might be because in the smaller injury, the reduction in cytoxin generated (less metronidazole prodrug per cell) might delay cell death and cell clearance compared to the other injury. For the cell death quantification using TUNEL, we specifically quantified the number of TUNEL+ in the ONL including colabelled DAPI nuclear staining (Supplementary Figure 3, arrowheads). Of note, some of the TUNEL+ may also be indicative of non-nuclear labeling associated with necrotic cell death or cell lysis (Moore et al., 2001). Signal in the green channel is also observed in the outer segments in the photoreceptors due to background signal in this layer often observed and seen both in uninjured and injured retinal sections.

With this approach, we quantified the percentage of Müller glia positive for proliferating cell nuclear antigen (PCNA) in our Tg(gfap:GFP) zebrafish crossed to either Tg(lws2:nfsb-mCherry), Tg(sws2:nfsb-mCherry), or Tg(xops:nfsb-mCherry) at 48 and 72 h post-injury (hpi), as a readout of the quiescent Müller glia that were activated to drive the regenerative neurogenic response (Figures 3A–H). Following widespread Lws2-ablation (Figures 3A,B), the majority of PCNA-positive, proliferating glia were found in the central sector (24.4 ± 2.8% at 48 hpi, 26.0 ± 3.5% at 72 hpi) and dorsal sector (19.8 ± 2.9% at 48 hpi, 22.1 ± 3.10% at 72 hpi), with no significant differences between these sectors at either time point recorded. In contrast, the proportion of proliferating Müller glia in the ventral sector was significantly reduced when compared to dorsal (p < 0.001) and central (p < 0.0001) Müller glia at both 48 (2.62 ± 1.70%) and 72 hpi (4.88 ± 2.94%), despite Lws2 cones being abundant in this region as well. This same pattern was observed following our smaller Lws2-ablation (Figures 3C,D); dorsal (8.4 ± 2.1%; p = 0.03) and central (12.7 ± 3.1%, p < 0.001) Müller glia showed significantly higher proliferation than ventral Müller glia (1.2 ± 0.9%) at 48 hpi. No significant difference in proliferation was observed at 72 hpi between Müller glia in dorsal (5.2 ± 2.2%), central (7.5 ± 1.7%), and ventral (1.3 ± 0.8%) domains.

Minimal Müller glia proliferation was observed following Sws2-ablation (Figures 3E,F), which has been previously reported in zebrafish (D'Orazi et al., 2020). PCNA-positive Müller glia were only observed at 48 hpi (2.3 ± 0.9% and 2.2 ± 1.1%) and 72 hpi (0.8 ± 0.5% and 1.0 ± 0.7%) in dorsal and central sectors of the retina, respectively. Ventrally located PCNA-positive Müller glia were not observed following Sws2-ablation. After Xops ablation (Figures 3G,H), which are more abundant in the ventral and dorsal regions of the retina at these larval stages, Müller glia activation occurred predominantly at the site of neural cell death consistent with previously published results (Lenkowski and Raymond, 2014). We observed only very minimal proliferation of central Müller glia (1.9 ± 1.0% and 2.4 ± 1.0%) at 48 and 72 hpi, respectively. Indeed, most of the proliferating Müller glia were detected in the dorsal domain (9.4 ± 1.8% at 48 hpi and 11.6 ± 2.1% at 72 hpi). However, ventrally located Müller glia proliferated, at a reduced capacity compared to dorsally located Müller glia (4.8 ± 2.1% at 48 hpi; p = 0.1 and 2.7 ± 1.4% at 72 hpi; p < 0.001).

Following rod ablation, proliferative Müller glia were predominantly located in the peripheral (dorsal and ventral) retina, which matched the overall distribution of rod photoreceptors (mCherry-positive) at this age (Figures 2D,H). Thus, we assessed if the differences in capacity for activation of Müller glia are influenced by signals local to the site of injury. Specifically, we investigated phagocytosis of dying cells and their debris, which can influence behavior Müller glia following injury (Sakami et al., 2019; Lew et al., 2022) including Müller glia proliferation (Bailey et al., 2010; Nomura-Komoike et al., 2020). For this, we quantified the number of mCherry containing presumed phagocytic Müller glia and compared this to the Müller glia re-entering the cell cycle (PCNA-positive) after widespread Lws2, smaller Lws2 and rod ablation, as these resulted in a substantial Müller glia proliferative response (Figure 4). As shown, we found that more Müller glia contained mCherry-positive debris than expressed PCNA in the small Lws2 cone and rod ablation studies, yet more Müller glia re-entered the cell cycle in the widespread Lws2 cone ablation (Figures 4A–L). In all cases, we found a significantly higher proportion of phagocytic (mCherry-positive) Müller glia within the proliferative population compared to the total Müller glia cohort at 48 hpi (22% vs. 13% Lws2 big injury; 69% vs. 48% Lws2 small injury; 87% vs. 46% Xops injury; p < 0.0001) and 72hpi (25% vs. 14% Lws2 big injury; 62% vs. 33% Lws2 small injury; 76% vs. 42% Xops injury; p < 0.0001; Figures 4H,J,L). Thus, Müller glia that phagocytose substantial amounts of dying photoreceptors are significantly more likely to respond to injury and mount a regenerative response, consistent with phagocytosis being one of the cellular processes that may convey injury signals to recruit Müller glia for regeneration.

Additionally, we wanted to detect even minute levels of mCherry debris, which may not have been visible after standard immunofluorescence (and antigen retrieval) methodology. This was conducted through repeating analysis of mCherry-positive and PCNA-positive Müller glia using an mCherry antibody with all injury paradigms (Supplementary Figure 4). Surprisingly, we indeed found that the majority of Müller glia contained mCherry-positive debris when stained with antibodies against mCherry at 24 hpi (99.4 ± 0.4%, 80.7 ± 7.6%, 73.2 ± 11.2% and 74.8 ± 6.2%), 48 hpi (99.7 ± 0.3%, 88.9 ± 5.9%, 85.1 ± 5.3% and 81.3 ± 3.7%), 72 hpi (99.6 ± 0.4%, 85.5 ± 4.3%, 81.5 ± 5.5% and 81.2 ± 4.4%) and 96 hpi (95.2 ± 2.8%, 73.4 ± 5.0%, 78.8 ± 5.8% and 73.4 ± 4.1%) in our widespread Lws2, smaller Lws2, Sws2 and Xops-ablations, respectively, Supplementary Figures 4A–H. Thus, the proportion of proliferating (activated) and total (including quiescent) Müller glia that contained any mCherry debris was both very high, when including those with substantial debris (analyzed above) as well as those with debris that could only be visualized after antibody staining. Thus, while we did still observe that proliferating Müller glia were significantly more likely to contain mCherry debris in our smaller Lws2 ablation at 72 hpi (p = 0.04), and rod photoreceptor ablation models at 48 hpi (p = 0.03, Supplementary Figures 4I–K), phagocytosis does not seem to be sufficient for Müller glia activation. Correlating this to our scRNAseq data, we found that phagocytosis-associated markers are also upregulated following injury specifically in both quiescent and proliferating Müller glia cell clusters (Figures 4M–O), consistent with this cellular process preceding de-differentiation and cell cycle re-entry. Thus, phagocytosis alone does not seem to be sufficient as a functional process to activate the injury response in glia, but is significantly correlated to Müller glia activation.

Given that we saw spatially distinct activation of Müller glia that could not be explained by site of cell ablation or phagocytic activity alone, we assessed for any heterogeneity in gene expression patterns prior to injury in the mature, quiescent Müller glia population. To achieve this, we performed scRNA-seq of FACS-Müller glia from 6 days post fertilization (dpf) Tg(gfap:GFP) zebrafish larvae. UMAP construction and unsupervised clustering revealed 11 populations of Müller glia and Müller glia-derived cells. Of these, 6 neighboring clusters were very similar (Figure 5A) and showed consistent expression of common genes typical of quiescent Müller glia, such as high expression of glial markers glula, gfap, and rlbp1a (Figure 5C). The remaining 5 clusters did not represent the quiescent populations we were focusing on (marked in gray in Figure 5A). These cell clusters expressed markers of proliferation (pcna) or neural differentiation (crx, neurod4, elavl3, elavl4) or intermediate states (her4.1, igfbp5b; Figure 5B). These clusters likely encompass either young Müller glia exiting a progenitor state or a subset of Müller glia that sporadically enter the cell cycle to generate rod photoreceptors (Bernardos et al., 2007; Ng et al., 2014). Additionally, one cluster displayed high expression of genes enriched for terms suggestive of stress response (Supplementary Figure 5A), such as hsp90aa1.2, hsp70.3 and hsd11b2 (data not shown). From these results, we excluded these clusters for further analysis, and we focused on clusters representing mature quiescent glia (C1-C6). We performed hierarchical clustering analysis of top-ranking expressed genes clusters C1 to C6 and identified key similarities (reflected in their close proximity on the UMAP plot), as well as subtle, but robust differences between these cell clusters (Figure 5A; Supplementary Figure 5B). We found an unexpected diversity of quiescent Müller glia gene expression states within the uninjured retina (Figures 5D,E). While these six clusters were overall grouped together due to their overwhelming similarities as quiescent glia, robust differences were observed in the gene lists that led to unbiased sub clustering.

Cluster C1 displays enrichment for terms including developmental process, neurogenesis and kinase regulator activity, in addition to retinal homeostasis (oxidoreductase activity and sodium ion homeostasis). This cluster expresses high levels of the igf2b (insulin-like growth factor 2b). Among other highly expressed genes in this cluster includes tob1a involved in development of dorsal structures (Xiong et al., 2006), fibroblast growth factor 24 (fgf24), the negative regulators of receptor tyrosine kinase signaling spry2 and spry4 (Felfly and Klein, 2013), the RNA-binding protein rbpms2b, and gap junction protein cx43. Unique to cluster C2 is expression of the glial differentiation gene meteorin (metrn) and fatty acid binding protein 7a (fabp7a). Cluster C2 also expresses the ID signaling gene id1 and thymosin beta-4 (tmsb4x) and enrichment for terms relating to pattern specification process and regulation of neurogenesis. Clusters C3, C4 and C5 were enriched for terms relating to ribosome activity, including ribosome biogenesis, translation initiation and structural constituent of ribosome. Cluster C3 expressed genes eef1a1l1 (translation elongation factor) and rpl3 (60S ribosomal protein L3). Cluster C4 expressed growth factor signaling-associated genes fgf24 and spry4 that were shared with cluster C1, as well as unique genes encoding the amino acid transporter slc1a4, vascular endothelial growth factor Aa (vegfaa) and hypoxia inducible domain family member 1a (higd1a). Additionally, cluster C5 is enriched for terms related synthesis of biological compounds and metabolism, and cells in this cluster express the bHLH transcription factor olig2, igfbp5b (insulin-like growth factor binding protein 5b), her12 and ephrin receptor-binding ligand efnb2a. Despite containing a relatively smaller proportion of cells, the central cluster C6 is distinct in its gene expression profile, displaying enrichment for terms relating to development, regeneration and cellular signaling, expressing genes of the retinoic acid pathway (aldh1a3, rdh10a), and TGFβ superfamily (inhbab, bmpr1ba, rgmd), as well as the caveolar-associated gene cav2, basement membrane protein-encoding gene smoc1 and midkine-a (mdka). In summary, we have identified the following putative quiescent Müller glia subpopulations: Two subpopulations associated with fibroblast growth factor (Fgf) signaling, which can be distinguished based on the presence (C1) or absence (C4) of igf2b, a fabp7a-expressing subtype (C2), Müller glia strongly associated with protein production (C3 and C5) and finally a Müller glia subtype associated with retinoic acid signaling (C6). Whether these clusters represent stable specializations in the terminal differentiation of this cell population to support their diverse roles in maintaining retinal homeostasis, fluid states through which individual glia stochastically transition through, or are indicative of varying developmental ages of Müller glia remains an interesting avenue to explore.

Without spatial or temporal information on individual Müller glia, the scRNA-seq data cannot distinguish whether these differing clusters represent distinct Müller glia subpopulations. If the molecular heterogeneity was linked to specific functions or represent stochastic fluid transition states, we hypothesized to find glia of different gene expression clusters distributed either uniformly or stochastically across all retinal regions. In contrast, we postulated that if the clustering was dependent on the age of Müller glia, these clusters would manifest as a central-to-peripheral distribution pattern in the retina. The retina of fish and amphibians exhibits indeterminate eye growth, with new neurons and glia continuously added from the peripheral ciliary margin zone (CMZ) into the growing central retina (Centanin et al., 2011; Fischer et al., 2013; Centanin et al., 2014). Therefore, older Müller glia are found more centrally than those in peripheral regions. With these possibilities in mind, we assessed the distribution of the distinct clusters both in the developing and adult retina. We performed RNAscope in situ hybridization on candidate genes efnb2a (clusters C2 and C5), fgf24 (clusters C1 and C4), and rdh10a (cluster C6). Unexpectedly, we found the presence of these markers in distinct regions of the retina with efnb2a expression found in the dorsal (21% of retina), fgf24 expression in central (50% of the retina) and rdh10a expression ventral retina (19% of the retina) respectively (Figures 6A–D). These markers co-localised with Müller glia (gfap:EGFP positive) in these regions (Supplementary Figure 6). As expected, PCNA-positive cells were found in the peripheral ciliary margin zone, likely from which igfbp5b-positive, her4.1-positive immature Müller glia are derived from (Figure 6D). Furthermore, this distribution was maintained throughout maturation, with a similar pattern observed at 12 months post-fertilization (Supplementary Figure 7A). Thus, our clustering analysis in quiescent glia has led to the identification of genes that mark regional populations in a pattern that persists with ageing.

We next assessed the pattern of Müller glia recruitment following Lws2 cone photoreceptor ablation using Tg(lws2:nfsb-mCherry) in adulthood. We first characterized the distribution of Lws2 cones and Müller glia in the uninjured adult fish (Figure 7A). Consistent with our data on larvae, both Lws2 cones and Müller glia are evenly distributed, with some additional cells in the central retina (Figure 7B). Following ablation, we used fgf24 and rdh10a as regional markers of spatially segregated Müller glia subpopulations. Consistent with the larval retina (Figure 3), we saw a specific bias of Müller glia in the dorsal and central fgf24+ glia regions, with very little activation and cell cycle entry activity (marked by PCNA immunostaining) in the ventral rdh10a + regions (Figures 7C,D). Thus, even in the context of widespread photoreceptor injury across the retina, Müller glia show region-specific heterogeneity in their capacity to proliferate, and this heterogeneity persists from larval to adult stages. Thus, our clustering analysis has led to the identification of genes that mark regional populations of Müller glia with differing probability to contribute toward the regenerative response.

We found that Clusters C1 to C6 featured GO terms for mature Müller glia functions including neuroprotection, metabolism and ion homeostasis (Newman et al., 1984; Resta et al., 2007; Bringmann and Wiedemann, 2012; Furuya et al., 2012; Reichenbach and Bringmann, 2016). Yet, wanted to determine whether these populations differed in their responsiveness to injury or regenerative potential. Cluster C6 revealed many markers that are known to influence stem cell renewal (Figure 5E), including the bone morphogenetic protein receptor bmpr1ba (Bond et al., 2012; Wang et al., 2014), TGF-β family ligands inhbab and rgmd (Heldin et al., 2009), retinoic acid pathway proteins aldh1a3 and rdh10a (Blum and Begemann, 2012), and the cytokine mdka (Ang et al., 2020). Immunofluorescent labeling of Bmpr1ba/b showed localization in neurons (DAPI-positive, Gfap-negative cells) and Müller glia (DAPI-positive, Gfap-positive cells) in Tg(gfap:GFP) zebrafish (Supplementary Figure 8). We observed labeling in Müller glia in the ventral retina, as well as in the central sector, with minimal labeling present on Müller glia in the dorsal retina. This supports the association of cluster C6 with glia in the ventral domain of the retina. Furthermore, the restriction of these proliferation-associated genes to a subset of Müller glia motivated us to explore potential regional differences in Müller glia proliferation.

To explore this regional heterogeneity in Müller glia activation, we conducted integration and re-clustering analysis of cells collected from our no-ablation and Lws2-ablation scRNA-seq datasets (Figure 8A). There, we identified 4 main quiescent Müller glia populations marked by fgf24, efnb2a and rdh10a expression as genes most uniquely expressed for these quiescent Müller glia clusters, consistent with observations from the uninjured sample alone (Figures 8B,C). Notably, fgf24-positive Müller glia spanned two clusters distinguished by the expression of the insulin growth factor-encoding gene igf2b and the transmembrane glycoprotein-encoding gene cd99. Curiously, cd99 exhibited greater expression following photoreceptor ablation. Furthermore, levels for immune-associated genes also changed following photoreceptor ablation, including the chemokine genes cxcl18b and cxcl14, which were higher in expression following photoreceptor ablation when compared to the uninjured dataset in these quiescent Müller glia clusters (Supplementary Figure 9A), specifically in this grouping of cd99-positive cells (Supplementary Figure 9B). As such, certain Müller glia may be functionally specialized to signal to immune cells.

The range of molecular profiles of quiescent Müller glia remained following neural injury. Gene ontology analysis revealed that the differentially expressed genes in fgf24; igf2b-positive, fgf24; cd99-positive, efnb2a-positive, and rdh10a-positive clusters were enriched for biological processes relating to development, regeneration, or differentiation (Figures 8D,E). As the additional branches of activated Müller glia and differentiating progenitors in our photoreceptor ablation sample must come from the original population of quiescent Müller glia from the uninjured sample, we compared whether all the quiescent Müller glia populations were recruited equally. The greatest differences observed in the proportion of Müller glia between the injury and non-injured control was in the following order: fgf24; igf2b-positive (70.7%), fgf24; cd99-positive (60.6%) clusters, followed by similar differences in rdh10a-positive (51.6%) and efnb2a-positive (47.9%) clusters (Figure 8C). Thus, within the fgf24-positive population of glia, proliferative ability may be influenced by expression of igf2b.