Differential gene regulation in DAPT-treated Hydra reveals candidate direct Notch signalling targets

Notch signalling facilitates cell fate decisions and pattern formation by inducing terminal differentiation and mediating lateral inhibition, boundary formation and synchronization of developmental processes in animals. Well-studied examples of Notch-regulated processes include the differentiation of the wing margin and the specification of neurons from neuroectoderm in Drosophila embryos and somite formation during vertebrate development (Liao and Oates, 2017; Siebel and Lendahl, 2017). The core components of the Notch pathway include the Notch receptor, the Delta/Serrate/Lag-2 (DSL) ligands and recombining binding protein suppressor of hairless (RBPJ) transcription factors [also called CSL, for CBF1 in mammals, Su(H) in Drosophila and Lag-1 in Caenorhabditis] (Andersson et al., 2011). Both the DSL ligands and Notch receptors are transmembrane proteins, therefore signalling occurs between directly adjacent cells. Interactions between DSL ligands and Notch receptors result in cleavage of the Notch receptor by presenilin followed by nuclear translocation of the intracellular domain of Notch (NICD) (reviewed in Mumm and Kopan, 2000). NICD works as a transcriptional co-activator of CSL factors.

Direct target genes of Notch signalling have been identified previously (reviewed by Giaimo et al., 2021; Wang et al., 2015). Targets of Notch signalling are activated or repressed in different cell types depending on the composition of transcriptional complexes induced by Notch activity and the epigenetic status at the respective loci. A primary and evolutionarily conserved target of Notch is the Hey-Hes family of transcriptional repressors. Other context-dependent direct target genes of Notch signalling that have been identified include Myc, cyclin D1 and MEK5c in tumour cells (reviewed in Borggrefe and Oswald, 2009). In hematopoietic cells, GATA3, the master regulator for T-cell development, and several Hox genes are direct Notch targets (Fang et al., 2007). Genome-wide analysis in Drosophila has shown that genes of the epidermal growth factor receptor pathway are direct targets of Notch signalling and it showed that Notch targeted activators and repressors of certain genes at the same time (Krejci et al., 2009). Notch also induces transcription of its own inhibitors, for example, the small Notch-regulated ankyrin repeat protein NRARP (Jarrett et al., 2019).

To reveal the ancestral core regulatory network directed by the highly conserved Notch signalling pathway, we have focused on a cnidarian, the fresh water polyp Hydra. As a sister to bilaterian animals, cnidarians hold an informative phylogenetic position. Moreover, Hydra provides the unique opportunity to obtain an animal-wide picture of Notch target genes with cell-type resolution due to the recently available single-cell expression map (Siebert et al., 2019).

Hydra polyps have a simple body structure, representing a tube with an oral head structure and an aboral foot. The head consists of the hypostome, with a central mouth opening surrounded by a crest of tentacles. The foot consists of a peduncle, terminating in the basal disc. The body column of the polyp is composed of two epithelial monolayers, termed ectoderm and endoderm, separated by an acellular extracellular matrix, the mesoglea. Ectoderm and endoderm are self-renewing epithelial cell lineages. A third cell lineage, the interstitial cells, resides in interstitial spaces of both epithelia (David and Campbell, 1972; David and Gierer, 1974). It is supported by self-renewing multipotent stem cells, which provide a steady supply of neurons, gland cells and nematocytes. Nematocytes are cnidarian-specific sensory cells, which harbour the nematocyst or cnidocyst used for capturing prey. Epithelial cells divide along the entire body column of the polyps (Holstein et al., 1991) leading to the displacement of cells towards the oral and aboral ends, and into asexually produced buds. Cells arriving at the base of tentacles or at the basal disc cease cell division and induce differentiation into tentacle or basal disc cells. Buds develop into new polyps and are then released from the parent polyp. Sexual reproduction occurs when interstitial lineage derived germ cells develop into egg and sperm cells (Bosch and David, 1986). Ectodermal tentacle cells are battery cells, where each cell harbours several mature nematocytes. Older cells are shed at the tips of the tentacles and the foot. Owing to continual cell divisions, almost all Hydra cells are replaced approximately every 20 days (Otto and Campbell, 1977). Therefore, the homeostatic animal is in a constant state of development requiring the presence of signalling for patterning the body axis and direct cell fate specification (Steele, 2002).

The Hydra Notch pathway components include the receptor HvNotch (Hv for Hydra vulgaris), the ligand HyJagged (Hy for Hydra) and the CSL-homolog, HvSu(H). The basic mechanisms of Notch signalling are conserved in Hydra, including regulated intramembrane proteolysis (RIP) through presenilin, followed by nuclear translocation of the NICD (reviewed in Mumm and Kopan, 2000). Moreover, the promoter of the Hydra HES-family member HyHes can be activated by the HvNotch NICD indicating that HyHes is a direct target of Notch signalling (Käsbauer et al., 2007; Münder et al., 2010; Prexl et al., 2011).

The presenilin inhibitor DAPT efficiently blocks nuclear translocation of NICD and phenocopies Notch loss-of-function mutations in Drosophila and zebrafish (Geling et al., 2002; Micchelli et al., 2003). In the cnidarians Nematostella vectensis and Hydractinia echinata, morpholino-mediated knockdown or CRISPR-Cas-mediated mutagenesis of Notch results in comparable phenotypes to those seen upon DAPT treatment in both organisms, with them displaying defects in nematocyte differentiation and tentacle patterning (Gahan et al., 2017; Marlow et al., 2012; Richards and Rentzsch, 2015).

In Hydra, we have shown that DAPT treatment inhibits NICD translocation, which results in four strong effects. First, DAPT blocks post-mitotic differentiation in the nematoblast and germ cell lineages. Early differentiating nematocytes are genetically specified by the expression of the achaete-scute homolog CnASH (Cn for Cnidarian) (Grens et al., 1995; Lindgens et al., 2004) and morphologically by the presence of a post-Golgi vacuole as an element of capsule development. This cell state disappears in DAPT treated animals. Second, DAPT blocks post-mitotic differentiation of female germ cells causing proliferating germ cell precursors to form tumour-like growths (Alexandrova et al., 2005; Käsbauer et al., 2007). Third, DAPT impairs boundary formation at both parent–bud and body column–tentacle boundaries in such a way that the typically sharp gene expression border margins at these structures become diffuse. At the parent–bud boundary this misexpression of the Hydra FGF-R-homolog kringelchen leads to failure of bud foot formation and detachment (Münder et al., 2010; Sudhop et al., 2004). At the base of tentacles, HyAlx expression, which demarcates the tentacle boundaries (Smith et al., 2000), becomes diffuse and we observe malformations of the head structure (Münder et al., 2013). Fourth, DAPT inhibits Hydra head regeneration and regenerating tissue is not able to re-establish an oral organiser as evidenced by lack of Wnt-3 expression. This leads to failure in developing a properly patterned head with hypostome and evenly spaced tentacles (Münder et al., 2013).

To gain a better understanding of the underlying molecular causes of the Notch inhibition phenotypes, we aimed to identify the transcriptional target genes of Notch signalling. We identified 831 genes that were differentially expressed in response to 48 h of DAPT treatment; 75% of these were downregulated. Single-cell expression data were used to uncover the gene expression patterns at cell-state resolution for the Notch-responsive genes. We found that Notch-responsive genes were expressed in cell states such as differentiating nematocytes and oral cell types, which is consistent with the DAPT-induced phenotypes. To identify potential direct targets of Notch signalling, we also profiled the gene expression changes that occurred immediately after DAPT removal. Investigating the expression dynamics of Notch responsive genes and performing motif enrichment analysis enabled us to predict likely direct targets of Notch signalling in Hydra.

To identify targets of Notch signalling in Hydra, we elucidated transcriptional changes that occur in response to DAPT treatment. We expected that sustained DAPT treatment would result in the misregulation of both direct and indirect Notch targets. We furthermore predicted that direct targets would return to control expression levels after DAPT removal more quickly than indirect targets. We profiled gene expression changes immediately after 48 h of sustained DAPT treatment (0 h time point) to identify all Notch-affected genes. In addition, we profiled gene expression 3 and 6 h after DAPT removal to monitor the recovery of these Notch-affected genes.

To characterise the 3 and 6 h time points after DAPT removal we used reverse transcription quantitative real-time PCR (RT-qPCR) to monitor the expression levels of two genes: (1) HyHES, which is a known direct Notch target (Münder et al., 2010), and (2) CnASH, which is expressed in post-mitotic differentiating nematoblasts (Lindgens et al., 2004), a cell state that is lost in response to DAPT treatment (Käsbauer et al., 2007). Loss of CnASH expression is a secondary (or indirect) effect of Notch inhibition and re-establishment of CnASH expression will only occur after DAPT removal once nematogenesis is restored.

As expected, both HyHES and CnASH were downregulated after 48 h of DAPT treatment. HyHES expression returned to normal levels between 5.5 and 8 h after inhibitor removal, whereas CnASH expression was still downregulated after 24 h (Fig. S1). RNA-seq was therefore performed on tissue samples collected after 48 h DAPT treatment (0 h) and at the 3 h and 6 h time points after DAPT removal, since the 6 h intervals appeared sufficient to distinguish direct from indirect Notch-targets and the 3 h intervals were added to monitor earliest responses in gene expression after resuming NICD activity. The workflow for this experiment is illustrated in Fig. 1.

Genes that were differentially expressed after 48 h of DAPT treatment (time point 0 h) were referred to as Notch-responsive genes (NR genes). Of the 831 NR genes identified, 624 were downregulated (75%) and 207 were upregulated (25%) (Fig. 2A). Clustering NR genes according to their fold changes (Fig. 2B) at the three time points after DAPT removal (0 h, 3 h, and 6 h) revealed 279 genes (201 down, 78 up) genes with re-established expression levels at 3 h, including the confirmed Notch target HyHES. A total of 194 genes (143 down, 51 up) showed re-established expression by 6 h and 313 genes (243 down, 70 up) were still differentially expressed at 6 h, including CnASH. A total of 45 genes, including CnGSC, were differentially expressed at time points 0 h and 6 h, but not at 3 h (Fig. 2A, ‘Other’). In addition, 160 genes were found to be differentially expressed at 3 and/or 6 h, but not at 0 h. These were excluded from further analysis.

Overall, these data reveal changes in gene expression caused by inhibition of the Notch pathway, and uncover which changes are rapidly reversed upon relief of this inhibition. This allowed us to explore the cell type-specific effects of DAPT treatment and identify possible direct targets of Notch signalling [the full list of NR genes is available via Figshare (doi:10.6084/m9.figshare.14681343)].

Next, we elucidated NR gene expression patterns by exploring Hydra single-cell expression data, which were available for 666 (80%) NR genes (Fig. 1). We defined cell state and spatial expression of NR genes on the basis of published cell state annotations (Siebert et al., 2019). Hierarchical cluster analysis revealed groups of genes expressed in specific cell states (Fig. 3): nematoblasts/nematocytes (violet, red, blue and yellow clusters, 315 genes), ectodermal epithelial cells, including battery cells (black cluster, 90 genes), endodermal epithelial cells including tentacle cells (grey cluster, 80 genes), and genes more ubiquitously expressed across cell states (cyan cluster, 79 genes). An additional small subset comprising 102 NR genes included genes with restricted expression in several distinct cell states, such as specific neurons, gland cells, germline cells or ectodermal basal disc cells (green cluster, 102 genes). The majority of these 666 NR genes fell into two broad categories: (1) 47% that were specifically expressed in nematoblasts and nematocysts, and (2) 25% that were specifically expressed in epithelial cells (black and grey cluster, Fig. 3).

In addition, we performed non-negative matrix factorisation (NMF) on the NR gene set as an unbiased means to uncover modules of co-expressed genes (metagenes) and identified 23 metagenes. We then visualized metagene expression on the t-distributed stochastic neighbour embedding (tSNE) representation of selected clusterings from Siebert et al. (2019) (Fig. S2). The NMF analysis identified cell-state-specific modules that were consistent with the hierarchical clustering results (Fig. S2; Fig. 3). Interestingly, a single metagene was found to be expressed in female germ line cells, suggesting Notch function during female gametogenesis (Fig. S2G).

The largest fraction of NR genes have nematoblast- or nematocyte-specific expression. In Hydra, this lineage comprises four types of nematocytes, each of which harbours a single capsule (or nematocyst) of the atrichous isorhiza, holotrichous isorhiza, stenotele or desmoneme type. Nematocytes develop from interstitial stem cells via a proliferative amplification phase with incomplete cytokinesis that results in the formation of nests of 4, 8, 16, and 32 nematoblasts. The cells in these nests undergo a final mitosis and start capsule morphogenesis, a process that can be divided into five stages: (1) formation of a growing capsule primordium from a large post-Golgi vacuole, (2) growing of a tubule elongation of the capsule, (3) invagination of the tubule into the capsule, (4) formation of spines inside the invaginated tubule and (5) hardening of the capsule wall. Nests with mature nematocytes break up, and single nematocytes then get incorporated into the battery cells of the tentacles or into epithelial cells of the body column (David and Gierer, 1974; Engel et al., 2002).

As Notch inhibition by DAPT treatment results in a severe block of nematocyte differentiation, which occurs coincident with or immediately after mitotic exit of differentiating nematoblasts (Käsbauer et al., 2007), we sought to identify the exact differentiation step that was affected. We therefore performed hierarchical clustering for NR genes with expression in nematoblast or nematocyte cell states (Fig. 3, violet, red, blue and yellow clusters) using the Hydra single-cell data (Siebert et al., 2019). The single-cell data revealed four distinct nematocyte differentiation trajectories, and gene expression state changes were identified along these trajectories from stem cells to differentiated nematocytes. Moreover, the single-cell analysis revealed eight distinct nematoblast stages along these four trajectories (nb1 through nb8). Two of these trajectories are annotated as desmoneme and stenotele differentiation, based on marker gene expression (Siebert et al., 2019). In the present study, the clustering of NR genes expressed during nematogenesis revealed that the majority of those genes are strongly expressed in cell states nb4, nb5, nb6, nb7, nb8 and in differentiated nematocytes (nem), with no or much lower expression in the earlier cell states of interstitial stem cells (ISC.nb), nb1, nb2 and nb3 (Fig. 4A). This was also observed by plotting the expression of the NR gene modules onto the single-cell tSNE representation (Fig. S2) which indicates expression in all three nematocyte types including desmonemes (Fig. S2B), stenoteles (Fig. S2C) and isorhizas (Fig. S2D).

To identify the point in the trajectories in which differentiating nematoblasts transition from proliferating to post-mitotic nematoblasts, we looked at the expression profiles of two genes that mark proliferating nematoblasts: (1) proliferating cell nuclear antigen PCNA (t10355aep) and (2) the Zn-finger transcription factor gene zic/odd-paired homolog Hydra-zic (Hyzic) (t13359aep; Lindgens et al., 2004). PCNA expression is seen in states ISC.nb, nb1 and nb2 classifying them as proliferating nematoblasts (Fig. 5A). Nb1 and nb2 express HyZic, confirming that HyZic expression is restricted to proliferative nematoblast states as had been shown before (Lindgens et al., 2004). The absence of PCNA-expression in cell states nb3 and nb4 suggests that these are the earliest post-mitotic nematoblasts producing the nematocyst spine and inner wall protein spinalin (Käsbauer et al., 2007; Koch et al., 1998) and spinalin expression is clearly seen in these cells (Fig. 5A). Expression of the early differentiation marker genes NOWA and CnASH become detectable in differentiation states nb5 through nb8 when nematocyst capsules are formed (Fig. 5A,C).

As further evidence that HyZic and CnASH mark mitotic and post-mitotic stages of nematogenesis respectively, using immunofluorescence we show that CnASH protein is detected in the cytoplasm of nematoblasts that contain vacuoles, which were visualised with anti-NOWA antibody (Engel et al., 2002). By contrast, HyZic protein was detected in the cytoplasm of nematoblasts without visible vacuoles and not found in CnASH positive nematoblasts (Fig. 5B).

Of the 315 NR genes that are expressed in nematoblasts or nematocytes, 314 were downregulated upon Notch-inhibition (Fig. 4A). These downregulated genes include many genes expressed in developing nematocytes such as POU4 (t11335aep), Prdl-a (t21636aep; Gauchat et al., 2004), HyDickkopf 3 (t20111aep; similar to HyDKK3; Fedders et al., 2004), CnASH (t10853aep, Grens et al., 1995; Figs 4B and 5B), NOWA (t15237aep; Engel et al., 2002) and Spinalin (t38568aep; Koch et al., 1998), this gene has now three NCBI entries and encodes a longer protein than initially described [an alignment is available via FigShare (doi:10.6084/m9.figshare.14714169); Fig. 4B]. Using double in situ hybridization to detect POU4 and HyZic transcripts, we found mutually exclusive expression of these two genes in differentiating nematocytes, which demonstrates that POU4 is expressed in post-mitotic nematocytes (Fig. S3A). Using in situ hybridization, we also showed that HyZic-positive nematocytes were not affected by DAPT treatment, whereas CnASH and POU4 expression was lost (Fig. S3B,C). As DAPT treatment causes the disappearance of post-mitotic differentiating nematocytes, which are recognised by their forming of post-Golgi vacuoles (Käsbauer et al., 2007), the seeming loss of POU and CnAsh expression after DAPT treatment (Fig. S3B,C) is caused by a loss of the developing nematocyte cell states expressing these genes.

More than 50% of the nematoblast-specific NR genes remained downregulated and did not recover their normal expression level within 6 h after the Notch-inhibitor was removed (Fig. 4B). This again suggests that downregulation of nematogenesis genes reflects the loss of cell states and is mainly an indirect effect of the block in this process caused by Notch inhibition. By contrast, some nematocyte-specific putative transcription factors did recover their expression levels quickly after DAPT removal. These included a possible class I member of the HMG box superfamily, similar to SOXB3 from Hydractinia echinata [t23172aep, XP_012555836.1; an alignment is available via FigShare (doi:10.6084/m9.figshare.14714169)], a protein with a C-terminal bZIP-Jun-domain [t17964aep; an alignment is available via FigShare (doi:10.6084/m9.figshare.14714169)] and a predicted forkhead box protein I1c-like [t1972aep, an alignment is available via FigShare (doi:10.6084/m9.figshare.14714169)], and (Fig. 4B). Given the rapid recovery of their expression after DAPT removal, these genes may be directly targeted by Notch signalling and possibly play a major role in driving nematogenesis.

In conclusion, our differential gene expression analysis confirms that inhibition of Notch signalling causes a decrease in gene expression in differentiating nematoblast stages, coinciding with the loss of this cell type within 48 h of DAPT treatment (Käsbauer et al., 2007). The gene expression of proliferating nematoblasts remained undisturbed.

About 25% of NR genes for which expression patterns were available had enriched expression in epithelial cells (Fig. 3, black and grey cluster) while the remaining non-nematoblast NR genes showed either sporadic or ubiquitous expression (Fig. 3, cyan and green clusters).

Since previous Notch inhibition studies demonstrated severe malformations of the Hydra head structure (Münder et al., 2013), we aimed to elucidate the effect of DAPT treatment on epithelial body column cells and their derivatives (e.g. specialized head and foot cells). Hierarchical clustering using single-cell data for epithelial cells revealed genes that were expressed (1) in all endodermal and ectodermal epithelial cell types along the oral-aboral axis (Fig. 6, grey cluster), (2) mainly in ectodermal epithelial cells (Fig. 6, cyan cluster) and (3) mainly in endodermal epithelial cells (Fig. 6, green cluster). The majority of these epithelial genes from the grey, cyan and green clusters were upregulated in response to DAPT treatment (Fig. 6). These included 36 genes associated with ER, Golgi and endosomal proteins, such as proteins involved in glycosylation like the oligosaccharyl transferase DAD1(t14233aep|DAD1), a negative regulator of cell death (Roboti and High, 2012). Some were involved in redox regulation and unfolded protein response and some were chaperones [see tables on FigShare; full list of NR genes (doi:10.6084/m9.figshare.14681343) and functional annotation of NR genes: (doi:10.6084/m9.figshare.14681319)]. Moreover, membrane proteins, including 12 G-protein coupled receptors, caspase D (t7281aep; Lasi et al., 2010) and the homolog of the ubiquitin-ligase and Notch-modulator mind bomb were also upregulated (t3105aep). Thus, many of the upregulated epithelial genes seem to be involved in stress responses to DAPT treatment. In contrast, Sp5, which is involved in Hydra head patterning (t29291aep; Vogg et al., 2019) was downregulated (Fig. 7A).

Two sets of genes comprised tentacle genes expressed in endodermal tentacle cells (Fig. 6, red cluster) and in ectodermal battery cells (Fig. 6, black cluster). In both sets, the majority of genes were downregulated upon Notch inhibition (88% of battery cell genes and 71% of endodermal tentacle genes). These included a gene encoding a Na⁺ channel in battery cells (t18364aep; Golubovic et al., 2007) and the collagen gene Hcol1 (t14477; Deutzmann et al., 2000) in endodermal tentacle cells (Fig. 7B). The endodermal matrix metalloprotease gene HMMP was upregulated (t16424aep; Leontovich et al., 2000). Both extracellular matrix genes, HMMP and Hcol1, recovered their expression levels within 3 h (Fig. 7B).

Furthermore, small sets of NR genes were specifically expressed in (1) ectodermal head cells, (2) endodermal head cells, and (3) endodermal foot cells (Fig. 6, included in cyan and green clusters). Another NR gene cluster was expressed in ectodermal basal disc cells (Fig. 6, yellow cluster). These expression patterns also could be seen on tSNE plots after NMF analysis (see Fig. S2I,L,M).

The NR genes expressed in endodermal and ectodermal head cells were largely downregulated and several of these have known functions in head patterning. Of note, HyALX (t16456aep; Smith et al., 2000) is expressed at tentacle boundaries and previous work demonstrated that HvNotch is needed to maintain this expression pattern (Münder et al., 2013). Furthermore, several potential head organizer genes including Wnt7 (t28874aep; Lengfeld et al., 2009), the transcription factor gene TCF (t11826aep; Hobmayer et al., 2000), an Otx-related homeodomain protein (t33622aep), an FGF homolog (t8338aep; annotation confirmed by Monika Hassel, Marburg, Germany) and CnGSC (t1216aep; Broun et al., 1999), were among this downregulated set of head-specific genes (Fig. 7B). Of those, HyALX, CNGSC, Wnt7, FGF and HyTCF recover their normal expression levels within 3 h making these genes candidates for direct targets of Notch signalling. The organizer gene CNGSC was also downregulated and recovered expression after 3 h. However, it was then downregulated again at 6 h. This unusual expression behaviour might indicate the presence of an inhibitory feedback mechanism responding to Notch signalling.

By contrast, the NR genes that were specifically expressed in endodermal foot cells and in ectodermal basal disc cells were largely upregulated in response to Notch inhibition. These include TGF-4 (t25624aep; Watanabe et al., 2014) and a predicted secreted Wnt inhibitor APCDD1 (t11061aep). Thus, Notch inhibition by DAPT resulted in reciprocal regulation of foot and head genes in Hydra, with genes normally expressed at the oral end being downregulated and genes normally expressed at the aboral end being upregulated (Fig. 7B).

These data indicate that Notch signalling regulates gene expression in battery cells and further head patterning genes, including the canonical Wnt signalling components HyWnt7 and HyTCF, whereas the BMP pathway component TGF-4 as well as a secreted Wnt inhibitor, both expressed in the foot, appeared to be negatively regulated by Notch.

The differential gene expression analysis revealed sets of genes that showed shared behaviour after Notch inhibition and re-activation after DAPT removal. This suggests shared regulation, and hence we performed a motif enrichment analysis to uncover respective regulatory elements in genes with similar expression dynamics. This analysis was done for the following gene sets: (1) downregulated only at 0 h, (2) downregulated at 0 and 3 h, (3) downregulated at 0, 3 and 6 h, (4) upregulated only at 0 h, (5) upregulated at 0 and 3 h, and (6) upregulated at 0, 3 and 6 h. Regions of open chromatin, as identified by previously published ATAC-seq data, within 5 kb upstream of each gene were considered in the enrichment analysis (see Fig. 8 and Materials and Methods for details) (Siebert et al., 2019).

The group of genes that were downregulated in response to DAPT treatment and then recovered normal expression by 3 h are the best candidates for being direct targets of Notch signalling. If genes are direct targets of Notch signalling, we would expect to find RBPJ-binding sites (Kopan and Ilagan, 2009). In line with our prediction, the RBPJ motif (Bailey and Posakony, 1995) was enriched in NR genes of this group (Table S1; Fig. 8B). Among the 21 genes with RBPJ-binding sites in their regulatory region, HyAlx (t16456aep) and Sp5 (t29291aep), with six RBPJ motifs, are the top candidates for direct targets of Notch signalling (Fig. 8D,E; Table S1). We also identified the transcription factors pituitary homeobox 1-like (specifically expressed in head cells of the endoderm, t5275aep) and a homeobox protein of the OTX-family (t33622aep). Putative RBPJ motifs were additionally present in genes encoding potential membrane or extracellular proteins, including a foot-specific secreted frizzled-related protein, a potential regulator of Wnt signalling (t15331aep, annotation provided by Bert Hobmayer, Innsbruck, Austria, personal communication).

In addition to the RBPJ-binding site, we found enrichment of further transcription factor-binding motifs belonging to 10 transcription factor families. Homeobox transcription factors were the most abundant motifs identified and several different HMG, forkhead and bHLH motifs were also found (Table S1). Interestingly, this corresponds with the downregulation of transcription factors that potentially bind to these domains, for example, HyHES (bHLH), Jun (bZIP), FoxP1 (Forkhead), HyAlx, (homeobox, t16456aep), OTX-related (homeobox, t33622aep) and PITX-related factors (homeobox, t5275aep) and three SOX-related proteins (HMG-boxes, t23837aep, t23172aep, t5528aep) [see Tables S1, S2, S3; an alignment and phylogeny is available via FigShare (doi:10.6084/m9.figshare.14714169)]. Therefore, these data reveal the possible components of a gene regulatory network influenced by Notch signalling.

For the group of genes downregulated at all three time points, six enriched motifs were found, most notably POU- and PAX-binding motifs (Fig. 8C; Table S1). The POU-gene has previously been implicated in nematocyte differentiation and was found enriched in genes expressed at late stages of nematogenesis (Siebert et al., 2019). HyPOU4TF-2 like was downregulated by Notch inhibition at 0 and 3 h (Table S2). The three predicted Hydra-Pax-genes [t9974aep, t6559aep and t11467aep, see table (doi:10.6084/m9.figshare.14681343) and alignment and phylogeny (doi:10.6084/m9.figshare.14714169) on FigShare] were not amongst the NR genes.

In the group of genes that were upregulated at 0 h, but recovered by 3 h, only the IRF9-binding motif (interferon regulatory factor) was found enriched. For genes that were upregulated at 0 and 3 h, but recovered expression at 6 h, many bZIP-factor-binding motifs were enriched.

As expected, this analysis uncovered RBPJ-binding sites in several NR genes. Of those, and in accordance with their quick recovery after DAPT removal, HyAlx and HySp5, each with six putative RBPJ-sites appeared to be the strongest candidates for direct transcriptional targets of Notch signalling, followed by a negative regulator of Wnt signalling (secreted frizzled-related protein) and the potential transcriptional repressor MAD, which has a hypothetical function in regulating proliferation in epithelia cells.

Inhibiting Notch signalling induces a block in nematocyte differentiation and disrupts head patterning in Hydra (Münder et al., 2013). Comparable Notch effects have also been described in two other cnidarian model organisms, namely Nematostella and Hydractinia (Gahan et al., 2017; Layden and Martindale, 2014; Marlow et al., 2012; Richards and Rentzsch, 2015). In this study, we have identified Notch-regulated (NR) genes by analysing RNA-seq data obtained at different timepoints after treatment of Hydra polyps with DAPT. Exploration of the Hydra single-cell gene expression atlas (Siebert et al., 2019) revealed sets of genes that were expressed in cell states consistent with observed inhibition phenotypes. Moreover, in many NR genes, we detected binding sites for DNA-binding protein RBPJ – the principal effector of Notch signalling.

Unexpectedly, we detected upregulation of genes encoding heat shock proteins and proteins involved in apoptosis (see table ‘Functional annotation of NR genes’ at doi:10.6084/m9.figshare.14681319), hinting at a stress response of the animals to the treatment, which we can attribute to DAPT, as effects of the solvent DMSO should be hidden in our experimental design. Promoter analysis of those upregulated genes revealed enrichment of a Trp-cluster motif (IRF9, Table S1). This motif is targeted by the interferon regulatory factor as part of a stress and anti-viral defence pathway in mammals (Jefferies, 2019). Its occurrence and function in Hydra genes should be elucidated in the future.

Strikingly, almost half of the Notch-responsive genes were expressed in cells of the nematocyte lineage; 99% of those were downregulated and expressed in post-mitotic nematoblast stages. This reflected Notch regulation of a gene module specific to post-mitotic nematoblasts that are in the process of capsule formation. In this module, we find genes that have been previously shown to encode structural capsule proteins, such as minicollagens (Engel et al., 2001), spinalin (Koch et al., 1998), NOWA (Engel et al., 2002), nematogalectin (Zhang et al., 2019), N-col15 (Adamczyk et al., 2008), nematocilin (Hwang et al., 2008) and others. Moreover, we find transcription factors like HyPOU (Siebert et al., 2019) and CnASH (Grens et al., 1995). This is consistent with Notch-inhibition blocking differentiation and the initiation of a transcription programme for capsule formation. We propose that HyPOU is involved in executing this programme, but that it is not a direct Notch target since expression is not re-established within 3 h. In accordance with this hypothesis, motif enrichment analysis identified POU4F DNA-binding motifs in the putative regulatory region of genes that did not recover from DAPT within 6 h. Similarly, CnASH expression does not recover within 6 h, confirming it as an indirect or secondary target (Fig. S1). In contrast, the potential Hydra Jun gene, which encodes a C-terminal bZIP-Jun-domain protein, did recover expression levels fast (Table S2) whereas bZIP-binding motifs appeared to be enriched in promoters of downregulated genes that remained downregulated for 6 h post treatment (Table S1). Recent work in Nematostella had indicated the Jun-homolog Cnido-Jun is involved in driving nematogenesis, as knockdown of Cnido-Jun resulted in loss of expression of NvNcol3 (Nv is Nematostella vectensis) and defects in nematocyte morphology (Sunagar et al., 2018). The mammalian transcription factor JUN is part of the activator-protein 1 (AP-1) transcription factor complex, which responds to numerous extracellular signals including MAP-kinase and cytokine signalling as a reaction to environmental signals. AP-1 is also known as a driver of differentiation in the immune system (reviewed in Katagiri et al., 2021). It is tempting to speculate that Notch controls AP-1-like transcriptional regulation in Hydra (Fig. S4B). This could be a mediator for adjusting nematogenesis in adult animals to nutrient-dependent requirements for nematocyte production from precursors. In starving animals, for instance, mature nematocytes are not used, thus turnover is low and this governs replenishment (Yaross and Bode, 1978).

Genes that are expressed in nematoblast precursors, including HyZIC (Lindgens et al., 2004) and the Hydra Pax-2A homolog (t9974aep) were not affected by DAPT. This is also true for the Hydra homolog of Myc, Hymyc1 (Hartl et al., 2010), the human homolog of which is a Notch target gene in mammals (Giaimo et al., 2021).

Genes that are potentially directly targeted by the NICD would not only be expected to recover their expression level quickly when DAPT treatment is removed and NICD is allowed to enter the nucleus but would be expected to also contain RBDJ sites in their promoter regions. Such binding motifs have been detected in a number of nematocyte-specific genes with unknown function (Table S3). These genes do not encode transcription factors, suggesting that the NICD directly activates the nematocyte differentiation gene complex. Future studies will reveal their role during nematoblast differentiation and also whether they can account for the missing differentiation cue that is directly blocked with DAPT.

However, as an alternative explanation, failure to carry out the nematoblast differentiation programme in our experiments could be caused by missing patterning signals from the Hydra head. This hypothesis is suggested because of the strong head phenotypes that we had previously described after DAPT inhibition. The first observable phenotype after 48 h of Notch inhibition was a substantial shortening of the tentacles. Moreover, transplantation experiments with GFP-labelled body column tissue indicated that, during the time of Notch inhibition, cells did not cross the boundary between body column and tentacles (Münder et al., 2013). A ‘neck-like’ structure appeared underneath the tentacle zone, where cells had ceased proliferating, but also did not differentiate into battery cells. In this study, we reveal a cluster of downregulated head-specific genes among the genes that are dysregulated in response to Notch inhibition.

Of particular interest is the aristaless-related gene HyALX, which has six potential RBPJ sites in its putative regulatory region. Our study strongly suggests that Notch signalling directly activates HyALX expression. HyALX has previously been proposed to instruct the specification of tentacle tissue (Smith et al., 2000), and we suggest that HyAlx could play this key role in directing tentacle fate by activating genes with homeobox transcription factor-binding motifs. In support of this, we found that the homeobox motif was enriched in the NR genes downregulated by DAPT that recover their expression quickly after DAPT removal (Table S1). Another potential Notch target with five RBPJ sites in its promoter region is the Max-dimerisation domain (MAD)-encoding MAD gene. This is part of the MYC/Max/MAD network of transcription factors that are involved in the regulation of cell proliferation. MAD forms heterodimers with the bHLH transcription factor MAX, which often mediates repression of proliferative gene activity (Lüscher, 2012). This could also play a role at the tentacle boundary where proliferation of epithelial cells is stopped when they pass into tentacles and become battery cells.

HyALX is expressed in evenly spaced rings at the body column–tentacle boundaries. After release of DAPT inhibition it rapidly recovers expression levels, yet it does not recover a regular expression pattern but becomes expressed in irregular rings, and in extreme cases is in only one ring surrounding the whole animal (Münder et al., 2013). Assuming that NICD acts as a direct activator of HyALX, this indicates that Notch signalling is resumed in the wrong places. Therefore, a feedback mechanism can be suggested, where the Notch signalling pattern depends on the head organizer, which in turn is co-instructed by Notch signalling.

The potential Hydra head organizer gene CnGSC (Broun et al., 1999) is downregulated by Notch inhibition and does not recover its activity after 6 h. Furthermore, Wnt7, TCF and Sp5, genes implicated in the canonical Wnt signalling pathway (Broun et al., 1999; Lengfeld et al., 2009; Vogg et al., 2019), are also downregulated. In contrast, we found a small cluster with foot and peduncle genes that were upregulated, including the BMP-pathway gene TGF-4 (Watanabe et al., 2014). Together, these data may guide the uncovering of molecular pathways responsible for the irregularly shaped heads that develop in polyps after a 48 h period of DAPT treatment (Münder et al., 2013). They also confirm a role of Notch signalling in establishing and maintaining the Hydra head organizer, which was previously discovered in transplantation experiments, where the organizer capacity of regenerating Hydra head tissue had been inhibited by DAPT (Münder et al., 2013).

The role of Notch signalling for the maintenance of tentacle boundaries can be explained when HyAlx and Sp5 are direct targets for activation by NICD. Expression of HyAlx would then always be maintained by a strong Notch signal at the tentacle boundary. Sp5 would also be expressed in response to this Notch signal to block the activity and expression of canonical Wnts at the boundary. All canonical Wnt genes are expressed in the head outside the tentacle zone. Non-canonical Wnt signalling, on the other hand, does not appear to be affected by NICD and therefore the PCP pathway is active at the boundary and guides movements of cells into tentacles (Fig. S4A).

This study suggests target genes of Notch signalling in Hydra, and provides a resource for the investigation of molecular mechanisms by which HvNotch affects patterning, maintenance of the head organizer and post-mitotic nematocyte differentiation. The expression of the only direct HvNotch target gene, for which experimental evidence is available, HyHes, was also found among NR genes, which quickly recovered original expression levels after DAPT removal. We have identified HyAlx and Sp5 as prime candidates for further direct HvNotch targets involved in head patterning due to their quick recovery after DAPT relief and the presence of RBPJ sites in their promoter regions. A candidate for a direct HvNotch target gene expressed in differentiating nematoblast states and quickly recovering from DAPT treatment is HyJun. As a component of the AP1 transcription complex, it might synchronise nematocyte differentiation as a response to the demand for mature nematocytes depending on usage. Moreover, the impact of HvNotch on regulation of this differentiation step might be conveyed by inducing expression of genes encoding proteins other than transcription factors, for instance genes that are required to form the post-Golgi vacuole. It has to be considered that many genes with as-yet-unknown functions are amongst potential direct Notch targets in nematoblasts.

Animals of the strain Hydra vulgaris (Basel) were grown in Hydra medium (0.1 mM KCl, 1 mM NaCl, 0.1 mM MgSO₄, 1 mM Tris and 1 mM CaCl₂) at 18°C and fed regularly with freshly hatched Artemia nauplii.

Regularly fed animals were starved for 24 h and incubated in either 20 μM DAPT with 1% DMSO in Hydra medium or in only 1% DMSO in Hydra medium (control sample) for 48 h. DAPT and DMSO were renewed every 12 h. Animals were collected, and total RNA was isolated at three different time points: directly at the end of 48 h (0 h), 3 h after DAPT removal (3 h) and 6 h after DAPT removal (6 h). After 48 h incubation, DAPT was removed and replaced with 1% DMSO in Hydra medium for the samples 3 h and 6 h. About 25 animals were collected per sample. Six biological replicates were analysed for RNA-seq respectively, two biological replicates (with three technical replicates each) were used for qPCR (Fig. S1).

For each sample, total RNA was extracted from 25 whole animals using the RNeasy Mini kit Plus (Qiagen) according to the manufacturer's protocol at time points of 0 h (48 h of DAPT treatment), 2 h, 5.5 h, 8 h, 10 h and 24 h after DAPT removal, for both DAPT-treated and control (1% DMSO only) animals. RNA quality and quantity were assessed using an Agilent Bioanalyzer. RNA with a RIN value of at least 8 was used for cDNA synthesis using the iScript cDNA synthesis kit (BioRad) according to manufacturer's protocol. A non-RNA and non-reverse transcriptase control were included.

Primers for qPCR were designed using the NCBI primer designing tool (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) (HyHES, fw 5′-CCCACCACCTAGTCCTTCTTC-3′, rev 5′-TTCTGCTTGGGCAAGTTTGG-3′; CnASH, fw 5′-AGACGTTCTAGTCATAGTGTTGTC-3′, rev 5′-AGCCATCATTGACCTGCTTTAC-3′) and tested to ensure they amplified the correct fragment from cDNA by gel electrophoresis. Gene-specific primer pairs that yielded one melt peak and a linear standard curve were used for qPCR quantification.

cDNA was diluted 1:25 to ensure the used concentration was within the standard curve of the primers. qPCR with SYBR green detection was performed using an CFX96 Touch Real-Time PCR Detection System (BioRad). Each measurement was performed in three technical replicates. A no template control (NTC) was included. The genes RPL13, EF1α and PPIB served as housekeeping genes and their geomean was used for normalization. The samples were analysed on a 96-well plate in a CFX96 Touch Real-Time PCR Detection System from BioRadRelative expression was calculated as 2^-(dCt(test sample)−dCt(reference sample)). The error bars represent the s.e.m. (Fig. S1).

Animals were briefly (1–2 min) relaxed in 2% urethane in Hydra medium and fixed immediately after in 2% paraformaldehyde in Hydra medium for 1 h. Animals were washed with PBS, permeabilized with 0.5% Triton-X-100/PBS (15 min) and blocked with 0.1% Triton-X-100/1% BSA/PBS (20 min). Primary antibodies were applied overnight at 4°C. After a PBS-wash, animals were incubated with secondary antibodies (2 h), washed again with PBS, counterstained for DNA with DAPI (Sigma, 1 µg/ml) and mounted on slides in Vectashield mounting medium (Axxora).

RNA in situ hybridization experiments were carried out as previously described (Grens et al., 1995) using digoxigenin labelled RNA probes (Roche) and substrates NBT/BCIP or BM Purple (Roche).

RNA-seq libraries were prepared for six biological replicates for each experimental condition. cDNA libraries were synthesized from total RNA using the strand-specific SENSE mRNA-Seq Library Prep Kit V2 for Illumina (Lexogen) and the Purification Module with Magnetic Beads (Lexogen). The samples were multiplexed and sequenced on three lanes on Illumina Hiseq2000 with a 100 bp paired end sequencing strategy. Downstream analyses were performed using the Galaxy platform and within R [RStudio Team (2016); version 1.1.463; RCode provided at doi:10.6084/m9.figshare.14681310]. Illumina adapters and polyA sequences were trimmed and splice leader sequences (Stover and Steele, 2001) were removed from both forward and reverse reads. The tool ‘fastqfilter’ was used to ensure the paired nature of the filtered dataset, to filter out reads with a quality score lower than 20 and to exclude reads with a read length shorter than 30 bp. Reads that contained ‘N’s were also removed from the dataset.

All forward and all reverse reads from all sequencing libraries were concatenated. The two resulting files were then used as input to the Trinity (version 2.8.4; Grabherr et al., 2011) de novo transcriptome assembler. The assembly was run with the following parameters: -strand-specific library, in silico normalisation, -min_contig_length 300, -min_kmer_cov 1, no genome guided mode and no Jaccard Clip options. The resulting reference transcriptome resulted in 62,419 transcripts (43,481 genes) with an average transcript length of 1008b and a median length of 588 bp. A total of 10% of the genes have an average length of 5017 bp and 50% of them are 1495 bp in average. The average GC content of all genes was 34.7%. Transcripts that belonged to the same gene were joined to form SuperTranscripts (tool ‘Generate SuperTranscripts from a Trinity assembly’; Galaxy version 2.8.4), which were then used for local Blast search. These were treated as genes models in downstream analyses (see below).

The processed reads of the 36 RNA-seq libraries (timepoints 0, 3 and 6 h, DAPT and control samples, six replicates) were separately mapped to the de novo assembled transcriptome reference, within the Galaxy platform. The reads were mapped as strand-specific and with a maximum insert size of 800. RSEM (Li and Dewey, 2011) – with Bowtie2 (Langmead and Salzberg, 2012; Langmead et al., 2009) as alignment method – was used as the abundance estimation method. The overall alignment rate was >95% for all samples except for one control sample, which had an alignment rate of only 88%. The transcript alignment files of all samples and the gene_to_transcript_map were used as input to generate an expression matrix for all 43,481 assembled genes.

The raw counts were used as input for differential expression (DE) analysis by DESeq2 (version 1.18.1). Genes that were not detected in all 36 samples were excluded from this analysis. DE analysis was performed for each time point separately by comparing the DAPT treatment replicates with those from the control animals (0 h DAPT versus 0 h DMSO, 3 h DAPT versus 3 h DMSO and 6 h DAPT versus 6 h DMSO). Differentially expressed genes at time point 0 h were selected according to their P-adjusted value [Padj(FDR) <0.01]. We refer to this gene set as Notch-regulated genes (NR genes). For each of these NR genes, we investigated whether DE was also identified at time points 3 h and 6 h, thereby applying the same cutoff for DE [Padj(FDR) <0.01].

Several blast searches were performed to annotate NR genes. The NCBI Hydra vulgaris protein database (on 2020.02.24) was interrogated using blastx. Sequences with no blast hit or a blast hit with an E-value >10×10⁻¹⁰⁰ were blasted manually. Three types of manual blast searches were performed, NCBI blastn and blastx and NCBI smartBLAST (full list of NR genes on Figshare; doi:10.6084/m9.figshare.14681343). Sequences for which no blast hits were found in either blast search were denoted with ‘no blast hit’. This was also the case for sequences, for which a blast hit was found but with an E-value >10×10⁻²⁰. For the genes that were blasted manually, the NCBI description and accession number was replaced by those of the blast hit with the highest E-value and query cover (for example, if the manual blastn search yielded a better hit than the local blast to the NCBI protein database). The PubMed accession number was added for known Hydra genes. Uniprot was used to search for information about the function and the compartment of the identified sequences, these were denoted as ‘unclear’ in cases it was unclear or unknown. Multiple alignments were performed for genes with a similar TrinityID and genes with similar/same NCBI description.

To make use of the available Hydra single-cell data, we first identified NR genes within the single-cell transcriptome reference using blastn (Siebert et al., 2019) (Transcriptome Shotgun Assembly project; GHHG01000000; see https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM4009036). Duplicated hits were removed by keeping the alignments with the highest blast score. Existing Seurat data objects were used to retrieve expression data and cell state annotations (Siebert et al., 2019; and see related data on Dryad at https://doi.org/10.5061/dryad.v5r6077). For hierarchical clustering approaches, average cluster expression was calculated for each cell state (Seurat_2.3.4::AverageExpression). Seurat objectes were then subsetted to the NR gene set and expression was scaled from 0 to 1. Hierarchical clustering was performed using functions stats::dist(‘euclidian’) and stats::hclust(‘ward.D’). A heatmap (gplots_3.0.0::heatmap.2) with the scaled average expression was generated.

Normalized expression information was extracted from the whole transcriptome Seurat object for each DE gene with an AEP reference and used for non-negative matrix factorization (NMF) analysis. This analysis was performed as described by Siebert and colleagues (Siebert et al., 2019).

To identify putative promoter regions of NR genes, we used a previously published ATAC-seq dataset generated from whole wild-type Hydra (Siebert et al., 2019) to locate regions of accessible chromatin (i.e. peaks) within 5 kb upstream of NR gene transcription start sites. We then grouped these NR promoter regions based on the expression dynamics of their putative target genes in our DAPT-treated RNA-seq time course. A total of six sets of NR genes were considered for downstream motif enrichment analyses: (1) genes that were downregulated but recovered by 3 h post-treatment, (2) genes that were downregulated but recovered by 6 h post treatment, (3) genes that were downregulated and remained downregulated at 6 h post treatment, (4) genes that were upregulated but recovered by 3 h post treatment, (5) genes that were upregulated but recovered by 6 h post-treatment, and (6) genes that were upregulated and remained upregulated at 6 h post treatment.

For our motif enrichment analysis, we used a curated set of known transcription factor binding motifs provided by the JASPAR database (Fornes et al., 2020). Specifically, we used position weight matrices from the non-redundant vertebrate, insect, nematode and urochordate JASPAR datasets. JASPAR-formatted position weight matrices were converted to HOMER-formatted motifs using the HOMER parseJasparMatrix function. HOMER-formatted motifs require the specification of a score threshold that is used for identifying true motif hits in a query sequence. No such score threshold is included in JASPAR-formatted motifs, so we manually set the threshold to be 40% of the maximum possible score (i.e. the score that would be received by a sequence that perfectly matches the canonical binding sequence) for each motif.

We then used this custom set of HOMER motifs to identify transcription factor-binding motifs that were significantly enriched in each of the six abovementioned NR peak sets. We did this by comparing the NR peak sets to non-NR peaks using a binomial enrichment test as implemented in the HOMER findMotifsGenome function. Motif enrichment results were then filtered using a false discovery rate threshold of ≤0.05.

We found that our raw HOMER results included numerous enriched motifs with highly similar sequences. To simplify these results, we sought to identify and remove redundant motifs from the results tables. To accomplish this, we first generated a matrix of pairwise similarity scores for all motifs in our custom motif set using the HOMER compareMotifs function. These similarity scores were then used to perform hierarchical clustering to identify groups of highly similar motifs. We then reduced the redundancy of our enrichment results by including only the most significantly enriched motif from each motif cluster in the final results table.

To identify putative RBPJ-binding sites in NR promoter regions, we used the HOMER scanMotifGenomeWide function to find sequences that matched the RBPJ and Su(H) binding motifs (JASPAR matrix IDs MA1116.1 and MA0085.1, respectively). In addition, we also made use of a custom Su(H) motif based on a previously reported description of the Su(H) consensus binding site (Bailey and Posakony, 1995). The custom HOMER Su(H) motif was generated using the HOMER seq2profile function; the score threshold was set to be 40% of the maximum possible score.

Plots of ATAC-seq read density and predicted RBPJ-binding sites were generated using the R Gviz package (Hahne and Ivanek, 2016). ATAC-seq reads from individual biological replicates were pooled before generating read density plots.

The Galaxy platform of Blum's group at the Gencentre Munich was used for read processing.

The peer review history is available online at https://journals.biologists.com/jcs/article-lookup/doi/10.1242/jcs.258768