Gli3R-mediated inhibition of hedgehog signaling alters the embryonic transcriptome in zebrafish

TL;DR

This study shows that inhibiting hedgehog signaling in zebrafish embryos affects the expression of genes in various tissues, including a previously unrecognized impact on ionocytes in the skin.

Contribution

The study reveals a novel role for Gli3-mediated hedgehog signaling in regulating ionocyte specification in zebrafish embryos.

Findings

Inhibiting hedgehog signaling reduces expression of genes in lateral floor plate and motor neurons.

Gene set analysis shows downregulation of ionocyte marker genes, suggesting a new role for hedgehog signaling in ionocyte specification.

Reduced expression of foxi3a and ionocyte transporter genes indicates altered development of osmoregulatory cells.

Abstract

Hedgehog signaling is a conserved developmental pathway that patterns diverse tissues during vertebrate embryogenesis. In zebrafish, disruptions to the hedgehog pathway cause well-characterized defects in specific cell types including neurons and glia derived from the ventral neural tube. We inhibited hedgehog signaling by overexpressing the Gli3 repressor ubiquitously and performed bulk RNA sequencing of 30 h postfertilization zebrafish embryos. Consistent with known roles of hedgehog signaling, we observed reduced expression of genes marking lateral floor plate, motor neurons, Kolmer–Agduhr cells, dopaminergic neurons, slow muscle cells, and anterior pituitary. Gene set enrichment analysis using marker genes derived from the Daniocell atlas also revealed downregulation of genes marking H+-ATPase-rich ionocytes, which are located in the embryonic skin and are responsible for osmotic homeostasis. Reduced expression of ionocyte-specific transporter genes and the transcription factor foxi3a suggests that Gli activity may play a previously unrecognized role in the specification of this cell type.

Figures

Morphological and transcriptomic phenotypes of Tg(ubb:zGli3R) embryos. a) Representative transgenic and wildtype animals at 2 dpf and 3 dpf. Scale bar = 200 μm. b) Volcano plot comparing gene expression in 30 hpf Tg(ubb:zGli3R) embryos and wildtype siblings. Genes with −log10(padj) >10 are plotted at 10. Transcripts encoded by chromosome 4 are denoted with X symbols and transcripts encoded by all other chromosomes are denoted with circles. n = 4 replicates of 4 pooled embryos for transgenic and 4 replicates of 4 pooled embryos for wildtype siblings. c) Heat map displaying expression of selected marker genes. From left to right, genes mark hedgehog pathway components, progenitor domains and interneurons in the ventral spinal cord, Kolmer–Agduhr cells, motor neurons, radial glia/astrocytes, dopaminergic neurons, slow muscle cells, and anterior pituitary. d) Cell types identified with GSEA using gene sets from the Daniocell single-cell atlas (Sur et al. 2023). Red bars represent terms with a negative normalized enrichment score and blue bars represent terms with a positive score. e) Bar plots of normalized counts for genes expressed in ionocytes. P values and padj for all genes are available in Supplementary Table 1. f) Representative transgenic and wildtype larvae labeled with Concanavalin A at 3 dpf. Scale bar = 200 μm. g) Quantification of Concanavalin A-positive foci in the yolk sac and yolk sac extension of wild type and transgenic larvae. n = 8 wildtype larvae and 8 transgenic larvae from 1 clutch. Statistics, two-tailed Welch's t-test where P ≤ 0.0001.

Full text

Introduction

The hedgehog signaling pathway plays a critical role in vertebrate patterning (Ingham and McMahon 2001) and is also implicated in regeneration (Imokawa and Yoshizato 1997) and oncogenesis (Taipale et al. 2000). Depending on the context, hedgehog can act as a morphogen (Heemskerk and DiNardo 1994) or mitogen (Wechsler-Reya and Scott 1999) and specifies multiple distinct cell types across organ systems. In the canonical pathway, hedgehog ligands (Sonic, Indian, and Desert hedgehog) (Echelard et al. 1993) bind to the receptor Patched (Marigo et al. 1996; Stone et al. 1996), which relieves inhibition of the transducer Smoothened (Chen et al. 2001; Zhang et al. 2001). Activation of Smoothened shifts the balance between repressor and activator forms of the Gli transcription factors, which in turn modulates the transcription of target genes (Ruiz i Altaba 1999; Wang et al. 2000).

Experimental manipulation of hedgehog in the zebrafish model system has revealed essential roles in development (Hammerschmidt et al. 1996; Huang and Schier 2009; Shen et al. 2013). In the presumptive spinal cord, hedgehog secreted by the notochord and medial floor plate exposes the neighboring lateral floor plate (p3), pMN, and p0-p2 progenitor domains to a gradient of hedgehog levels, which contributes to specification of neural cell types including interneurons, motor neurons, Kolmer–Agduhr cells, oligodendrocytes, and astroglia (reviewed in Danesin and Soula 2017). In addition to gross morphological phenotypes like cyclopia, ventral tail curvature, and rounded somites, inhibition of hedgehog in zebrafish can expand or reduce these ventral spinal cord progenitor domains and their derivative cells (Brand et al. 1996; Karlstrom et al. 1996; Chen et al. 2001).

However, subtle differences in phenotypes result from the specific method used to disrupt hedgehog signaling (i.e. mutations, morpholinos, overexpression of Gli repressors, and chemical inhibitors). These phenotypes may vary due to compensation between paralogs, incomplete or off-target effects of antagonists, maternal contribution of mRNA in homozygous mutants, and noncanonical or hedgehog ligand-independent Gli activity. For example, knockdown of gli1 or gli2b in maternal-zygotic smoothened mutants reduces motor neuron number compared with mutants, suggesting that basal Gli activity contributes to motoneurogenesis (Mich and Chen 2011). Similarly, overexpression of the Gli3 repressor (Gli3R) in embryos treated with the Smoothened antagonist cyclopamine reduces the number of motor neurons compared with smoothened mutants (England et al. 2011). These experiments highlight how changes in the expression and processing of Gli can modify specification of cells in the ventral spinal cord even in the absence of upstream hedgehog signaling.

While previous work has manipulated Gli activity in targeted phenotypic assays, the transcriptome-wide consequences of Gli3R overexpression in zebrafish remain unexplored. To address this gap, we overexpressed a zebrafish Gli3R (zGli3R) transgene (Huang and Schier 2009) ubiquitously and performed bulk RNA sequencing (RNA-seq) of 30 h postfertilization (hpf) zebrafish embryos. Our findings suggest that Gli activity may be involved in the development of ionocytes. We anticipate that our dataset will be most relevant to zebrafish and developmental biologists who are interested in hedgehog signaling and cell-type specification during vertebrate development.

Materials and methods

Zebrafish husbandry and transgenesis

Animal experiments were approved by the UMass Chan Institutional Animal Care and Use Committee (IACUC protocol 202300000053). Transgenic animals were generated in an Ekkwill (EK)-based wildtype strain, and both larvae and adult animals were maintained on a 14 h/10 h light/dark cycle at 28 °C. The zGli3R transgene was subcloned from pT2-hsp-zGli3R-EGFP into Gateway pDONR221 (Invitrogen) and then into a Gateway-compatible Tol2 destination vector containing the ubb promoter (Mosimann et al. 2011) and the myl7:GFP transgenesis marker (Bradford et al. 2022). The annotated plasmid map including full sequence is available in Supplementary File 1. To produce transgenic animals using Tol2, one-cell embryos were injected with 1 nL of 10 ng/µL Tol2 mRNA (Kawakami et al. 2004) and 20 ng/µL plasmid. F0 embryos with expression of the myl7:GFP transgenesis marker were grown to adults.

Sample collection for RNA-seq

A mosaic F0 male was crossed to a wildtype female, and F1 embryos from 1 clutch were collected at 30 hpf for RNA-seq. Transgenic embryos and wildtype control siblings were sorted based on the myl7:GFP transgenesis marker. The insertion site was not mapped, and transgenic embryos may possess multiple insertions of the ubb:zGli3R construct. Four dechorionated embryos were pooled and frozen on dry ice for each biological replicate, and 4 biological replicates were collected per genotype.

RNA-seq

RNA extraction, quality control, library preparation, and sequencing were conducted at Azenta Life Sciences (South Plainfield, NJ, USA) as follows: Total RNA was extracted from frozen tissue samples using Qiagen RNeasy Plus Universal mini kit following manufacturer's instructions (Qiagen, Hilden, Germany). RNA samples were quantified using Qubit 3.0 Fluorometer (ThermoFisher Scientific, Waltham, MA, USA) and RNA integrity was checked with 4200 TapeStation (Agilent Technologies, Palo Alto, CA, USA).

RNA-seq libraries were prepared using the NEBNext Ultra II RNA Library Prep Kit for Illumina using manufacturer's instructions (NEB, Ipswich, MA, USA). Briefly, mRNAs were initially enriched with Oligo d(T) beads. Enriched mRNAs were fragmented for 15 min at 94 °C. First strand and second strand cDNA were subsequently synthesized. cDNA fragments were end repaired and adenylated at 3′ ends, and universal adapters were ligated to cDNA fragments, followed by index addition and library enrichment by PCR with limited cycles. The sequencing library was validated on the Agilent TapeStation (Agilent Technologies, Palo Alto, CA, USA), and quantified by using Qubit 3.0 Fluorometer (Invitrogen, Carlsbad, CA) as well as by quantitative PCR (KAPA Biosystems, Wilmington, MA, USA).

The sequencing libraries were multiplexed and clustered onto a flowcell on the Illumina NovaSeq instrument according to manufacturer's instructions. The samples were sequenced using a 2×150 bp paired end configuration. Image analysis and base calling were conducted by the NovaSeq Control Software. Raw sequence data (.bcl files) generated from Illumina NovaSeq was converted into fastq files and demultiplexed using Illumina bcl2fastq 2.20 software. One mismatch was allowed for index sequence identification.

RNA-seq analysis

Reads were aligned with the STAR aligner (2.7.3a-GCC 6.4.0-2.28) (Dobin et al. 2013) to GRCz11 release 104 using the Zebrafish Transcriptome Annotation version 4.3.2 (Lawson et al. 2020). Uniquely mapped reads ranged from 21.5 to 29.9 M per sample. Transcripts with zero counts in 4 or more of 8 total samples were removed, and raw counts were normalized using rlog counts in DESeq2 (Love et al. 2014). The default DESeq2 method was used for determining significance, and raw and adjusted (Benjamini–Hochberg) P values are in Supplementary Table 1. Sample scripts for differential gene expression analysis of RNA-seq data are available on GitHub: https://github.com/thymelab/BulkRNASeq.

Downstream analysis of bulk RNA-seq data was conducted using R (Version 4.4.2). Gene set enrichment analysis (GSEA) was performed using the “GSEA” function from clusterProfiler (Version 4.14.6) (Wu et al. 2021). For chromosomal location GSEA, the gene set was generated using RefSeq Genes and Gene Predictions from the GRCz11 assembly, and genes were sorted by abs(log2 fold change) before GSEA. For Daniocell GSEA, the gene set was generated as previously described (Dang et al. 2025) and is available in Supplementary Table 4. Gene sets were derived from clusters including neural, spinal cord/glia, eye, ionocyte/mucus, otic/lateral line, hematopoietic/vasculature, and taste/olfactory. Genes were sorted by log2 fold change before Daniocell GSEA. Volcano plot, GSEA bar plot, and individual gene bar plots were created using ggplot2 (Version 3.5.2) (Wickham 2016). Heatmap was generated using the “heatmap.2” function from gplots (Version 3.2.0) (Warnes et al. 2025). Code to generate figures is available in Supplementary File 2.

Live imaging

A mosaic F0 male was crossed to a wildtype female, and F1 embryos from 1 clutch were collected for bright-field imaging and from a second clutch for Concanavalin A staining. For bright-field imaging, 2 days postfertilization (dpf) and 3 dpf animals were anesthetized with MS-222 and live mounted in 2% low-melt agarose (Fisher Bioreagents, BP165-25) made in embryo medium with methylene blue. For Concanavalin A staining, Alexa Fluor 488-conjugated Concanavalin A (Invitrogen, C21401) was reconstituted to 5 mg/mL in 0.1 M sodium bicarbonate then diluted to 500 μg/mL in embryo medium without methylene blue. The 3 dpf larvae were incubated in Concanavalin A for 60 min at room temperature, washed in embryo medium, anesthetized with MS-222, and live mounted in 2% low-melt agarose made in embryo medium. Images were obtained with a Leica M165 FC stereo microscope and K3C color camera. For quantification of Concanavalin A labeling, foci in the yolk sac and yolk sac extension were counted manually using ImageJ (Version 1.53t). Significance was determined using a two-tailed Welch's t-test in GraphPad Prism (Version 10.5.0).

Results

Differential gene expression in Tg(ubb:zGli3R) embryos

We cloned zGli3R downstream of the ubiquitin B (ubb) promoter and coinjected F0 embryos with Tol2 to produce transgenic lines. Transgenic F1 offspring possessed severe morphological defects and died at larval stages (Fig. 1a). We isolated a mosaic F0 parent with high germline transmission and performed bulk RNA-seq on 30 hpf embryos pooled from a single clutch. Analysis of differentially expressed genes identified 888 upregulated genes and 1,018 downregulated genes at a threshold of P adjusted (padj) < 0.05 (Fig. 1b; Supplementary Table 1). Transgene insertion can cause enrichment of differentially expressed genes based on chromosomal location (White et al. 2022), and we observed a significant enrichment of genes on chromosome 4 (Supplementary Table 2).

*Morphological and transcriptomic phenotypes of Tg(ubb:zGli3R) embryos. a) Representative transgenic and wildtype animals at 2 dpf and 3 dpf. Scale bar = 200 μm. b) Volcano plot comparing gene expression in 30 hpf Tg(ubb:zGli3R) embryos and wildtype siblings. Genes with −log10(padj) >10 are plotted at 10. Transcripts encoded by chromosome 4 are denoted with X symbols and transcripts encoded by all other chromosomes are denoted with circles. n = 4 replicates of 4 pooled embryos for transgenic and 4 replicates of 4 pooled embryos for wildtype siblings. c) Heat map displaying expression of selected marker genes. From left to right, genes mark hedgehog pathway components, progenitor domains and interneurons in the ventral spinal cord, Kolmer–Agduhr cells, motor neurons, radial glia/astrocytes, dopaminergic neurons, slow muscle cells, and anterior pituitary. d) Cell types identified with GSEA using gene sets from the Daniocell single-cell atlas (Sur et al. 2023). Red bars represent terms with a negative normalized enrichment score and blue bars represent terms with a positive score. e) Bar plots of normalized counts for genes expressed in ionocytes. P values and padj for all genes are available in Supplementary Table 1. f) Representative transgenic and wildtype larvae labeled with Concanavalin A at 3 dpf. Scale bar = 200 μm. g) Quantification of Concanavalin A-positive foci in the yolk sac and yolk sac extension of wild type and transgenic larvae. n = 8 wildtype larvae and 8 transgenic larvae from 1 clutch. Statistics, two-tailed Welch's t-test where ***P ≤ 0.0001.

We next examined the expression of marker genes for cell types known to involve hedgehog signaling (Fig. 1c; Supplementary Table 1). Progenitor domains in the ventral spinal cord give rise to specific subtypes of neurons and glia, including interneurons, motor neurons, Kolmer–Agduhr cells, oligodendrocyte precursor cells (OPCs), and astrocytes (Reviewed in Danesin and Soula 2017; Cucun et al. 2024). Genes marking most of these cell populations were downregulated or unchanged, but gfap, which is expressed in radial glia and astrocyte-like cells (Marcus and Easter 1995; Agius et al. 2004), was upregulated. Consistent with the literature, marker genes for dopaminergic neurons (Chen et al. 2001; Holzschuh et al. 2003), slow muscle cells (Blagden et al. 1997), and anterior pituitary (adenohypophysis) (Herzog et al. 2003) were also downregulated.

Functional enrichment analysis of cell types using Daniocell single-cell atlas

GSEA using gene sets derived from single-cell expression data can reveal potential changes in cell composition (Capps et al. 2025; Dang et al. 2025; Moyer et al. 2025). GSEA with gene sets derived from the Daniocell single-cell atlas (Sur et al. 2023) highlighted 28 clusters with negative or positive enrichment compared with control (Fig. 1d; Supplementary Table 3). Several of the identified clusters with negative enrichment scores, including floor plate (glia.15 and glia.2), Kolmer–Agduhr cells (glia.30), spinal OPCs (glia.20 and glia.26), dopaminergic neurons of the diencephalon (neur.34), and anterior pituitary (tast.10) are consistent with the known roles of hedgehog during zebrafish development. Other cell types with negative enrichment scores included ventrolateral mesoderm (hema.20), hemangioblasts (hema.23), type III taste buds (tast.2), and multiple subtypes of H+-ATPase-rich ionocytes and ionocyte progenitors (iono.17, iono.15, iono.10, and iono.12). In contrast, olfactory sensory neurons (tast.9), lymphoid hematopoietic stem cells (hema.15), and primitive erythroblasts (hema.32) showed positive enrichment scores. Further investigation of known markers of ionocytes revealed downregulation of the master regulators foxi3a (Hsiao et al. 2007) and gcm2 (Esaki et al. 2009; Shono et al. 2011) as well as genes expressed in H+-ATPase-rich ionocytes (slc4a1b, atp6v0a1a, ca15a, ceacam1, and si:dkey-192d15.2), and Na+-K+-ATPase-rich ionocytes (trpv6 and s100a11) (Fig. 1e) (Sur et al. 2023). Concanavalin A binds H+-ATPase-rich ionocytes but not Na+-K+-ATPase-rich cells (Lin et al. 2006), and live labeling of 3 dpf larvae with Concanavlin A revealed a significant reduction of fluorescent foci in the yolk sac and yolk sac extension of transgenic larvae compared with wildtype siblings (Fig. 1, f and g).

Discussion

Although activation of hedgehog signaling removes Gli3 repressor activity in some contexts (Litingtung and Chiang 2000; Litingtung et al. 2002), gli3 does not appear to be expressed before 12 hpf in zebrafish (England et al. 2011; Sur et al. 2023), and activator forms of Gli may possess basal levels of activity that are hedgehog ligand independent (Karlstrom et al. 2003; Huang and Schier 2009; Mich and Chen 2011). In this study, we explored the transcriptomic consequences of overexpressing zGli3R ubiquitously from the earliest stages of embryonic development. In contrast to experimental manipulations that inhibit hedgehog signaling, such as mutations in smoothened or cyclopamine treatment, ectopic expression of zGli3R has the potential to shift the balance of Gli activator and repressor forms even in the absence of hedgehog signaling.

We first examined marker genes for cell types previously shown to depend on hedgehog signaling (Fig. 1c) and then applied unbiased GSEA with Daniocell terms to discover new cell types that may be regulated by Gli (Fig. 1d). The validity of our experimental approach is supported by misexpression of the hedgehog pathway component ptch2, which is a transcriptional target of hedgehog signaling (Shen et al. 2013), and of known class II genes, which are activated by hedgehog signaling and include nkx6.1, nkx6.2, and olig2 (Fig. 1c) (Briscoe et al. 2000; Guner and Karlstrom 2007; England et al. 2011). Genes marking Kolmer–Agduhr cells (sst1.1, urp1, and urp2), motor neurons (phox2a and isl2a), dopaminergic neurons (th and slc6a3), slow muscle cells (notum2 and tnnt2d), and anterior pituitary (pomcb) were also significantly downregulated (Sur et al. 2023). GSEA using Daniocell gene sets identified several of the same populations, including anterior pituitary (tast.10), dopaminergic neurons (neur.34), and Kolmer–Agduhr cells (glia.30). However, GSEA also nominated cell types with less literature support for hedgehog involvement in zebrafish, including taste buds (tast.2), hemangioblasts (hema.23), and H+-ATPase-rich ionocytes (iono.17, iono.15, iono.10, and iono.12).

Ionocytes are specialized cells found in the embryonic skin that respond to osmotic changes in the aquatic environment and are functionally analogous to amniote renal tubular cells (Reviewed in Hwang and Chou 2013). Progenitors derived from non-neural ectoderm give rise to both ionocytes and keratinocytes (Jänicke et al. 2007). Expression of the transcription factors foxi3a and foxi3b in ionocyte progenitors activates Notch1 in neighboring cells, which in turn suppresses ionocyte cell fate and promotes keratinocyte development (Hsiao et al. 2007; Jänicke et al. 2007; Esaki et al. 2009). The transcription factor gcm2 also plays a role in the specification of ionocytes (Chang et al. 2009), and we observed reduced expression of foxi3a and gcm2 but not of foxi3b (Fig. 1e). Staining with Concanavalin A validated the reduction of H+-ATPase-rich ionocytes predicted by GSEA (Fig. 1, f and g). Given that foxi3a and foxi3b are differentially expressed in H+-ATPase-rich and Na+-K+-ATPase-rich ionocytes (Hsiao et al. 2007), it is possible that activation of hedgehog signaling promotes the differentiation of H+-ATPase-rich ionocytes at the expense of Na+-K+-ATPase-rich ionocytes or other ionocyte subtypes.

Although hedgehog signaling has not been implicated in the specification of zebrafish ionocytes, it does regulate gcm2 expression during development of the amniote parathyroid (Moore-Scott and Manley 2005; Grevellec et al. 2011). Moreover, the Daniocell single-cell atlas shows detectable expression of hedgehog components smo, ptch2, gli1, gli2a, gli2b, and gli3 in non-neural ectoderm (neur.26). However, a recent study that used single-cell RNA-seq to track changes in cell abundance reported that cyclopamine treatment increases the abundance of non-keratinocyte epidermal progenitors and Na+-K+-ATPase-rich ionocytes (Barkan et al. 2025). Future hypothesis-based approaches will determine whether zGli3R overexpression affects ionocyte development via hedgehog signaling or by suppressing basal Gli activity.

A limitation of our study is that overexpression of the zGli3R transgene causes major morphological phenotypes at larval stages (Fig. 1a), which may contribute to changes in gene expression due to developmental delay or allometric shifts in the distribution of tissue types. We attempted to minimize these effects by performing RNA-seq at an earlier time point (30 hpf), and GSEA using a total of 228 gene sets derived from the Daniocell single-cell atlas revealed only 28 gene sets with significant enrichment, indicating that the tissue distributions remain relatively normal. In the future, RNA-seq on dissected tissues rather than whole embryos may help to disentangle the primary and secondary effects of zGli3R overexpression.

Supplementary Material

jkaf306_Supplementary_Data