Some like it hot: population-specific adaptations in venom production to abiotic stressors in a widely distributed cnidarian

Abstract

Background

In cnidarians, antagonistic interactions with predators and prey are mediated by their venom, whose synthesis may be metabolically expensive. The potentially high cost of venom production has been hypothesized to drive population-specific variation in venom expression due to differences in abiotic conditions. However, the effects of environmental factors on venom production have been rarely demonstrated in animals. Here, we explore the impact of specific abiotic stresses on venom production of distinct populations of the sea anemone Nematostella vectensis (Actiniaria, Cnidaria) inhabiting estuaries over a broad geographic range where environmental conditions such as temperatures and salinity vary widely.

Results

We challenged Nematostella polyps with heat, salinity, UV light stressors, and a combination of all three factors to determine how abiotic stressors impact toxin expression for individuals collected across this species’ range. Transcriptomics and proteomics revealed that the highly abundant toxin Nv1 was the most downregulated gene under heat stress conditions in multiple populations. Physiological measurements demonstrated that venom is metabolically costly to produce. Strikingly, under a range of abiotic stressors, individuals from different geographic locations along this latitudinal cline modulate differently their venom production levels.

Conclusions

We demonstrate that abiotic stress results in venom regulation in Nematostella. Together with anecdotal observations from other cnidarian species, our results suggest this might be a universal phenomenon in Cnidaria. The decrease in venom production under stress conditions across species coupled with the evidence for its high metabolic cost in Nematostella suggests downregulation of venom production under certain conditions may be highly advantageous and adaptive. Furthermore, our results point towards local adaptation of this mechanism in Nematostella populations along a latitudinal cline, possibly resulting from distinct genetics and significant environmental differences between their habitats.

Background

Nematostella vectensis is a burrowing sea anemone which specializes in estuarine environments with a unique role as an infaunal predator [1]. These brackish habitats are characterized by variable daily and seasonal abiotic conditions, particularly temperature, salinity, and ultraviolet (UV) light [2,3,4,5,6]. Nematostella has a broad geographic range along the Atlantic and Pacific coasts of the USA and southern Canada where the extent of variation in environmental conditions changes by latitude or location within the estuary [1]. Nematostella tolerates salinities from ~ 8.96 to 51.54‰ and reported temperatures from − 1.5 to 28.5 °C [1, 7] in the field, and laboratory experiments have shown it can acclimate to even broader ranges. Like many coastal invertebrates, Nematostella exhibits extensive genetic diversity with significant population genetic structure throughout its range [8, 9]. Current evidence of genetic structure and a life history that likely reduces dispersal between locations (collective egg masses, demersal larvae) is consistent with limited gene flow between estuaries. The combination of geographically structured genetic variation and differences in environmental conditions is the context where we may expect that populations might be adapted to different ranges of environmental parameters [8, 10, 11].

Previous research with Nematostella adults from different geographic locations has shown evidence consistent with local adaptation for particular phenotypes or genetic loci [12]. Nematostella from locations along the Atlantic coast have temperature-dependent growth rates and thermotolerance consistent with a thermal gradient from low to high latitudes [8]. Similarly, different genotypes vary in their tolerance to oxidative stress, which may be related to genetic variation in the transcription factor NF-κB [13] or superoxide dismutase [14]. Although whole genome comparisons have not yet been completed to identify additional loci where genetic variation is structured between populations, a survey of allelic variation for a subset of gene coding loci suggests it may be pervasive [15]. The diversity of genes potentially involved in adaptation to abiotic variation (such as heat, salinity, and oxidative stresses) could be large. For example, Nematostella has large numbers of heat shock proteins (HSPs) and antioxidant genes [16, 17] likely to be instrumental in mounting a cellular stress response.

Species also may adapt to their biotic environment based on the distribution of their prey and predators. Consistent with most cnidarians, Nematostella relies on its venom system to mediate antagonistic interactions with predators and prey. In particular, sea anemones have a complex decentralized venom system to produce multiple toxins localized to different tissue and cell types [18,19,20,21,22]. Strong evidence supports that in Nematostella, venom production is regulated throughout the life cycle and correlates with the type of interactions it is exposed to. Specifically, at the larval non-feeding stage, it produces toxins to deter fish predators while as an adult it combines both defensive and offensive toxins [18, 23]. The Nv1 toxin is the major component of the adult venom and is produced by ectodermal gland cells at very high levels as it is encoded by multiple gene copies [24, 25]. Nv1 is lethal even at low doses for crustaceans, which could be either predators or prey [18]. Under controlled diel light conditions, the expression of Nv1 is significantly higher during the day when compared with night [26]. In fact, Nv1 is among the most differentially expressed genes in Nematostella. The higher expression during the day correlates with an antipredator response to visual predators like shrimp and fish.

Venoms evolved in many animal lineages. Nevertheless, what they have in common despite their independent origins is chemical complexity and very high production levels. These characteristics led to a hypothesis of high metabolic cost of venom biosynthesis [27]. However, this has only been studied in scorpions and snakes [28,29,30], with no studies in cnidarians. Moreover, even in snakes, there are contrasting views regarding the cost of venom to individuals [31]. The metabolic cost of venom production affects the selection pressure imposed on the venom system and evolutionary mechanisms acting on the venom. Addressing the venom cost in Nematostella is important because it has become a powerful model for studying general principles of venom evolution and ecology.

Regulation of the venom production across Nematostella’s life cycle supports a hypothesis that tuning venom composition based on ecological requirements may be advantageous [18]. Otherwise, anemones would be expected to maintain high venom expression regardless of stage or environment. Stress response mechanisms require additional metabolic cost as they involve production of high levels of specialized proteins, such as chaperones, depending on the severity of the physiological stress imposed by the environment.

In this study, we hypothesized that there might be an advantage to downregulate venom biosynthesis depending on the physiological status of the animal and its exposure to environmental stressors in order to prioritize metabolic resources for supporting the stress response. To test this hypothesis, we combined analyses of physiology, gene expression, and proteomics of Nematostella adults challenged with combinations of environmental stressors to determine the response in the production of venom. Further, we also compared responses between individuals collected from different locations along this species’ native range (i.e., Atlantic coast of North America) to determine to what extent tradeoffs in stress response and venom production differ between these geographically and genetically isolated populations. Our approach allowed us to test these hypotheses related to venom production cost in Nematostella and how local environmental conditions correlate with venom adaptations.

Results

Venom production is metabolically expensive

To study whether venom production is metabolically costly, we depleted Nematostella venom reserves by mechanical stimulation and measured time-dependent changes in respiration rate, which correlates with metabolic rate. A similar method has been used to assess metabolic cost of venom production in snakes and scorpions [28, 29]. Oxygen intake was measured every hour for 5 h following the stimulation (Fig. 1a; Additional file 1: Data S1). After 2 h, oxygen consumption increased by more than 34% and it persisted for the next 2 h with a slight decrease by the fifth hour. Oxygen consumption rates did not increase within 3 h for anemones that did not receive a physical stimulus (Additional file 2: Fig. S1).

Fig. 1

Venom production is metabolically expensive. a Changes in oxygen consumption following fishing line treatment measured by closed-chamber respirometry (Additional file 1: Data S1). The error bars represent standard errors. b Percentage differences of expression for genes encoding toxins (solid lines) and nematocyst structural proteins (dashed lines) following fishing line treatment. c Absolute differences of expression for toxin-coding genes (Nv1, Nematolysin_1b, NveSKT1), nematocyst structural protein NCol3 (dashed line), stress response (Catalase), and housekeeping (HKG4 and HKG5) genes following fishing line treatment. Noticeably, Nv1, NveSkT1b, and Nematolysin 1b are the genes with the largest increase in expression in absolute units. d Gene expression at the 3-h time point represented as a percentage chart. The toxins Nv1, NveSkT1b, and Nematolysin 1b are contributing to a total of 77% of the transcripts among all the genes measured in this experiment. In b, c, and d, gene expression is represented as normalized fluorescence units measured by nCounter technology (Additional file 1: Data S2); the error bars represent standard deviations

To confirm that the elevated metabolic rate involved the biosynthesis of venom, we repeated the venom depletion treatment and measured expression levels of genes encoding toxins produced by both adult males and females, nematocyst structural proteins, genes involved in general stress responses (e.g., cytoplasmic heat shock proteins and superoxide dismutases), and several housekeeping proteins by nCounter technology (Additional file 1: Data S2). In addition to previously characterized venom components, we also included the putative toxins NEP8-like and NveSkT1 that have not been reported earlier (Table 1, Additional file 2: Fig. S2). After 2 h, we observed an increase in the expression of all the genes measured (Fig. 1b). At the 3–5-h time points, expression levels had diverse dynamics: several genes encoding toxins and nematocyst structural proteins (NveSkT1, NEP6, NEP3-like, NEP8, NR2, NCol3) kept increasing while others (Nv1, NEP14) showed wave-like fluctuations. Because the basal expression levels of the Nv1, Nematolysin1b, and NveSkT1 toxins were comparatively high, this twofold increase in transcription contributed to the overall transcriptional activity much more than any other gene measured (they account for 77% of the transcripts among the genes measured at the 3-h time point; Fig. 1c, d). Thus, increased expression of venom components and nematocyst structural genes positively correlates with the increase in the respiration rate.

Table 1 Venom components studied in the present work. Gene models are taken from https://figshare.com/articles/Nematostella_vectensis_transcriptome_and_gene_models_v2_0/807696

Study of toxin production under stress conditions

Nematostella polyps were sampled from five locations distributed from the north to the south along the North American Atlantic coast (Nova Scotia (NS), Maine (ME), New Hampshire (NH), Massachusetts (MA), and North Carolina (NC); Fig. 2a, b) and then cultured in the laboratory under common garden conditions. To determine the regulation of toxin production in response to abiotic stresses, we subjected adults from each location to the heat stress of 28 °C and 36 °C, UV light, and low salinity (5‰) and high salinity (45‰; MA and NC only) water. Expression levels of toxin and stress response genes were measured by the nCounter platform (Additional file 1: Data S3, 4S). To assess the biological significance of the expression changes in this experiment, we set a threshold of 2.4-fold change (corresponding to the maximum change in the expression of the normalizing housekeeping genes (HKG) [18]). In Table S1 (Additional file 3), changes in gene expression that are both biologically and statistically (p < 0.05, Student’s t test) significant are highlighted in bold red.

Fig. 2

Nematostella populations from different climatic conditions respond differently to heat stress. a Distribution of Nematostella populations along the East coast of North America. b Average high July (red) and average low January (blue) temperatures across the populations from the North to the South. The data originate from publicly available sources (https://www.usclimatedata.com/climate, https://en.climate-data.org/north-america/) for the cities closest to the sampling points (Halifax, NS; Portland, ME; Manchester, NH; New Bedford, MA; Wilmington, NC). Solid line—control temperature (20 °C), dashed line—heat stress temperature (36 °C) used in our experiments. c, d Gene expression dynamics under heat stress (+ 36 °C) (c) and low salinity stress (5 ‰) (d) among the populations from the North to the South (Additional file 3: Table S1). If the change in expression of a gene is greater than 2.4 times and p < 0.05 (Student’s t test), the corresponding data point is outlined in bold and labeled by a red asterisk. Only genes with expression change (increase or decrease) higher than 2.4-fold and with statistically significant (p < 0.05, Student’s t test) difference between control and treatment in at least two populations are shown. Error bars represent standard deviations accounting for error propagation of the treatment/control ratios. e, f Expression of toxins in Massachusetts (e) and North Carolina (f) following heat stress. Gene expression is represented as Log₁₀ of normalized fluorescence units measured by nCounter technology (Additional file 3: Table S1). p values calculated by Student’s t test are shown for each gene. If the change in expression of a gene is greater than 2.4 times and p < 0.05, the corresponding p value is shown in red. Error bars represent standard deviations

The lower temperature of 28 °C, high salinity, and UV light exposure did not result in a substantial change in toxin expression (Additional file 3: Table S1; Additional file 2: Fig. S3). However, under the high heat stress of 36 °C, we observed changes in the expression of toxins and HSPs (Fig. 2c). In the ME, NH, and MA populations, Nv1 production decreased beyond the threshold and up to 25-fold in the NH population compared to the control conditions (20 °C); however, in the ME population, the difference in Nv1 expression was not statistically significant (p > 0.05, Student’s t test). Additionally, NEP14, NEP16, NEP4, NEP8-like, and Nv1 toxin genes exhibited reduced expression under heat stress in the MA population (Fig. 2e, f). In the NS and NC populations, the expression level of Nv1 and NEP toxins did not change more than the threshold of 2.4-fold. HSP expression showed interesting dynamics (Fig. 2c): the expression change was less pronounced in the most northern NS population (most of the changes were not statistically significant; only HSP90A expression significantly increased compared to the control), then increased in ME, NH, and MA and decreased in NC. Thus, HSP expression showed an opposite trend to Nv1 expression: the highest increase in HSP expression in the NH and MA populations corresponds to the largest decrease in Nv1 expression.

Physiological stress induced by low salinity resulted in decreased Nv1 expression in NS, NH, and MA populations with the highest change in NS sample (10-fold) (Fig. 2d); however, other toxin genes were not affected (apart of small change in NEP16 in MA). The expression of HSP70A and HSP90A increased in NS, NH, and MA.

To dissect regulation of venom production further, we focused on the MA and NC populations, which showed significantly different responses to temperature and had similar Nv1 expression levels under control conditions (Additional file 3: Table S1A). According to temperature data from publicly available sources, MA and NC also have quite different temperature profiles, including, for example, the difference of 3.7 °C between high average July temperatures (Fig. 2b). The different thermal regimes between MA and NC habitats are supported by our temperature monitoring data acquired by loggers which had been placed into ponds in March 2016 and measured temperature until December 2016 (Fig. 3a; Additional file 1: Data S5). In the NC location, temperature reached 36 °C or higher on 110 days while in MA only on 38 days during the monitored period (Fig. 3b).

Fig. 3

Thermal regimes differ between Nematostella habitats in Massachusetts and North Carolina. a Water temperature in native Nematostella habitats recorded in March–December 2016. Measurements were taken every 20 min by a temperature logger (Additional file 1: Data S6). b Number of days when the water temperature reached 36 °C or higher in March–December 2016

Transcriptome annotation and Gene Ontology analysis of the stress response

Nematostella in its natural environment, like most organisms, experiences stressors in combinations. We compared these former responses to temperature, salinity, or UV light with a combination treatment of all the three abiotic conditions. Adults from NC and MA were subjected to high temperature (36 °C, 24 h), UV light (6 h), and high salinity (40‰, 24 h). Control conditions were 20 °C, 15‰, and darkness. Animals were snap frozen (4 per sample), RNA was extracted, and Illumina RNA-seq was performed in parallel for stress and control conditions. A principal component analysis (PCA) revealed that > 70% of the variation observed across all the samples could be attributed to location and stress exposure (Additional file 2: Fig. S4).

Our differential expression analysis identified 3637 transcripts that were differently expressed between controls and experimental groups for MA and NC populations (Fig. 4b). Under the combined stress conditions, 738 genes were differentially expressed in both populations, 542 genes in NC, and 592 genes in MA (Fig. 4a, Additional file 3: Tables S2–S4). Within 1280 transcripts that were differentially expressed in Nematostella from NC, 457 contained Gene Ontology (GO) information comprised of 4667 different GO groups across these transcripts. For the MA population, there were 1330 transcripts that were differentially expressed, 513 of which contained some GO information comprised of 5405 different GO groups across these transcripts.

Fig. 4

Differential gene expression among Massachusetts and North Carolina populations after “combined” stress. a Venn diagram of differentially expressed genes that are shared and unique for the MA and NC populations (Additional file 3: Tables S2 – S4). b Heatmap of differentially expressed genes across treatments and populations. c Subclusters of differentially expressed genes: SC1 (blue)—genes upregulated in NC, SC2 (green)—genes downregulated following stress, SC3 (purple)—genes upregulated in MA, SC4 (pink)—genes upregulated following stress. d Heatmap of toxin genes and stress response genes. The latter are the most differentially expressed genes between the control and stressed conditions in both populations

We used comparative GO analysis to identify functionally important GO groups in relation to environmental stress response. However, our analysis did not identify any characterized functionally important groups more highly represented than other GO groups (Additional file 3: Table S5). The GSEA analysis was able to more accurately identify GO groups that were variable across treatments. In both populations, the “response to heat” GO term (GO:0009408) was enriched (Additional file 3: Table S6). The most differentially expressed genes between the control and stressed conditions for both populations can be traced back to stress response genes such as HSP20 and HSP70 (Fig. 4c). Beyond these stress response genes, our GSEA also revealed an enrichment in the upregulation of transcription factors (TFs) in both populations following stress response. Differences in upregulation of specific transcription factors were observed among the populations, with 30 TFs upregulated in both populations, 27 TFs upregulated only in MA, and six TFs upregulated only in NC (Additional file 3: Table S7). Clustering of differentially expressed genes generated four distinct subclusters (SC1–SC4; Fig. 4c) sharing similar expression dynamics. These four subclusters can be defined as genes upregulated following stress (SC4), genes downregulated following stress (SC2), genes upregulated in NC (SC1), and genes upregulated in MA (SC3; Additional file 3: Table S8). GSEA performed on subclusters confirmed “response to heat” GO term (GO:0009408) to be enriched in the subcluster consisting of genes upregulated following stress (SC4).

In the control groups, many of the differentially expressed transcripts had relatively low levels; therefore, even a small change in gene expression resulted in a dramatic fold change. When considering both LogCPM (log counts per million) and TPM (transcripts per million), the change in Nv1 becomes more apparent. The highly expressed Nv1 toxin gene (TPM = 1018.6 ± 339.5; TMM [weighted trimmed mean of M values] = 847 ± 243 in MA; TPM = 4470.7 ± 1397.7; TMM = 3824 ± 1112 in NC) for the control is the most downregulated locus in the MA population (TPM = 133.1 ± 65.5; TMM = 150.6 ± 69.2) (Additional file 3: Table S3), but not in the NC population (TPM = 2126.4 ± 303.5; TMM = 2109.8 ± 156.7) (Additional file 3: Table S4).

Northern (MA) but not southern (NC) Nv1 expression reacts to combined stress conditions

Among toxin genes, only Nv1 in MA showed a substantial change under stress conditions, where it decreased 5.5 times, p = 0.01 (Fig. 4c). The expression of NveSkT1, Nematolysin1b, nematocyst toxins (NEP3, NEP4, NEP3-like, NEP8, NEP14, NEP16, NEP6), and structural proteins (NCol3 and NR2) remained more stable in both populations as changes were in the range of 0.4–2.1-fold and in most cases were not statistically significant. Among the stress response genes that we analyzed under individual stress conditions, heat shock protein genes (HSP70B, HSP70A, HSP90A) exhibited an increased expression level of 11- to 821-fold. In contrast, the expression of the CuZnSod3 superoxide dismutase gene remained stable while MnSod1 decreased approximately twofold (Fig. 4d).

To study the effects of the stress conditions on Nv1 production at the protein level, we used LC-MS/MS (Fig. 5, Suppl Data 4). It was noticeable that control label-free quantification (LFQ) values decreased between our measurements apparently due to degradation or aggregation of Nv1 peptides. Thus, to control for this technical variation, the samples were run in pairs of control and stress for each population in triplicates. In each pair from MA, the LFQ value for the stress sample was approximately twice lower than that in the control. On the other hand, in NC, stress LFQ values were slightly higher than those in the control. Thus, the observed trend supports our findings on the transcriptomics level.

Fig. 5

Nv1 production at the protein level under control and “combined” stress conditions in Massachusetts (a) and North Carolina (b) populations. The measurements were done by LC-MS/MS (Additional file 1: Data S5); the numbers correspond to technical replicates

Discussion

In the current study, we provide a suite of molecular and physiological evidence that toxins are metabolically costly to produce for a venomous species and that individuals from different populations along a thermal gradient exhibit unique modulations in toxin expression when experiencing a range of abiotic stressors.

One of the significant findings for this research is the first substantial evidence for a cost of venom production in a cnidarian. We have identified that an increase in expression of venom components and nematocyst structural genes correlates with an increase in respiration rate, a measure used to assess the general metabolic rate. This increased metabolism might be an evidence for elevated biosynthesis rates, including the biosynthesis of venom components. Similar approaches have been applied to study the metabolic cost of venom production in other animals [28, 29]. Toxin proteins have been hypothesized to be energetically expensive to produce, and thus, venomous species should be under strong selection to modulate the quantity and composition of venom cocktails [27, 36]. Studies in reptiles and arachnids [27] support this hypothesis where individuals modulate feeding preference or behavior depending on their energetic state or the relative venom quantities (where starved individuals produce less venom). This hypothesized tradeoff may be supported by stony corals, which strongly downregulate their small cysteine-rich peptides (SCRiPs)—later identified as neurotoxins [37]—under heat and UV stress [38], and the sea anemone Anthopleura elegantissima, which downregulates the Nv1 homologs Anthopleurin-C (called “Toxin PCR 7” in that study) and APETx1 under similar stresses [39]. Thus, adjusting venom production to the metabolic status may be an evolutionary adaptation in cnidarians in general.

Our data from Nematostella indeed show that among venom components, the highly abundant toxin Nv1 is the most downregulated under heat stress conditions in multiple populations. These massive shifts in Nv1 expression may be driven in part by its high copy number resulting in the production of this essential toxin to be highly dynamic. Similar patterns may be observed in other venomous lineages that also exhibit highly expanded gene families encoding toxins, such as metalloproteinases in rattlesnakes [40]. The investment of energy into heat response (e.g., production of numerous heat shock proteins, HSPs) appears to be traded-off with other high-cost physiological processes that do not contribute to survival under heat stress (e.g., venom production). It is hypothesized in both sea anemones [27] and snakes [40] that selection is acting on the expansion of gene families encoding toxins to drive increased protein production. Additionally, this mechanism may also allow for drastic and rapid shifts in toxin expression to meet biotic and abiotic factors.

In sea anemones, venom release involves secretion of venom components from ectodermal gland cells [24] and discharge of single-use nematocysts through mechanical and/or chemosensory cues [41]. In our experiment, after the mechanical stimulation, venom components need to be replenished and also new nematocytes have to maturate and produce nematocysts. Therefore, multiple processes involved in venom regeneration may explain the response where venom-related genes as well as other genes were upregulated after 2 h following the treatment.

Another significant finding in the current study is the result that Nematostella polyps from different geographic locations differ in their production of toxins when experiencing the same environmental stressors. While we measured significant changes in expression of multiple toxins under different exposures to high temperature, low and high salinity, and UV light, the largest population-specific impact was at the higher temperature exposure. Responses to the heat shock followed an inverted bell-shaped profile along a latitudinal temperature gradient: a statistically significant decrease in Nv1 expression was observed in the NH and MA populations but not in the NS and NC populations (Fig. 2). One of the possible explanations would be adaptation to the local temperature conditions for these locations. In NC, where the average high July temperature is 32.1 °C and average low January temperature does not go below 0 °C, Nematostella would be exposed to heat shock of 36 °C regularly. On the other hand, in NS, the average high temperature is 22.6 °C meaning that Nematostella is very unlikely to experience a heat shock of 36 °C in this habitat and would be rather adapted to a lower temperature range. The lack of a gene expression response in NS animals may be due to a lack of adaptations to such an extreme temperature compared with typical high temperatures causing severe stress and might be explained by a general metabolic shutdown. In ME, NH, and MA, the average high temperatures are between 26.2 and 28.4 °C and low temperatures are below 0 °C indicating that Nematostella would be exposed to both high and low temperatures and thus is adapted to more diverse temperature regime. Our temperature monitoring data in Nematostella habitats confirm that NC population may be exposed to 36 °C for 6 months each year while the MA population only during 4 months. We hypothesize that heat response mechanisms in all populations are inducible based on the expression of HSP genes, but there is only a statistically significant tradeoff for physiological processes in NH and MA populations. The NC anemones have potentially evolved other physiological mechanisms to cope with higher temperatures and to avoid the tradeoff with the production of toxin proteins. However, it is possible that for this population a tradeoff might be occurring at even higher temperatures. A similar mechanism of plasticity in environmental stress response was described in the coral Porites astreoides [42].

Interestingly, the combined stress conditions (heat, UV light, and high salinity) provoked a similar response as heat stress alone: Nv1 was downregulated in MA but not in NC. This result indicates that single stressors (temperature) and combined stressors induce a similar transcriptional response in the most abundant venom component, Nv1. Analysis of the changes in functional Gene Ontology groups under stress between MA and NC populations revealed that, in concordance with the previous lines of evidence, heat shock proteins underwent significant upregulation in both populations. Differences in the regulation of Nv1 biosynthesis may be explained by a number of transcription factors differentially upregulated between the populations under stress. Factors underlying adaptation of the NC population to retain high venom production levels under stressful conditions are yet to be discovered. Overall population structure needs to be examined in order to determine if divergence in toxin gene expression levels correspond to genetic divergence throughout their distribution. Future work may identify unique Nv1 coding regions which correspond to certain geographic locations or isolated localities.

Despite the fact that multiple genes were differentially expressed between the populations, our analysis of enriched GO terms did not reveal any specific functional pathways upregulated in the NC population and potentially keeping venom production beneficial or neutral under the stress. One of the reasons might be incomplete annotation of cnidarian genomes and GO terms due to high divergence from bilaterian organisms with well-annotated GO terms such as vertebrates, Drosophila melanogaster and Caenorhabditis elegans.

Our earlier and present work have demonstrated that Nv1 toxin is regulated by diverse environmental and endogenous factors: abiotic stress, light/dark cycles, and life stages. Unlike other venom components, it is encoded by multiple gene copies potentially underlying the high expression level and providing higher regulatory flexibility. Because responses to the abiotic stress differ between populations experiencing a gradient of environmental conditions, our data suggest that regulation of venom expression is adjusted and potentially “optimized” to conditions reflective of their local environments. Thus, modulation of venom production appears to be one of the evolutionary adaptations in Nematostella. Adaptation of this trait is plausible as we showed that production of venom and venom-delivery cells is metabolically expensive and the balance between the benefits and costs of venom production might change dramatically between habitats due to different abiotic conditions.

Conclusions

Venom production is metabolically costly in the cnidarian Nematostella. The abiotic stress requiring additional energy investments results in venom downregulation in Nematostella. While we tested these phenomena under several stress conditions, the largest population-specific impact was observed at higher temperature exposure. Under heat stress, the highly abundant toxin Nv1 is the most downregulated gene in multiple populations. It appears that the investment of energy into heat response is traded-off with other high-cost physiological processes that do not contribute to survival under heat stress such as venom production. Interestingly, Nematostella polyps from different geographic locations differ in their production of toxins when experiencing the same environmental stressors possibly due to distinct genetics and significant environmental differences between their habitats. Taking together, our results with anecdotal observations from other cnidarians suggest that this might be a universal phenomenon in Cnidaria with important implications for adaptation.

Methods

Venom discharge and oxygen consumption

To provoke venom discharge, Nematostella polyps from lab population were subjected to mechanical stimulation with a gelatin-coated fishing line imitating interaction with predators and prey following an approach described by [43]. Briefly, 25% gelatin-coated probes were dried for 24 h, rehydrated in 15 ppt saltwater, and pulled through the tentacles of adult Nematostella polyps in one single motion. Closed-chamber respirometry is an effective means to measure oxygen consumption to determine shifts in metabolic rate for animals in different conditions. We used comparative respirometry to measure changes in oxygen consumption in anemones following mechanical stimulation with a fishing line. Anemones were placed in 3-ml water-jacketed respiration glass chamber (1 anemone/chamber), given approximately 30 min to open their tentacle crown, and then the tentacles were probed with the monofilament line. The respiration chambers were then closed and measured continuously for 5 h. Oxygen uptake by individual anemones was measured using Clarke-type oxygen electrodes (YSI, USA). Two-point calibration of electrodes was performed before each day, and continuous data acquisition of oxygen concentrations was made using a BIOPAC Data acquisition system (BIOPAC, USA). Time effects on respiration were calculated relative to the baseline respiration for each individual during the time course and compared with anemones that received no physical stimulus (control). The experiment was repeated with 4 biological replicates for the physical stimulus and 2 replicates for the control.

To confirm that respiration changes correlated with changes in expression of venom genes, venom discharge treatment was repeated and animals were sampled and frozen in liquid nitrogen every hour for 5 h for RNA extraction. Every time point was sampled in triplicates, 3 animals/replicate.

Stress conditions

Animals originating from wild populations (NS, ME, NH, MA, and NC) were maintained under controlled laboratory conditions (room temperature, 15‰ artificial sea water (ASW)) for at least several months after collection from the field. Before treatments, the animals were acclimatized at 20 °C for 24 h in the dark in 15‰ ASW. For the individual stress conditions, animals were exposed to 28 °C, 36 °C, 5‰ ASW, UV light, 45‰ ASW (only MA and NC), or control conditions (20 °C, 15‰, dark) for 24 h. Animals were placed into tubes and frozen to obtain 3 replicates for each condition, 2 animals/replicate, total 87 samples. For the combined stress treatment, half of the acclimatized animals from MA and NC populations were transferred into 40‰ ASW and exposed to 36 °C for 24 h. For the first 6 h of treatment, the UV lamp was on. Control animals were kept at 20 °C for 24 h in dark in 15‰ ASW. After the treatment, animals were placed into plastic tubes and frozen in liquid nitrogen to make three replicates for each condition and 4 animals/replicate, total 12 samples.

Quantification of gene expression

Total RNA was extracted from the frozen samples using RNeasy Mini Kit (Qiagen, Germany). RNA quality was assayed with a Bioanalyzer NanoChip (Agilent Technologies, USA).

For the samples collected following venom discharge and individual stress samples, the nCounter platform (NanoString Technologies, USA; performed by Agentek Ltd., Israel, and MOgene, USA) was used. One hundred base pair probes (Additional file 3: Table S9) were designed to specifically bind to transcripts encoding toxins, nematocyst structural proteins, and stress response genes. For normalization of expression, a geometric mean of expression levels of five genes with stable expression across development was used (similarly to [18, 23]). Fold change and absolute change relative to the 0-h time point were calculated for the samples collected at 1-h, 2-h, 3-h, 4-h, and 5-h time points after venom discharge. Additionally, at the 3-h time point, the expression level relative to the total expression of all the genes measured in this experiment (100%) was calculated and represented as a pie chart (Fig. 1d). For the 87 individual stress samples, the measurement was performed in two nCounter batches. For each population (for each nCounter batch separately), every stress condition was compared to the control and p values (Student’s t test) and fold changes were calculated.

For the 12 “combined stress” samples, cDNA libraries were constructed by the Sense kit (Lexogen, Austria), pooled together, and sequenced by the NextSeq500 platform (Illumina, USA) with 400 million read depth and 40-bp paired end reads (performed at the Center for Genomic Technologies, The Hebrew University of Jerusalem). Raw reads were submitted to the NCBI Sequence Read Archive database (BioProject ID: PRJNA601530). The raw reads were mapped to Nematostella gene models [44] by Star aligner [45]. Across the 25,721 predicted gene models, the number of unique mapped reads varied per samples (24,161,373–34,332,469), resulting in an average of 1166 reads mapped per predicted gene. Differential expression analysis was performed by EdgeR Bioconductor package [46]. This allowed us to identify transcripts with significant deviations in expression level. In our EdgeR analysis, we identified transcripts that were differentially expressed across all treatments in combination, as well as within each focal locality. Metrics of fold change (TPM, TMM, and LogCPM) were used as a threshold to screen potentially informative transcripts for further evaluation. A PCA analysis was conducted using the built-in Trinity toolkit script PtR to ensure no outliers [47].

Proteomics

For measurements of toxin production at the proteomics level, only females from MA and NC populations were used. After the “combined stress” treatment, the animals were frozen (3 animals/tube) and lysed in 8 M urea and 400 mM ammonium bicarbonate. Females were used to control for sex-specific variation observed in Nv1 expression [12]. The lysates were centrifuged (22,000×g, 20 min, 4 °C), protein concentration was measured with BCA Protein Assay Kit (Thermo Fisher Scientific), and 10 μg aliquots of protein were sent for LC-MS/MS analysis by a Q Exactive Plus mass spectrometer (Thermo Fisher Scientific) at the Proteomics Center of the Alexander Silberman Life Sciences Institute, The Hebrew University of Jerusalem. All the procedures were performed as described by [18]. The samples were analyzed in technical triplicates. The mass spectrometry data have been deposited to the ProteomeXchange Consortium [48] via the PRIDE [49] partner repository with the dataset identifier PXD016943 and are also included into Additional file 1 (Data S5).

Transcriptome annotation and GO analysis

Using information derived from the previously predicted and annotated gene models [44], we evaluated whether differentially expressed transcripts may have carried both functional information based on TMM normalized fold change values as well as LogCPM and TPM values. TMM normalization aides when scaling for variation in library size as well as transcript diversity being sampled [50]. The use of TMM normalization fold change alone aided in identifying transcripts expressed at very low levels in the control groups with potentially high functional importance. LogCPM is also a measure of log counts per million contrasting control and treatment, providing more weight to differences in gene expression when they are more highly expressed across all treatments. In contrast, TPM is an expression value that considers transcript length variation when measuring the number of reads per one million being mapped back to a transcript [51]. This measurement is more stable across samples, permitting comparisons between samples and treatments.

Beyond individual transcripts, we evaluated changes in functional Gene Ontology groups (GO) between the different populations across control and stressed environments. This was done by two different approaches: first, we used custom Python scripts as previously described in [22], which divided TMM values across associated GO terms. The GO groups were then grouped based on semantic similarity and identified using REVIGO [52]. The second approach we performed was gene-set enrichment analysis (GSEA) on the differentially expressed genes using GOseq [53]. This required genes to first be annotated against the Swiss-Prot database (accessed 27 October 2019) using BLASTp and gene ontologies mapped. Differentially expressed genes were partitioned into clusters by cutting the hierarchical tree at a height of 60%. Clusters were then manually curated to generate four distinct subclusters and GOseq performed on each to identify enriched and depleted GO terms.