Active Notch signaling is required for arm regeneration in a brittle star.

Cell signaling pathways play key roles in coordinating cellular events in development. The Notch signaling pathway is highly conserved across all multicellular animals and is known to coordinate a multitude of diverse cellular events, including proliferation, differentiation, fate specification, and cell death. Specific functions of the pathway are, however, highly context-dependent and are not well characterized in post-traumatic regeneration. Here, we use a small-molecule inhibitor of the pathway (DAPT) to demonstrate that Notch signaling is required for proper arm regeneration in the brittle star Ophioderma brevispina, a highly regenerative member of the phylum Echinodermata. We also employ a transcriptome-wide gene expression analysis (RNA-seq) to characterize the downstream genes controlled by the Notch pathway in the brittle star regeneration. We demonstrate that arm regeneration involves an extensive cross-talk between the Notch pathway and other cell signaling pathways. In the regrowing arm, Notch regulates the composition of the extracellular matrix, cell migration, proliferation, and apoptosis, as well as components of the innate immune response. We also show for the first time that Notch signaling regulates the activity of several transposable elements. Our data also suggests that one of the possible mechanisms through which Notch sustains its activity in the regenerating tissues is via suppression of Neuralized1.


Introduction
In metazoans, a surprisingly small number of signaling pathways are required to control animal development [1]. Knowledge of how these pathways work and interact in different contexts is key to gaining a mechanistic understanding of processes, including development, growth and post-traumatic regeneration in adults. Post-traumatic regeneration often requires dynamic changes in the balance between undifferentiated progenitors and specialized differentiated cells. In response to trauma, the cells of the damaged adult tissue have to be activated and instructed to re-enter the cell cycle, engage in migratory behavior with coordinated spatial rearrangements, and eventually differentiate into specialized cell types of a new body part. As

PLOS ONE
PLOS ONE | https://doi.org/10.1371/journal.pone.0232981 May 12, 2020 1 / 19 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 [18,19]. The only functional study of the Notch signaling pathway in the context of adult echinoderm regeneration was performed in the sea urchin Lytechinus variegatus [20]. This work demonstrated the requirement of the functional Notch signaling for the proper outgrowth of amputated external appendages, such as spines and podia. The cellular and molecular processes regulated by Notch signaling in echinoderm regeneration remain unknown. In addition, echinoid spines and podia are relatively simple structures. The role of the Notch signaling pathway in the regeneration of more complex organ systems and appendages in adult echinoderms has yet to be addressed. Our aim in this study is to establish the functional role of the Notch signaling in arm regeneration in the brittle star Ophioderma brevispina (Say, 1825) and identify the target genes that are regulated by the pathway. Brittle star arms are segmented body appendages with complex internal anatomy. Each brittle star arm contains a calcareous endoskeleton composed of serial vertebral ossicles and several peripheral elements. Associated with the skeleton, the brittle star arm has a system of muscles and ligaments, two systems of coelomic canals, and a complex nervous system including a radial nerve and numerous peripheral nerves [21,22]. Brittle stars have emerged as important models in regenerative biology. They have been used in studies of skeletogenesis and biomineralization [23,24], morphogenesis, and regulation of growth and differentiation [25]. Here, we show that exposing regenerating brittle stars to the Notch pathway antagonist DAPT significantly impairs regeneration. We also identified genes regulated, directly or indirectly, by the pathway by performing a transcriptome-wide gene expression analysis (RNA-Seq). We show that Notch affects a multitude of biological processes involved in arm regeneration, including the extracellular matrix composition and remodeling, cell proliferation, death and migration, activity of mobile genetic elements, and the innate immune response. Our data also indicates an extensive cross-talk between Notch signaling and other key cell signaling pathways, such as Wnt, TGF-β, Toll, JNK, and others, suggesting that regeneration depends on a complex regulatory environment.

Results
The Notch signaling pathway in the brittle star As expected, genes corresponding to the major components of the Notch pathway (i.e., the Notch receptor, the Delta and Serrate ligands, the transcriptional regulator RBPJ, two Notch target genes of the Hes family, and the Notch signaling modulator), as well as a Notchless protein homolog that plays a role in the Notch signaling pathway, were identified in the complete transcriptome of O. brevispina. Details are provided on a general Feature Format version 3 (GFF3) file, available in the supplements to this article (S1 File).

Pharmacological inhibition of the Notch signaling slows down arm regeneration
On day 14 post-autotomy, after continuous exposure of the treatment and control cohorts to the DAPT reagent (3 μM) or DMSO (vehicle), respectively, two arms from each individual were processed for histological analysis. Six individuals per experimental group were used. As compared to the control individuals, the length of the regenerated arm (outgrowth) in the DAPT-treated cohort was on average 37% shorter (967 μm vs 605 μm, T-test p-value = 0.002) (Figs 1 and 2A). Likewise, the newly regenerated arm portion in the DAPT-treated animals had, on average, two fewer segments then the arms in the control group (5.5 vs 7.5, T-test pvalue = 0.047) (Fig 2B). Even though the new arms in the DAPT-treated individuals were smaller and less segmented, they still contained all major radial organs, including the radial A-A 0 @: Aboral view of regenerating arms from four representative control animals exposed to DMSO (vehicle). Two arms from each animal are shown. B-B 0 @: Aboral view of regenerating arms from four representative animals exposed to 3 μM DAPT. Two arms from each animal are shown. C and D: Representative sagittal sections through the regenerating arm of a control (DMSO-treated) individual (C) and a DAPT-treated individual (D). Hematoxylin and eosin staining. Arrowheads show the arm outgrowth (regenerate). Dashed lines show the position of the autotomy plane. Abbreviations: ac-arm coelom; rnc-radial nerve cord; rvc-radial canal of the water-vascular system. https://doi.org/10.1371/journal.pone.0232981.g001 nerve cord, the radial canal of the water-vascular system, and the arm coelom (compare Fig 1C  and 1D).

De novo transcriptome assembly
We are generating genomic and transcriptomic resources for the brittle star O. brevispina. However, since at the time of writing the assembly and annotation of the full genome for this species is still in preparation, we chose to use a de novo assembled transcriptome as a reference to characterize the Notch pathway target genes.
The de novo transcriptome was generated from 17,318,775 MiSeq and 832,245,006 HiSeq quality filtered and adapter trimmed reads. The single MiSeq library represented pooled samples from intact and regenerating arms at different states of regeneration, whereas six HiSeq libraries corresponded to three control (DMSO-treated) and three DAPT-treated regenerating individuals on day 14 post-autotomy (see Methods). Sequence reads were assembled with Trinity [26,27] into 2,463,269 contigs (1,169,021 Trinity "genes") with the average/median contig length of 421.6/260 nt and contig N50 of 527 nt. The key assembly metrics are listed in Table 1.
To assess the quality of the assembled transcriptome, we performed a series of benchmark tests. First, to assess the representation of reads in the transcriptome, all cleaned Illumina reads were aligned back to the assembled contigs with Bowtie 2 [28]. The vast majority of the reads (94.09%) mapped back to the assembly. Of these, 75.42% of the reads aligned as proper pairs one (21.62%) or more (53.80%) times.
Second, to determine the representation of completely or almost completely (>80% of the length) assembled protein-coding transcripts, we compared our contigs to the reference proteome of the sea urchin S. purpuratus [29], the echinoderm species with the best-annotated genome to date. This analysis showed that 7,397 sea urchin orthologs (out of 35,786) are represented in our transcriptome by full-length and nearly full-length transcripts.
Third, the completeness of the assembly in terms of protein-coding gene content was assessed using BUSCO [30] and the conserved metazoan gene dataset. Out of 978 genes (or 98.7%) in the metazoan database, 966 genes were recovered in the assembled transcriptome as "complete" (i.e., their length fell within two standard deviations of the BUSCO group mean length). Of these complete genes, 282 matched a single contig, whereas multiple copies represented the remaining 684. The high number of "duplicated" genes is a known phenomenon in de novo transcriptome assembly, as even in the absence of any sequencing errors, inherent biological complexity of the transcriptome (e.g., single nucleotide polymorphism and alternative splicing) makes assembly algorithms report multiple isoforms for individual genes [31].

Identification of the genes affected by the Notch pathway perturbation
To identify genes, whose expression changes in response to the Notch signalling perturbation, we performed a transcriptome-wide gene expression comparison between the DAPT-treated cohort and the control (vehicle-treated) animals. These differentially expressed genes comprise potential direct and indirect targets of the pathway. To reduce the redundancy of the de novo assembled transcriptome and to be able to perform transcript expression quantification at the "gene" level, we clustered the 2,463,269 Trinity contigs into 1,012,954 "clusters" using the Corset tool [32]. Corset was applied after the reads representing the libraries from the control and DAPT-treated individuals were aligned back to the assembled transcriptome [33]. Differential gene expression analysis was performed with the DESeq2 package [34], which identified 1,978 significantly up-regulated transcripts and 2,434 significantly down-regulated in the DAPTtreated cohort, as compared to control (DMSO-treated) individuals. For the purposes of this study, deferentially expressed genes were defined as those whose expression in response to DAPT treatment changed by at least the factor of 1.5 in either direction and the associated pvalue of the statistical test adjusted for multiple comparisons was less than 0.05 (Fig 3).
To extract the biological meaning behind these extensive lists of genes, we used the DAVID functional annotation resource [35,36]. We annotated the transcripts representing the Corset clusters by comparing them against the Uniprot database. This approach yielded lists of Uniprot IDs for up-and down-regulated genes, containing 180 and 142 entries, respectively. As the "gene population background" [35], we used a list of 59,110 entries yielded by annotation of all Corset clusters. The percentage of the gene identifiers that were successfully "mapped" to the internal IDs within the DAVID knowledge base were 89%, 88%, and 92% for the list of down-regulated, up-regulated and "background" genes, respectively. In DAVID, we then used the "mapped" genes to perform functional annotation clustering by classifying the input genes into groups based on measuring relationships among annotation terms associated with these genes. Each cluster is assigned an enrichment score, which is defined as the geometric mean of the enrichment p-values associated with each annotation term in the group [35]. The results of the DAVID cluster analysis are shown in Tables 2 and 3. The few genes, whose identifiers failed to map to the DAVID database, were annotated manually (Tables 4 and 5). The down-regulated genes clustered by DAVID include genes that (a) code for extracellular matrix proteins and mediate mediate interactions between cells and the extracellular matrix (e.g., cell adhesion), (b) regulation of multiple signaling pathways (including TGF-β, retinoic acid, JNK, p38, NF-kappa-B, JUN, TNF, Wnt, Toll-like receptor signaling, and EFGR signaling), (c) innate immune response (antiviral response, phagocytosis, lectin pathway), (d) protein ubiquitination involved in protein homeostasis, quality control, immunity, cell death, and regulation of cell signaling transduction, (e) cell survival, proliferation, and differentiation; (f) semaphorin signaling, cytoskeleton organization/remodeling and cell migration, (g) function of the endomembrane system, including the endoplasmic reticulum; (h) cellular homeostasis, The significantly up-regulated transcripts (red) are to the right, and significantly down-regulated transcripts (blue) are to the left. Each gene is represented with a dot. Grey dots represent the genes whose expression level did not change significantly in response to DAPT treatment. The differentially expressed contigs were defined as those whose associated adjusted p-value was less than 0.05 and the log2 fold change in expression exceeded ±0.58. Volcano plot showing the log2 fold change on the X-axis and -log10 of adjusted p-values on the Y-axis. Each gene is represented with a dot. The 1,978 significantly up-regulated transcripts (red) are to the right, and the 2,434 significantly downregulated transcripts (blue) are to the left. Grey dots represent the genes whose expression level did not change significantly in response to DAPT treatment. The differentially expressed contigs were defined as those whose associated adjusted p-value was less than 0.05 and the log2 fold change in expression exceeded ±0.58.). Boxed labels represent UniProt IDs (with label copy number within parenthesis) shown in Tables 2 and 3. including regulation of intracellular ion concentration and ribosomal structure and function ( Table 2).
Out of 14 down-regulated genes that did not get mapped to the DAVID database, seven represented proteins derived from retrotransposon genes (e.g., DNA-polymerase, reverse transcriptase, Pol). The remaining genes encode for three lectins, an antibacterial protein neuromacin and the antiviral tripartite motif-containing protein 5 ( Table 4).
The up-regulated genes in DAVID functional clusters are known to be involved in (a) innate immune activity (antiviral and antibacterial), (b) cell signaling pathways (Toll-like, NFkappa-B, JNK, p38, MAPK, VEGF, TOR, TGF-beta, BMP, Wnt), (c) positive regulation of neuronal and synaptic differentiation; (d) cell proliferation and differentiation; (e) cell death and removal of apoptotic cells by phagocytosis; (f) stress response (e.g., hypoxia, mitochondrial stress, negative regulation of the heat shock response); (g) epigenetic regulation (e.g., maintaining the transcriptionally repressing state by Polycomb proteins), (f) cell adhesion and migration; (g) RNA maturation, processing and splicing, (h) mitochondrial function in cellular respiration and cell death ( Table 3). One of the genes that are up-regulated in response to inhibition of the Notch pathway is E3 ubiquitin-protein ligase Neurl1, an inhibitor of the Notch pathway, suggesting a control mechanism through mutual inhibition. Among 18 up-regulated genes that did not have a match in the DAVID database, seven represented transcripts derived from transposable elements, including a transposase from a Tc3-like DNA transposon, reverse transcriptase genes from three retroelements, and two putative retrotransposon integrase genes. The remaining unmapped genes are involved in innate immune response, cell signaling, cell adhesion, lipoprotein endocytosis, amino acid metabolism. In addition, one of the transcripts matched and ERG-like transcription factor (known protooncogene) ( Table 5).

Discussion
Here we show that the proper function of the Notch signaling pathway is required for arm regeneration in a brittle star. Pharmacological inhibition of the pathway with DAPT antagonist during the first 14 days of regeneration resulted in a significantly reduced overall length of the outgrowth and a smaller number of segments in the new arm. A similar effect of Notch pathway inhibition on regeneration was previously shown for sea urchin spines and podia. Treatment of DAPT resulted in a dosage-dependent reduction of regrowth of those structures post-amputation [20]. This suggests that the involvement of Notch signaling in the regeneration of body appendages is conserved among echinoderm classes. The transcriptome-wide quantitative comparison of gene expression between the control individuals and the individuals with the inhibited Notch pathway provided us with an insight into the biological roles of the pathway in the brittle star arm regeneration. Our analyses show that the Notch signaling pathway directly or indirectly affects more than two thousand downstream genes.
Functional annotation of the genes that significantly changed their expression in response to the pathway perturbation suggested that Notch affects a wide array of biological phenomena in regeneration. Additionally, this work provides evidence suggesting that the activity of some of the mobile DNA elements (transposons) is, at least to some extent, regulated by Notch signaling.
First, multiple components of crucial cell signaling pathways changed their expression in both positive and negative ways in response to Notch pathway inhibition. This data indicates that individual pathways are not working in isolation during brittle star arm regeneration, but rather involved in a complex cross-coordination. This conclusion is consistent with previous studies that established the interaction of Notch signaling with other signaling pathway including, for example, the Bmp and Wnt pathways in the vertebrate heart development and repair [10], interaction between the Notch and Wnt signaling systems in Drosophila development and human cancers [37], and with the Wnt, BMP and Yap/Taz pathways in regulation of neuronal stem cells [38].
Second, a number of studies have previously demonstrated a crucial role of the extracellular matrix (ECM) remodeling in echinoderm regeneration [24,[39][40][41]. However, the upstream molecular mechanisms that initiate and coordinate these changes have heretofore not been described. Our data here suggests that in brittle star arm regeneration, Notch signaling is involved in both positive and negative regulation of different aspects of ECM composition and cell-ECM interaction, including cell migration. Third, perturbation of Notch signalling results in downregulation of innate immune pathways, suggesting that they are positively regulated by the pathway in regenerating tissues. These include lectins, antiviral response, and phagocytosis. This elevated immune response might represent a reaction to the invasion of pathogens and the accumulation of damaged or necrotic cells at the site of autonomy. Alternatively, molecular components of the innate immune response are also known to function as key modulators of regeneration, including wound healing and cell division [42]. For example, innate immune mechanisms are known to facilitate changes in epigenetic states and activate mechanisms required for nuclear reprogramming and cell dedifferentiation [43].
Fourth, the Notch signaling performs its classical functions in regulating cell proliferation, survival, programmed cell death, and differentiation in brittle star arm regeneration. All these cell events are prominently involved in regeneration in O. brevispinum (Mashanov et al., in preparation), and in other echinoderms [44]. Fifth, the Notch pathway controls some of the cellular housekeeping functions, such as regulation of the endoplasmic reticulum, ion balance, ribosomes, RNA processing, and mitochondrial function.
Sixth, Notch affects transcriptional activity (both positively and negatively) of mobile genetic elements (RNA and DNA transposons). Transposons, especially RNA transposons, are known to be active not only in the germ line, but also in somatic cells. Transposons affect the host cell gene expression via various mechanisms, including 1) providing promoter and enhancer sites that change the expression of the host genes, 2) creating splice and polyadenilation sites, and 3) through RNA interference [45]. Although, in some cases, transposon activation is known to cause genetic disorders and autoimmune reactions [45], they do also play important "positive" roles in adult and developing tissues. For example, in mammalian neurogenesis retrotransposition activity of LINE-1 elements contributes to neuronal diversity [46,47]. We have also recently shown that retrotransposons are differentially expressed in the sea cucumber neural regeneration. In response to injury, some of the retroelements increased their levels of expression, but there were others that were suppressed. The cells with activated retroelement activity remained alive and contributed to the newly regenerated tissues [48,49]. We have hypothesized that mobile DNA plays a role in post-traumatic regeneration in echinoderms. The host molecular mechanisms that deferentially regulate the transcriptional activity of mobile genetic elements are largely unknown. This work provides evidence that the activity of mobile DNA is controlled by a cell signaling pathway. Our data also provide further mechanistic insight into how Notch can positively regulated expression of retroelements. One of the genes repressed by active Notch signaling is deoxynucleoside triphosphate triphosphohydrolase SAMHD1 (Table 3), an antiviral protein that is also known to inhibit endogenous retroelements such as LINE-1 [50]. Therefore, in the brittle star arm regeneration, the Notch pathway can activate retrotransposons through the double negative regulation mechanism.
Our data also provide an insight into how the Notch signaling, once activated, might sustain its function in the brittle star arm regeneration. One of the genes that we found to be suppressed by the pathway is a homolog of Neuralized1, which is an antagonist of the Notch pathway [51].
We fully understand advantages and limitations of the whole-transcriptome RNA-seq gene expression analysis approach implemented in this study. On the one hand, this strategy provides us with an insight into an array of direct and indirect tagets of the Notch pathway, including discovery of new roles for genes in regeneration, such transposons. On the other hand, this approach lacks the resolution to distinguish between the immediate targets of the pathways and those that may be several levels farther down the gene regulatory network hierarchy. Nevertheless, the present study yields a valuable resource in the form of a catalog of Notch-affected genes involved in regeneration. Ongoing and future functional genomics work will further unravel the details of the functional relationships between the individual genetic components and their specific functional involvement in post-traumatic tissue formation.

Conclusions
• The activity of the Notch signaling pathway is required for proper post-autotomy regeneration of the brittle star arm. Inhibition of the pathway results in a significant reduction of the outgrowth rate.
• There is evidence of an extensive cross-talk between the Notch pathway and other signaling pathways in the cell.
• Inhibition of the Notch pathway affects a wide range of cellular and molecular processes in the regenerating tissues, including: • composition and remodeling of the extracellular matrix • cell adhesion and migration • innate immune response • cell proliferation, survival, differentiation, and programmed cell death • "housekeeping functions" (functions of the endoplasmic reticulum, ribosomes, mitochondria, RNA processing, ion balance).
• The Notch pathway directly or indirectly regulates (positively or negatively) the activity of several RNA and DNA transposable elements.
• A possible mechanism by which Notch sustains its activity in the regenerating tissues is suppression of Neuralized1.

Animal maintenance and pharmacological treatment
Adult individuals of the brittle star Ophioderma brevispina (Say, 1825) were purchased from the Marine Biological Laboratory (Woods Hole, MA) and were allowed to acclimate overnight in aerated sea water before being subjected to experimental procedures. The DAPT reagent, N-[N-(3,5-Difluorophenacetyl)-L-alanyl]-S-phenylglycine t-butyl ester, was purchased from Sigma-Aldrich (catalog number D5942) and dissolved in DMSO (Dimethyl sulfoxide) to make a 20 mM stock solution.
The animals were divided into two cohorts (5 animals in each)-the DAPT treatment group and the control group. The DAPT treatment cohort was continuously exposed to 3 μM DAPT prepared by diluting the stock solution in filtered seawater. The control group was exposed to the matching concentration of DMSO (vehicle). The animals were kept in glass vials submerged in 200 ml of the DAPT or DMSO solution, respectively. The solutions were changed daily.
The treatments started when the animals still have intact arms; then, after 24 hours, the arms were autotomized as described below and the treatments continued until 14 days postautotomy.
In each individual, all five arms were autotomized by squeezing a single arm segment with a fine forceps. Since the rate of regeneration in serpent stars is known to vary with the position of the autotomy plane along the proximodistal axis of the arm [25], we kept the injury paradigm consistent by autotomizing all arms at the level of the 15th segment (counting from the disk).
On day 14 post-autotomy, three of the five arms of each animal were pooled and used for total RNA extraction. The two remaining arms from each individual were processed for histological analysis.

Histology and image analysis
Tissue samples were fixed overnight in 4% paraformaldehyde prepared in 0.01M PBS (pH 7.4) at 4˚C. Fixed tissues were washed in PBS, decalcified in 10% EDTA, cryoprotected in graded sucrose solutions, and embedded in the Tissue-Tek OCT medium (Sakura). Serial longitudinal sections (10 μm thick) were cut with a Leica CM1860 cryostat, collected on gelatin-covered slides, and incubated at 42˚C overnight. They were then stained with Hematoxylin and Eosin and mounted in the DPX Medium (Electron Microscopy Sciences). The sections were photographed using an Olympus BX60 compound microscope equipped with a SPOT RT camera. The length of the regenerating arm was measured in six animals from each treatment group using the Fiji/ImageJ software [52,53] in calibrated micrographs. The means between the two conditions were compared using Student's T-test in R [54].

Total RNA extraction
The tissue samples for RNA extraction included the most distal segment of the stump and the entire outgrowth (regenerate). The samples from three arms of each individual were pooled together and quickly homogenized in 500 μl of ice-cold TRI reagent (Sigma-Aldrich, 93289) using a disposable Micro Tissue Homogenizer (Kimble, K7496250030). The homogenized samples were briefly vortexed, an additional 500 μl of the TRI reagent were added to each tube, and the subsequent steps of RNA purification were performed as directed by the manufacturer's protocol.

Sequencing
Illumina HiSeq. Three total RNA samples from DAPT-treated individuals (each sample representing a different individual) and three samples from control individuals (DMSOtreated) were used for sequencing on the Illumina HiSeq 2000/2500 platform. The libraries were barcoded, pooled, and sequenced in a single flow cell in the Rapid Mode (2 × 100 bp). These reads were used both for the de novo transcriptome assembly and for gene expression analysis, as reads corresponding to each individual library could be identified and analyzed separately.
Illumina MiSeq. To facilitate the de novo assembly, we also prepared the following single combined sample to maximize the representation of expressed genes in the assembled reference transcriptome. Total RNA was extracted as above from intact animals and from untreated regenerating individuals on days 1, 3, 5, 18, and 30 post-autotomy. Each condition was represented by seven individuals. These 42 RNA samples (6 conditions × 7 individuals) were mixed in equimolar quantities and sequenced on the Illumina MiSeq platform (2 × 250 bp). The final stages of library preparation and sequencing were outsourced to the Duke Center for Genomic and Computational Biology. The raw sequencing reads were deposited at the NCBI Sequence Read Archive (SRA) under the accession number GSE142391 (https://www.ncbi.nlm.nih.gov/ geo/query/acc.cgi?acc=GSE142391).

De novo transcriptome assembly and quality assessment
The code used for the transcriptome assembly and subsequent RNAseq data analysis can be found in the accompanying additional file (S2 File). The raw Illumina reads were trimmed with Trim Galore [55] to remove adapter sequences and low-quality bases at the 3'-end (with the base quality < 20). If a cleaned read was shorter than 20 nt, the entire read pair was discarded. This procedure yielded 17,318,775 MiSeq read pairs (combined library with pooled samples from different arm regeneration stages) with the read length ranging between 20-250 nt. Individual HiSeq libraries representing control and DAPT-treated individuals contained 34,676,875±5,960,462 (mean± standard deviation) read pairs with sequence length ranging between 20-100 nt. The total number of cleaned HiSeq read pairs was 832,245,006. All cleaned Illumina reads, from both the MiSeq and HiSeq technologies, were pooled and assembled with Trinity [26,27] with the minimum reported contig length of 150 nt. The assembled transcriptome is available at the NCBI Gene Expression Omnibus (GEO) database (accession number GSE142391). The quality of the assembly was assessed by running several tests. First, to assess the representation of reads in the de novo transcriptome, we mapped all quality trimmed Illumina reads back to the assembled contigs using Bowtie 2 (Version 2.1.0) [28]. Second, we determined the number of completely or almost completely (>80% of the length) assembled protein-coding transcripts by comparing our contigs to the reference proteome of the sea urchin S. purpuratus [29]. Third, the completeness of protein-coding gene representation in the transcriptome was assessed with BUSCO (v3.0.2) [30] run in the "transcriptome mode" against the evolutionary conserved metazoan gene set (metazoa_odb9, creation date: 2016-02-13, number of species: 65, number of BUSCOs: 978).

Identification of the main components of the Notch signaling pathway
To identify the main components of the Notch signaling pathway in the transcriptome of O. brevispina, we performed BLAST analysis (tblastn, E-value cutoff of 1e − 5; [56]) using Notchrelated genes from the published transcriptomes of H. glaberrima [17] and S. purpuratus [5], keeping all hits. Each transcript identified this way was considered a putative Notch-related gene. We them used the NCBI's conserved domain database (CDD; [57]) to search for conserved domains that are categorize the target genes, including the Notch receptor, ligands (Delta and Serrate), the transcriptional regulator RBPJ, two Notch target genes of the Hes family (orthologous to HES-1 and HES-4), and the Notch signaling modulator Numb.

Differential gene expression analysis
The cleaned reads from each of the six Illumina HiSeq libraries representing the control and DAPT-treated animals were aligned back to the assembled contigs using the Salmon tool(version v0.13.1) [33] with the -dumpEq, -validateMappings, and -hardFilter options. We then used Corset (version 1.08) [32] to mitigate the redundancy issue inherent to de novo transcriptome assemblies. Corset clusters isoforms into "genes" based on the proportion of shared reads and expression patterns and thus outputs gene-level counts.
Quantitative gene expression analysis was performed using the DESeq R package [34]. To optimize memory usage, the raw count matrix was pre-filtered by only keeping rows that contained 10 reads or more. After normalizing the data and fitting the model, the independent hypothesis weighting approach for p-value adjustment implemented in the IHW R package [58] was applied to optimize the statistical power of the analysis. To improve fold change estimates for the gene expression data, a DESeq2 log2 fold change shrinking algorithm was applied. It uses information from all genes to improve the estimates for genes expressed at low levels with high dispersion values but does not change the total number of genes that are identified as deferentially expressed [34]. Genes were considered deferentially expressed if their adjusted p-value was below 0.05 and the fold change in expression level exceeded 1.5 in either direction. The gene expression data were submitted to the GEO database (accession number GSE142391).
Functional annotation of deferentially expressed genes was performed using DAVID (Database for Annotation, Visualization, and Integrated Discovery) [35], which was accessed from within R interface using the RDAVIDWebService package [36]. To prepare data for DAVID, we annotated the Trinity transcripts representing all Corset clusters by matching them to the Uniprot database using BLASTX with the cutoff e-value threshold of 10 −5 . The resulting list of Uniprot accession IDs was used as the background gene set for DAVID. The up-regulated and down-regulated genes were annotated in the similar manner before being imported into DAVID. For functional annotation, we set the following annotation categories: "GOTERM_BP_ALL", "GOTERM_MF_ALL", "GOTERM_CC_ALL", "KEGG_PATHWAY", "BIOCARTA", and "PFAM". DAVID cluster reports were parsed and analyzed using a custom R function (S2 File).