Cellular pathways during spawning induction in the starlet sea anemone Nematostella vectensis

In cnidarians, long-term ecological success relies on sexual reproduction. The sea anemone Nematostella vectensis, which has emerged as an important model organism for developmental studies, can be induced for spawning by temperature elevation and light exposure. To uncover molecular mechanisms and pathways underlying spawning, we characterized the transcriptome of Nematostella females before and during spawning induction. We identified an array of processes involving numerous receptors, circadian clock components, cytoskeleton, and extracellular transcripts that are upregulated upon spawning induction. Concurrently, processes related to the cell cycle, fatty acid metabolism, and other housekeeping functions are downregulated. Real-time qPCR revealed that light exposure has a minor effect on expression levels of most examined transcripts, implying that temperature change is a stronger inducer for spawning in Nematostella. Our findings reveal the potential mechanisms that may enable the mesenteries to serve as a gonad-like tissue for the developing oocytes and expand our understanding of sexual reproduction in cnidarians.


Results
Transcriptomic analysis during spawning induction. To understand the molecular pathways leading to the release of oocytes in Nematostella, females were induced by temperature elevation and light exposure (see "Methods") and then sampled for RNA-Seq analysis at 1, 2, 5, and 8 h during induction in addition to a pre-induction control sampling (Fig. 1A). Transcriptome sequencing produced an average of 16.5 million raw reads per sample, filtered and mapped to the N. vectensis genome (NCBI genome assembly GCA_000209225.1 ASM20922v1). Based on FPKM (fragments per kilobase million) normalized read counts per transcript per sample, we estimated the distance between samples using principal component analysis (PCA) ordination, which demonstrates that although individual variance at each time point exists among the samples, the distance between the baseline and the treated samples grows with time during induction (Fig. 1B). Generalized linear model in DESeq2 revealed ~ 1500 transcripts that were found to be differentially expressed (FDR < 0.05) during spawning induction, as compared to pre-induction levels (Fig. 1C, Supplementary Table S1). MA plots ( Fig. 1D-G) demonstrate that during the induction process, both the number of differentially expressed transcripts and their fold-change magnitudes are gradually increased. This gradual effect is also seen in a heatmap, which additionally shows two distinct groups with opposite trends (Fig. 1C). The first group includes 608 transcripts that were upregulated 8 h after induction began, whereas the second group includes 628 transcripts that were downregulated 8 h into induction.
Gene ontology (GO) enrichment analysis. Next, we analyzed the enriched pathways of the differentially expressed transcripts during the experiment timeline. GO analysis revealed that 88% of the significantly enriched upregulated biological processes and molecular function categories were either specific to the first hour of induction (n = 117, 42%) or shared between two or more time points (n = 126, 46%) (Fig. 2, Supplementary  Table S2). Processes overrepresented among upregulated transcripts included terms related to light reception and chemical stimuli involving numerous receptors, as well as processes linked to the ECM and the nervous system. Among the significantly enriched downregulated biological processes and molecular function categories, 47% (n = 168) were shared between time points and included terms related to cell cycle, oxidation-reduction, and various metabolic processes (Fig. 2, Supplementary Table S2).
Processes related to light, ECM, and cytoskeleton are upregulated during spawning induction. Since light is one of the triggers of spawning induction in Nematostella 26 , we analyzed the GO term sensory perception of light stimulus (GO:0050953), which was upregulated during the first hour of induction (Figs. 2, 3A). The analysis revealed that transcripts associated with this term are expressed mainly at 1 h and 8 h into induction and that many are receptors, specifically G protein-coupled receptors (GPCRs) (Fig. 3A). Moreover, 5 and 8 h into induction, there was an upregulation of transcripts involved in photoreceptor activity (GO:0009881) (Figs. 2, 3B), and 8 h into induction, there was also an upregulation of the photoperiodism process (GO:0009648), though not upholding significance (p = 0.05). Together, these two processes regulated 11 transcripts, five of which are core components of the circadian clock. Among the circadian clock transcripts, Clock and three cryptochromes (Cry1a, Cry1b, and Cry) were downregulated before induction and their expression elevated during the process, peaking after 8 h. The transcript Cry1 displayed an opposite expression trend, as it was upregulated before induction started and decreases during the process (Fig. 3B). Other transcripts that are related to photoreceptor activity, such as melanopsin, rhodopsin, protease usp6 n-terminal-like, rhodopsinlike GPCR and opsin BCRH2-like, displayed a similar expression trend of upregulation during the induction process.
We also analyzed the term ligand-gated ion channel activity (GO:0015276), which was common to all sampling time points (Figs. 2, 3C). The analysis revealed four different GABA A receptor subunits and three transcripts related to the excitatory neurotransmitter glutamate. Two of the glutamate-related transcripts (171792, 132356) displayed similar patterns of upregulation, mainly at the first hour of induction, whereas the third transcript (138860) was strongly upregulated at the eighth hour of induction (Fig. 3C).
We then used STRING to analyze the cytoskeleton and ECM proteins related to the GO terms actin filamentbased movement (GO:0030048), extracellular matrix structural constituent (GO:0005201), and positive regulation of actin filament polymerization (GO:0030838), which were upregulated at the first, second, and eighth hour of the induction process, respectively. The analysis revealed a complex network that can be divided into three main groups: ECM, actin cytoskeleton, and microtubule cytoskeleton (Fig. 3D). Some of the microtubule isoforms were upregulated, whereas others were downregulated and displayed a weak connection to isoforms of motor proteins, such as dynein and kinesin (Fig. 3D). The myosin motor proteins are connected to actin, which showed strong connections to three GTPases: RhoA, Cdc42, and Rac1. While the first two were upregulated, Rac1 was downregulated. Additionally, connected to actin is the integrin receptor Itgb1, which was downregulated and strongly connected to different isoforms of collagens comprising the ECM network, all found to be upregulated (Fig. 3D). Overall, the upregulated processes indicate that spawning induction is perceived by different receptors, some of which are light-sensitive receptors that initiate signaling cascades and cytoskeletal rearrangement. www.nature.com/scientificreports/ Metabolic processes are downregulated during spawning induction. A large group of downregulated processes included seven terms associated with metabolism. Only three of these processes were shared among time points, two of which are related to lipid metabolism, namely fatty acid metabolic process and neutral lipid metabolic process (Fig. 2). KEGG analysis of the term fatty acid metabolic process (GO:0006631) revealed that nine transcripts function in fatty acid degradation, and except one, all of them were downregulated (Supplementary Fig. S1). These results indicate that the degradation of fatty acid and consequently the formation of acetyl-CoA, which is required for ATP production, are reduced during spawning induction.  www.nature.com/scientificreports/ Light effect on the induction process. Next, we performed real-time qPCR on selected transcripts obtained from RNA-seq to test their transcriptional response to light during induction. For that, Nematostella females were induced by temperature elevation under either light or dark conditions and were sampled at the same time points as in the RNA-seq experiment ( Supplementary Fig. S2A). While in both groups spawning rates were about 80%, some transcripts displayed different expression patterns in response to light ( Supplementary  Fig. S2). Out of seven clock-related transcripts examined, five transcripts (Cry1a, Cry1b, Clock, BHLH transcription factor Helt, and glutamate receptor Grin1) were upregulated during light induction and at the eighth hour of the dark induction, whereas one transcript, thyrotroph embryonic factor (Tef), was downregulated as induction proceeded under both light and dark conditions ( Supplementary Fig. S2B-H). We additionally tested two cytoskeleton-related transcripts, tubulin 1 (Tuba1) and Wnt3, and the aquaporin (Aqp) water-selective channel.
With the exception of Aqp, no significant differences in expression levels were found between the two treatment groups ( Supplementary Fig. S2I-K). These results support the findings of light perception during spawning induction but also demonstrate that temperature raise is a stronger signal for spawning.

Discussion
To explore the molecular mechanisms and pathways underlying spawning induction in the basal metazoan Nematostella, we analyzed the transcriptome of females before and during the induction process. Our results revealed a diverse array of up-and-downregulated processes, which lead to successful spawning. Furthermore, a complementary light/dark induction experiment revealed that the expression levels of most tested transcripts changed significantly under both conditions, demonstrating that temperature change is a stronger signal for spawning. These molecular data support the predominance of temperature change over light, which was shown before in spawning experiments 26 .
In cnidarians, light has been shown to have a fundamental role in regulating sexual reproduction. The spawning of the hydromedusa Clytia hemisphaerica and Cladonema pacificum is triggered by light-dark transition and is mediated by opsin photoreceptors from the GPCR family [34][35][36] . In the corals Acropora millepora and Acropora digitifera, the photopigment melanopsin, which belongs to the opsin family, is upregulated during the spawning night, possibly triggering different signaling cascades contributing to gamete release 17,19 . We found upregulation of several GO terms related to light reception and transduction as well as to GPCR and ion channel activity, including melanopsin and rhodopsin receptors from the GPCR family and a glutamate receptor. This suggests that in Nematostella, light-sensitive proteins are responsible for receiving the light signal and initiate different cascades, such as GPCR signaling and glutamate release. Additionally, our qPCR analysis showed that the glutamate receptor Grin1, an ion-gated channel 37,38 , was upregulated at the eighth hour of both light and dark inductions. Previous studies have shown that the glutamate pathway interacts closely with the molecular feedback loops of the circadian clock 39,40 . In A. digitifera and A. millepora, upregulation of transcripts related to glutamate release by photoreceptors was observed during the spawning event 17,19 . Together, these results suggest a conserved role for GPCR and glutamate signal transduction during the reproduction process of Nematostella. However, the exact mechanism is yet to be uncovered.
The core components of the circadian clock, cryptochromes and Clock, have been identified in Nematostella and were shown to have diel rhythmicity. Moreover, there is evidence that Nematostella exhibits circadian behavior and physiology under diel light cycles that are maintained upon the removal of light 20,[41][42][43] . In A. millepora, two cryptochromes that displayed rhythmic gene expression under light:dark cycles were suggested to play a role in synchronizing spawning 18 . Our results showed that the expression of Cry1a, Cry1b, and Clock was upregulated during light induction, while in the dark induction, they were upregulated at the eighth hour only. Since cultivated cultures of Nematostella can be induced in light and dark conditions, our results suggest that Nematostella Crys' and Clock are not essential for initiating the process yet, may play non-circadian roles later during spawning induction. Another clock-related gene, Tef, which belongs to the PAR family of transcription factors, was shown to be expressed in vertebrates with a circadian rhythms 44 and to regulate downstream light-induced genes under circadian control in zebrafish 45 . Our results demonstrated that in Nematostella, Tef is not affected by light, as it was downregulated in both light and dark conditions. This is in line with the findings in A. millepora, where Tef was downregulated prior to and during spawning 17 . It is possible that Tef functions as a repressor of a yet unknown pathway and therefore its downregulation allows the pathway to operate. Although the circadian clock plays an important role in the reproduction of many organisms from different phyla [46][47][48][49] , our results indicate that spawning of Nematostella, during controlled laboratory inductions, is most likely not regulated by the endogenous clock. Notably, the circadian expression pattern of these clock genes is strongly dependent on consistent light:dark cycles 42,50 ; however, the studied sea anemones were maintained in constant temperature  www.nature.com/scientificreports/ and darkness for a long period and were not entrained by daily exogenous cues. Therefore, more investigation is required to understand the role of circadian clock components during spawning induction in Nematostella.
One of the expected changes during spawning induction is the rearrangement of cytoskeletal networks. In the corals Acropora gemmifera, A. millepora, and A. digifera, cytoskeletal reorganization was elevated during and after spawning 17,19,51 . Our results demonstrated that, likewise, extensive changes in the cytoskeleton and its related proteins occur during spawning induction in Nematostella. We have discovered transcripts that are involved in actin regulation, such as the small GTPases RhoA, Rac1, and Cdc42 (Fig. 4). Cdc42 promotes the formation of filopodia 52 and RhoA regulates the formation of linear actin network, focal adhesion and stress fibers, which are the major contractile structures required for cell adhesion during migration 53,54 . When RhoA is activated, it affects the actomyosin contractile network 55 . This stands in line with our results showing upregulation of RhoA and myosins. Rac1 controls the formation of branched actin network and the formation of lamellipodia in the leading edge of motile cells 54 . RhoA and Rac1 usually antagonize each other 56 and the balance between them was shown to control cell shape 57 . The finding of upregulation of RhoA and downregulation of Rac1 suggests that such balance may allow the mesenteries to change shape and create space for oocyte development. Additionally, it was shown that in migrating cells, branched actin network is localized in the leading edge of the cell, while stress fibers are present in the cell body, anchoring it to the surface through focal adhesions 58 . Our findings indicate a similar actin organization, suggesting that cell migration occurs during spawning induction.
The ECM serves as a substrate for cell migration and transfers signals into the cell from the surrounding environment through adhesion receptors such as integrins, which promote reorganization of the cytoskeleton 59,60 . In addition, ECM remodeling plays a crucial part in regulating diverse cellular processes such as cell shape and growth 61 . Our results showed upregulation of many different isoforms of collagen, the most abundant ECM proteins 62 (Fig. 4), possibly indicating ECM rearrangement.
Many of the downregulated GO terms we found are related to metabolic processes, such as fatty acid metabolism. This could imply that during induction, less energy is generated in the mesenteries. In contrast, fatty acid oxidation is increased in the oocytes 63 . We additionally found that transcripts related to cell cycle, an energy consuming process, were downregulated (Fig. 4). Furthermore, DNA, RNA, and protein-related processes were downregulated, suggesting that cell maintenance was reduced. Altogether, these results indicate that there is a decrease in metabolic processes that are not required for spawning. The mesenteries were previously shown to host the synthesis of the fatty acid-rich vitellogenin (Vtg), which is the main yolk protein in Nematostella oocytes 31,63 . In fish, vitellogenesis is accompanied by an increase in fatty acid in the liver, where the process occurs. Decrease in fatty acid in the liver and their increase in gonads were observed before spawning and  www.nature.com/scientificreports/ occurred due to transport of Vtg to the oocyte 64,65 . We, therefore, postulate that the observed decrease in fatty acid breakdown might be a result of Vtg synthesis. Together, these findings may possibly indicate that the mesenteries function as gonads during spawning in Nematostella.
To conclude, this study reveals molecular pathways underlying spawning induction in female sea anemones, N. vectensis. Our results demonstrate that upon the reception of light and temperature signals, processes involving numerous receptors, circadian clock components, the cytoskeleton, and ECM are upregulated, whereas housekeeping processes are mostly downregulated. As shown by the conceptual model in Fig. 4, we suggest that these changes may enable the role of the mesenteries as a gonad-like tissue for the developing oocyte and eventually lead to spawning. Our findings demystify the progression of spawning induction in basal metazoan and could assist in deeper understanding the evolution of spawning patterns of higher organisms that share a common ancestor with Cnidaria.

Methods
Nematostella cultivation. Nematostella were kept in artificial seawater (12.5 ppt, pH 8.0, Red Sea, Israel) at 18 °C in the dark. They were fed five times a week with freshly hatched brine shrimps (Artemia salina) with weekly seawater changes. Spawning was controlled and induced as described 26 . In short, on the day before spawning induction, the anemones were not fed, and spawning was conducted in an incubator under white light by raising the temperature from 18 to 25 °C for 12 h.

RNA-seq experimental design.
For the RNA-seq experiment, 25 well-fed female anemones were collected. As a control, five individual anemones were placed in a dark incubator at 18 °C, sampled before spawning induction, and immediately snap-frozen in liquid nitrogen. Spawning was then induced as described above. At 1, 2, 5, and 8 h during induction, five females per time point were sampled and immediately snap-frozen in liquid nitrogen. Samples were stored at − 80 °C until processing.  The process is triggered by light and temperature elevation. We propose that light stimulates melanopsinlike homologue, which activates rhodopsin and glutamate receptors and initiates GPCR signaling cascades. Other receptors are also activated and together they affect cell cycle, cytoskeleton and ECM organization, focal adhesion, fatty acid oxidation, and circadian clock components. Eventually, these processes lead to successful spawning. Upregulated transcripts are written in purple and downregulated transcripts are in green. www.nature.com/scientificreports/ the RNA-seq experiment, 30 females were induced in the same incubator that was covered to create dark conditions. In addition to a pre-induction sampling that served as a control for both groups, individual sea anemones from both treatment groups were sampled at the same time points as in the RNA-seq experiment at 1, 2, 5, and 8 h during induction and frozen immediately. To ensure a proper spawning process, 87 females were induced by temperature and light and 48 by temperature in the dark; 81% and 77% of them, respectively, released egg sacks.

RNA extraction.
Total RNA was extracted from frozen anemone females using Quick-RNA MiniPrep Kit (Zymo Research, Irvine, USA), according to the manufacturer's instructions with minor modifications. Lysis was performed using 3 mm glass beads in TissueLyser II (QIAGEN) and additional centrifugation at 1000g for 1 min was then performed to precipitate unbroken particles. Genomic DNA residues were removed by Turbo DNase treatment (ThermoFisher Scientific), and final RNA was eluted in 80 μl nuclease-free ultra-pure water.
Purified RNA concentrations were determined using NanoDrop 2000c spectrophotometer (Thermo Scientific). Extracted RNA from the RNA-seq experiment was quality assessed using TapeStation 2200 (Agilent Technologies) and extracted RNA from the qRT-PCR experiment was quality assessed using gel electrophoresis (1% agarose). RNA samples were stored at − 80 °C until use.
Next-generation sequencing and data processing. Total  Differential expression analysis. Differential expression analysis was conducted using generalized linear models (GLMs) in Bioconductor DESeq2 (in R 3.5.1), considering a single factor: induction time 68 . Differential expression was considered to be significant at p-value < 0.05 (adjusted using Benjamini-Hochberg correction), and transcript abundance is presented as FPKM values produced in DESeq2 (Supplementary Table S1). MA plots were generated with DESeq2 package in R language using DESeq2 lfcShrink function. Heatmaps of normalized FPKM values of all the significantly expressed transcripts were generated using the heatmap.2 function from package GPLOTS (v2.17.0) in R. Expected factor-dependent trends were verified using Principal component analysis (PCA) graphs in stats (v3.6.2).
Gene ontology (GO) enrichment analysis. Functional enrichment analysis was conducted using Bioconductor GoSeq (in R 3.5.1) 69 , with the Wallenius bias method, to allow correcting for gene-length biases in differential transcript expression. In addition, functional enrichment was detected using the score-based algorithm GSEA 70,71 considering log 2 fold changes as scores and GO terms were grouped to representative GOs using REVIGO 72 . Non-canonical database for GoSeq was first built using Trinotate 73 . Heatmaps of significantly expressed transcripts within selected terms were generated using the heatmap.2 function from package GPLOTS (v2.17.0). In addition, networks and connections between transcripts were produced using KEGG pathway map (https:// www. kegg. jp) 74 and STRING (https:// string-db. org/).
Real-time qPCR. cDNA was synthesized using High Capacity cDNA Reverse Transcription Kit (Applied Biosystems) following manufacturer's protocol without RNase inhibitor and modification of Step two in the PCR reaction program, which lasted 60 min. Primers were designed based on the exon-exon junction using Primer Express 3.0 software (Applied Biosystems) for a set of 10 transcripts of interest, which were significantly up-or down-regulated in the RNA-seq experiment (sequences were obtained from the published Nematostella JGI database https:// mycoc osm. jgi. doe. gov/ Nemve1). Primer sequences used in this study are listed in Supplementary Table S3. The qPCR reaction was performed using StepOnePluse (Applied Biosystems) 96-well machine with 5 μl Fast SYBR™ Green Master Mix (Applied Biosystems), 0.5 μl of 10 μM primers and 25 ng cDNA template in a 10 μl volume, using reagent samples without cDNA as negative controls as described before 75 . The thermal profile was set to 95 °C for 20 s, followed by 40 amplification cycles of 95 °C for 3 s, 60 °C for 30 s, dissociation cycle of 95 °C for 15 s and 60 °C for 1 min and then brought back to 95 °C. FUE (far upstream element-binding protein 3) transcript (244532) was used as an internal reference control. The constant expression of the FUE level was verified across treatments by qPCR analysis of reactions loaded with equal amounts of cDNA. The quantification analysis was done by the comparative CT method 76 , and all samples were quantified relatively to time 0 (pre-induction). Results are presented as the means with standard error of four biological replicates. While most data were normally distributed according to the Shapiro-Wilk test, the rest were log 10 transformed and retested for normal distribution. The significance of the results was determined using two-sample Student's t test or the nonparametric Mann-Whitney U-test for data that were not normally distributed after transformation (SPSS software v. 20). The results were considered statistically significant if the null hypothesis could be rejected at the p < 0.05 level. www.nature.com/scientificreports/