Transcriptomic responses to grazing reveal the metabolic pathway leading to the biosynthesis of domoic acid and highlight different defense strategies in diatoms

Background A major cause of phytoplankton mortality is predation by zooplankton. Strategies to avoid grazers have probably played a major role in the evolution of phytoplankton and impacted bloom dynamics and trophic energy transport. Certain species of the genus Pseudo-nitzschia produce the neurotoxin, domoic acid (DA), as a response to the presence of copepod grazers, suggesting that DA is a defense compound. The biosynthesis of DA comprises fusion of two precursors, a C10 isoprenoid geranyl pyrophosphate and l-glutamate. Geranyl pyrophosphate (GPP) may derive from the mevalonate isoprenoid (MEV) pathway in the cytosol or from the methyl-erythritol phosphate (MEP) pathway in the plastid. l-glutamate is suggested to derive from the citric acid cycle. Fragilariopsis, a phylogenetically related but nontoxic genus of diatoms, does not appear to possess a similar defense mechanism. We acquired information on genes involved in biosynthesis, precursor pathways and regulatory functions for DA production in the toxigenic Pseudo-nitzschia seriata, as well as genes involved in responses to grazers to resolve common responses for defense strategies in diatoms. Results Several genes are expressed in cells of Pseudo-nitzschia when these are exposed to predator cues. No genes are expressed in Fragilariopsis when treated similarly, indicating that the two taxa have evolved different strategies to avoid predation. Genes involved in signal transduction indicate that Pseudo-nitzschia cells receive signals from copepods that transduce cascading molecular precursors leading to the formation of DA. Five out of seven genes in the MEP pathway for synthesis of GPP are upregulated, but none in the conventional MEV pathway. Five genes with known or suggested functions in later steps of DA formation are upregulated. We conclude that no gene regulation supports that l-glutamate derives from the citric acid cycle, and we suggest the proline metabolism to be a downstream precursor. Conclusions Pseudo-nitzschia cells, but not Fragilariopsis, receive and respond to copepod cues. The cellular route for the C10 isoprenoid product for biosynthesis of DA arises from the MEP metabolic pathway and we suggest proline metabolism to be a downstream precursor for l-glutamate. We suggest 13 genes with unknown function to be involved in diatom responses to grazers. Electronic supplementary material The online version of this article (10.1186/s12867-019-0124-0) contains supplementary material, which is available to authorized users.


Background
Phytoplankton survival depends, among other factors, on the capability to defend the cells against grazing zooplankton. Predation is the main source of phytoplankton mortality [1] and diatoms have evolved a variety of strategies to reduce predation. Thick silicate frustules and spiny armors, formation of toxins, adaptable cell sizes and formation of long chains are examples of defense traits in diatoms [2]. Defense mechanisms induced by the presence of predators require cells to sense and recognize the predator from a distance. Diatoms exhibit complex signaling mechanisms that allow perception of environmental cues such as the presence of gametes of opposite sex and the presence of bacteria [3,4]. Diatoms of the genera Pseudo-nitzschia and Skeletonema respond to predator cues from copepods: at least two Pseudo-nitzschia species by inducing toxin production, and Skeletonema marinoi by shortening the chain length [5][6][7]. The induced defense responses are both diatom specific and predator specific [5][6][7]. Thus two phylogenetically related diatoms, Fragilariopsis cylindrus and Niztschia frigida, do not induce formation of domoic acid (DA). Pseudo-nitzschia seriata elicits induced DA production in response to predator cues from herbivorous copepods, but when exposed to a copepod carnivore, a response is not elicited [8]. Predator cues as referred to in this paper are info-chemical signals that mediate interactions between two individuals and result in an adaptive response in the receiver, here specifically copepodamides excreted by copepods [9]. Presently, 21 copepodamide compounds have been identified and their composition is copepod species specific [10]. Pseudo-nitzschia seriata increases DA production as a response to the naturally occurring mixture of copepods and to at least three isolated copepodamides in ecosystem-relevant doses [10,11]. Cellular signal transduction depends on the perception of external cues, but how the diatoms or other phytoplankton species perceive the predator cues is still poorly understood. Perception by signal transduction via G protein-coupled receptors is proposed for the dinoflagellate Alexandrium catenella exposed to a heterotrophic dinoflagellate grazer [12]. When exposed to different grazers (copepods), the induced defense response and the discrimination between copepod species were suggested to be regulated by protein kinases and calcium signaling [13]. To our knowledge, only one other study has investigated the molecular mechanism behind grazer-induced defense mechanisms in a diatom. Gene regulation and metabolomic response in Skeletonema marinoi as a response to grazers showed that transcripts linked to G protein-coupled receptors and to nitric oxide synthesis were differentially expressed, and that Skeletonema reacts the same way in the presence of two different copepod species [14].
Several species of the diatom genera Pseudo-nitzschia and Nitzschia pose a threat to organisms in the marine environment due to their ability to produce DA, a marine biotoxin known to cause severe harm up the food web, and amnesic shellfish poisoning in humans. To date 52 Pseudo-nitzschia species are described and exactly half of these (26) plus two Nitzschia species have been reported to produce DA [15]. Pseudo-nitzschia was demonstrated to be a source of DA already in 1988, but the ecological role of toxin production is still not fully understood. Pseudo-nitzschia cells react to several abiotic and biotic environmental factors such as changes in light, temperature, salinity, pH, pCO 2 and nutrient levels, as well as to the presence of bacteria and grazers by changing the production of toxins [16,17].
As a result of a major research effort which began during the early 1990s, the biosynthetic pathway of DA is gradually being uncovered. Using 13 C-and 14 C-labelled precursors to detect the major metabolic pathways for DA synthesis, it was suggested that DA originates from condensation of a C10 isoprenoid with a tricarboxylic acid (a product of the citric acid cycle) [18,19]. The precursor of the C10 isoprenoid was suggested to be geranyl pyrophosphate (GPP) and originate from acetyl CoA through the mevalonate (MEV) pathway [20] and the tricarboxylic acid proposedly derives from l-glutamate [18]. GPP is synthesized through the isopentenyl pyrophosphate pathway and originate either through (1) the proposed MEV pathway located in the cytosol and/or (2) an alternative route in the plastid, the methyl-erythritol phosphate metabolic (MEP) pathway [21] (Additional file 1: Figure S1). Both pathways are present in diatoms but the MEP pathway is not complete in all diatoms [21,22]. The early findings have been confirmed and it was demonstrated that the amino function in DA is generated by nucleophilic displacement of GPP by an unknown intermediate [23]. Recently, six compounds with potential functions in the condensation of GPP and l-glutamate were identified and their structure described [24]. The same study also confirmed N-geranyl-l-glutamic acid as a potential intermediate. Four candidate genes, DabA-D, that can code for the condensation of GPP and Keywords: Pseudo-nitzschia, Fragilariopsis, Grazer induced defense, Domoic acid, Gene expression, Methyl-erythritol phosphate metabolic pathway, Geranyl pyrophosphate, l-Glutamate, Proline L-glutamate towards DA were recently identified [25]. These four genes were upregulated under phosphate limitation or elevated pCO 2 , both conditions known to induce DA production [26]. DabA catalyzes the N-geranylation of L-glutamate to form N-geranyl-l-glutamic acid, while DabC and D catalyze dainic and isodomoic acids. The final isomerization reaction to complete the condensation to DA remains unknown. Another tricarboxylic acid, proline, a structural analog to DA, is suggested to be implicated in DA biosynthesis. Based on isotopic labeling pattern of DA, a coupling of DA production to the proline metabolism was suggested [18]. Proline could be an upstream precursor of DA, and a correlation between Pseudo-nitzschia cells with high accumulation of DA and low proline content has been found [27].
The principal aim of this study is (1) to acquire information on genes involved in the metabolic processes of DA synthesis in a species of Pseudo-nitzschia. For revealing whether the MEP or the MEV pathway is a part of the metabolic pathway, and whether we can confirm previous suggestions on genes involved in DA synthesis. (2) To acquire information on the molecular mechanism of perception of predator cues by studying transcriptional processes in the two diatoms Pseudo-nitzschia and Fragilariopsis. The target species tested are P. seriata, known to induce DA when exposed to predator cues, and a species of Fragilariopsis not known to induce this defense mechanism.

Domoic acid content and growth rates
The cellular toxin content in Pseudo-nitzschia increased by sevenfold (t-test P = 0.032) as a response to predator cues exposure, while no increased toxin production was found in the controls (t-test P = 0.19). Fragilariopsis did not produce DA neither in response to grazer cues nor in the control. Both species grew exponentially and Fragilariopsis had a higher division rate than Pseudo-nitzschia (t-test P < 0.05) 0.64 ± 0.01 day −1 and 0.12 ± 0.06 day −1 , respectively. (See details on DA levels and division rates in Table 1).

Overview of gene expression profiles in Pseudo-nitzschia
and Fragilariopsis when exposed to predator cues Of the two diatoms, only Pseudo-nitzschia showed differential gene expression when exposed to predator cues (Fig. 1). In total 1128 genes were differentially expressed, of which 814 genes were upregulated and 314 genes downregulated compared to the control (Additional file 2: Tables S1 and S2). Based on the KEGG and KOG databases, 285 of the regulated genes (25% of the total) have a functional annotation. The remaining 75% of the regulated genes have a BLAST hit but could not be assigned to any specific function. These contigs are comparable to sequences in the data bases but their function is unknown. The assembly information is given in Table 2.
Several genes involved in major metabolic pathways such as amino acid, carbohydrate and lipid metabolism are differentially expressed, with the majority (> 70%) of genes being upregulated compared to the control (Fig. 2). Differentially expressed genes assigned to the category of cellular information processing mainly showed a higher expression compared to the control. Genes were differentially regulated, both up and down, in the category of genetic information processing with functions in folding, sorting and degradation of RNA as well as replication and repair. In the category of environmental information processing, the 22 differentially expressed genes were equally either down-or upregulated (Fig. 2). Most importantly five upregulated genes have a function in the isopentenyl pyrophosphate pathway, assigned to the category of terpenoid and polyketide synthesis (highlighted green in Figs. 2 and 3), these are candidate genes for the C10 isoprenoid in the biosynthesis of DA. The majority of contigs could not be assigned to a function. All of the contigs with highest upregulation belong to this category.

Signal transduction
We searched the transcripts for genes involved in cellular signal transduction, as the perception of external cues must transduce cascading molecular responses. A vesicular neurotransmitter transporter with a known function as extra cellular membrane interaction transporter (PSN0026807) is downregulated. Translation initiation factor 4E (PSN0006232) and serum/glucocorticoid-regulated kinase 2 (PSN0000383) are upregulated as well as GTP-binding protein (PSN0015960) and a RING-type E3 ubiquitin transferase (PSN00338470). Downstream regulator for G-binding proteins, mitogen-activated protein kinase (PSN0000026) and a mitogen-activated protein kinase homolog (PSN001617) are downregulated. The results indicate that Pseudo-nitzschia transduce signals from the copepods via G protein pathways. Nitric oxide signaling may be involved in the response to grazers and a nitric oxide synthase-interacting protein (PSN0002048) is upregulated. Genes in the category of signal transduction differentially expressed in Pseudo-nitzschia as response grazers are listed in Table 3.

Candidate genes for stress
There is strong indication that Pseudo-nitzschia cells are under stress when exposed to grazers. Various genes known to be a response to stressors are upregulated. Heat shock proteins and transcription factors are among the highly upregulated genes, like the HSP90 and one of its homologs HSP90-7, HSP40, heat shock factor protein 2, heat stress transcription factor A-2, and cytochrome P450 genes. Genes indicating cellular stress in Pseudonitzschia as response to the presence of grazers are listed in Table 4.

Common grazing responses in diatoms
In Fragilariopsis, not a single gene is differentially expressed as a response to grazers (Fig. 1). When comparing the regulated nucleotide sequences in Pseudonitzschia with a similar study conducted on Skeletonema [14] we found congruence in upregulated genes. We found 13 genes, (out of only 34 and 55 matches, depending on time point in Skeletonema) that match upregulated genes in Skeletonema (Additional file 1: Figure S2; Additional file 3: Tables S5 and S6). Unfortunately, all these 13 genes lack annotations, which makes it difficult to reconstruct mechanisms for the response (Table 5). However, 10 of these 13 genes seem to be functionally related i.e. they show high similarities at the nucleotide level  (Additional file 1: Figure S3; Additional file 3: Table S9). The oxide synthase-interacting protein (PSN0002048) is upregulated, but this not a homolog to c9307_gl_il in [14] or MMMETSP104-20121108 1594 in [22].

Candidate genes for domoic acid biosynthesis
Based on evidence and indications in earlier studies of DA production in Pseudo-nitzschia, we searched specifically for genes involved in the regulation of the isopentenyl pyrophosphate pathway, such as acetyl-CoA, isopentenyl pyrophosphate (via both the MEV and MEP pathways), GPP biosynthesis, l-glutamate biosynthesis including proline degradation into l-glutamate, and the citric acid cycle. All genes involved in both pathways to isopentenyl pyrophosphate were confirmed to be present in the transcriptome of Pseudo-nitzschia and Fragilariopsis. Five out of seven genes coding for enzymes involved in the MEP pathway were upregulated under grazer-induced DA induction. The gene encoding for isopentenyl pyrophosphate isomerase responsible for the condensation into GPP was also upregulated, as was the DabA gene involved in the condensation of N-geranyll-glutamic acid ( Table 6, Fig. 3). None of the six genes in the proposed MEV pathway were significantly differentially expressed (Fig. 3). The non-significant regulation of the genes in the MEV pathway is randomly up or downregulated. Genes directly involved in the synthesis of l  Table 6), indicating proline to be the downstream precursor for l-glutamate.

Comparison with other studies on the biosynthetic pathway in Pseudo-nitzschia by gene expression analysis
Gene expression in Pseudo-nitzschia multiseries at low and high DA production induced by silicate limitation has been investigated [28]. When comparing exponential vs. stationary phase, the study suggests 12 genes to be involved in DA production. Ten of these are present in Pseudo-nitzschia seriata but were not differentially expressed at induced DA conditions. By investigating three Pseudo-nitzschia species, one of which is toxigenic, Pseudo-nitzschia multistriata, one gene involved in nitrite oxide synthesis was confirmed to be present only in the DA-producing strain [22]. Neither this gene sequence nor any other hit for a nitric oxide synthesis enzyme was present in P. seriata in our study. Gene expression in P. multistriata was investigated under two DA-inducing conditions, phosphate limitation and elevated pCO 2 [25]. We looked for regulation of genes involved in either of the isoprenoid pathways within the   Table 6 P. multistriata data set available on JGI, but did not find direct evidence for either pathway. One gene upregulated under all three conditions (grazer exposure, phosphate limitation and elevated pCO 2 ) is PSN000450, which encodes for an enzyme in the MEP pathway ( Table 6, Fig. 3). All four (DabA-D) genes from the [25] study are upregulated in P. seriata (Table 6, Fig. 3) and are not present in Fragilariopsis.
Among the 1028 regulated nuclear contigs in P. seriata, 441 match the 788 protein sequences in P. multistriata. Out of those, 23 are regulated in P. seriata exposed to grazers, and in P. multistriata under   Figure S2; Additional file 3: Tables S3 and S4).

Diatom response to grazers
We report here > 1000 genes differentially regulated in Pseudo-nitzschia as a response to copepod grazers while no differential gene expression was detected in Fragilariopsis under the same conditions. The two phylogenetically closely related phytoplankton species have  therefore evolved different response strategies to threats from grazers. The lack of response in gene expression levels in Fragilariopsis indicates that in this taxon, protection against grazing might have evolved as a constitutive resistance mechanism, a trait that is always present. Hence, the results indicate that Fragilariopsis has not developed a sensory system to detect the predator cues from C. finmarchicus in the way shown here for Pseudonitzschia and for Skeletonema [14]. Fragilariopsis species are known to be fast growing and able to adapt to various types of habitat, and they successfully prevail in polar regions, both in the pelagic and in sea ice, as well as being globally distributed [29,30]. A thick silica frustule as in Fragilariopsis, opposed to a thinner frustule in Pseudo-nitzschia, can provide mechanical protection-an armor that the copepod either cannot break or, if consumed, may pass undamaged through the gut [31]. The lack of response at the gene expression level does not indicate grazer defense to be an inducible trait, but does not exclude that other Fragilariopsis species may modify frustule thickness on demand. Inducible defenses in diatoms triggered by presence of copepods are still scarcely studied. To date this trait has been documented only in Skeletonema and Pseudo-nitzschia. We compared our results with the available data demonstrating gene expression in Skeletonema as a response to the same grazer species, C. finmarchicus, at two time points [14,31]. There are only 34 (time point 1) and 55 (time point 2) matches with the contigs from our Pseudo-nitzschia datasets, thus highlighting the genetic differences in the defense traits of the two diatoms (chain length versus toxin production). However, 13 of the genes are upregulated in Pseudo-nitzschia and Skeletonema at both time points. This strongly suggests those genes to be candidates for common grazer responses in diatoms. The response of Pseudo-nitzschia to the predator cues, i.e. the induced DA production, supports previous findings [6-8, 10, 32]. The induction of toxins as a response to predator cues implies that DA is an induced defense mechanism. However, the defense role is still speculative as DA-containing cells are ingested by their predators independent of their toxicity. Phycotoxins may, however, be effective after ingestion. Thus, mortality has been detected in copepods after feeding on toxic Pseudonitzschia [8], and reduced escape response has been seen in copepods having fed on toxic Pseudo-nitzschia [33]. No cost of DA production is detected in this study in terms of reduction in the growth rate of Pseudo-nitzschia and the same has been seen in [6,7,32]. Reduction of growth rate correlated, however, with high DA production in [8]. One of the prerequisites for inducible defenses is that they save the organism energy, but reports of costs in relation to induced defense are few and this aspect needs further attention. High costs of plasticity may over time be eradicated by evolution [34,35] and low costs (e.g. minor reduction in growth rate) can go undetected in laboratory experiments where factors such as nutrients and energy (light) are not limited [36]. Small differences may not be statistically significant but may have a strong ecological impact when scaled up to natural populations [37]. Consequently, growth as a physiological and metabolic prime parameter to estimate cost and the physiological status of the cells, might not be a sufficient way to detect the metabolic investment of unicellular organisms. Regulation of genes in the category cell growth and death were upregulated, indicating that the Pseudo-nitzschia cells are forced to relocate energy, possibly because of the allocation costs of the DA production. Cellular regulative processes and changes in metabolic pathways mirror metabolic costs in terms of energy equivalents such as consumption of ATP and NADPH. These are for example needed for the synthesis of new enzymes and metabolites. Various genes of the KEGG and KOG categories might be involved in costs of production of amino acids and carbohydrate metabolism that are essential for all protein biosynthesis, such as the provision of glucose for cellular respiration.

The role of MEP in the biosynthesis of domoic acid
Given that DA is synthesized by GPP and l-glutamate, the results of the present study reveal that the MEP pathway is a very likely precursor pathway for synthesis from GPP and further to DA. The majority of the genes in this pathway are upregulated in Pseudo-nitzschia when cells are triggered to produce DA compared to the control. The presence of both the MEP and the MEV pathways has been identified in various toxigenic and non-toxic Pseudo-nitzschia species, and in one Fragilariopsis species [22]. All genes involved in both pathways were confirmed to be present in the cDNA library of P. seriata and Fragilariopsis. However, we did not detect any differential gene expression of genes coding for enzymes in the MEV pathway as a response to predator cues. The MEP pathway was relatively recently fully described [38,39] and has the same central role as the MEV pathway, to form isoprenoids. Isoprenoids are a class of organic compounds essential for all plants, for forming sterols, brassinosteroids, cytokinins, phytols, plant hormones and carotenoids [21]. The two pathways are located in different compartments of the cell, the MEV pathway in the cytosol and the MEP pathway in the plastid. Our results thus suggest that the formation of DA occurs in the plastid.
We did not detect any regulation of genes involved in the processes of the citric cycle for the proposed precursor of l-glutamate. The only regulation detected that can give hints to the downstream precursors are two genes in the proline metabolism. This supports the hypothesis that proline is coupled with the synthesis of DA. In theory, proline can be metabolized to l-glutamate as illustrated on Fig. 3. However, proline is a stress-related amino acid in many plants and algae [40][41][42] and may be regulated for other purposes than DA synthesis.

Conclusions
We propose that the C10 isoprenoid product for biosynthesis of DA arises from the plastid MEP pathway rather than the cytosolic MEV pathway, as five out of seven genes involved in the MEP pathway are upregulated in Pseudo-nitzschia cells triggered to produce DA. We did not detect any gene-regulation that can be traced to the proposed l-glutamate from the citric acid cycle as precursor for DA, but suggest that the l-glutamate precursor could derive from the proline metabolism. Further, we demonstrate that the two phylogenetically closely related species have evolutionary distinctly different strategies for coping with grazer threats and only Pseudo-nitzschia responds with an induced defense when exposed to predator cues. Finally, we suggest 13 genes with unknown function to be involved in the responses of diatoms to grazers.

Organisms and experimental set up
The diatom strains were established from samples collected in Disko Bay, Greenland and permission for sampling was issued by The Government of Greenland, Naalakkersuisut. Pseudo-nitzschia seriata (strain Disko 8) was isolated, June 2013; Fragilariopsis sp. (strain A4-14) was isolated, May 2014. Biovolume of Disko 8 cells is 1695 µm 3 and of A4-14 is ~ 830 µm 3 . The strain of P. seriata Disko 8 has previously been shown to produce DA in the exponential growth phase as a response to predator cues from Calanus finmarchicus as well as to isolated compounds of predator cues, a group of polar lipids named copepodamides [9]. Both strains were cultured at the University of Copenhagen in L-medium [43] at 4 °C and a light:dark cycle of 16:8 m −2 s −1 . Approximately 3 weeks prior to the experiments, the cells were cultured in L-medium with additional NH 4 to avoid differences in concentrations among treatments due to NH 4 excreted by the copepods [44]. Calanus finmarchicus, originally collected from the Trondheim fjord and kept in culture at Biotrix Trondheim was used as source of predator cues. Diatom cells at a concentration of ~ 4000 cells mL −1 were placed into two-chambered incubators; the chambers were 720 mL each and connected via two apertures of 7.5 cm in diameter. The apertures were covered with a 5-μm mesh, which separated the organisms but allowed the predator cues to pass through [7]. The incubators were placed on a plankton wheel which rotated at 2 rounds m −1 for 24 h. They were sampled for cellular toxin content, diatom cell concentration and gene expression analysis using a 100-mL syringe. After a 24-h acclimation period to the experimental conditions, six copepods were placed in one of the chambers to produce predator cues during the experiment. The control was without exposure to animals. Calanus finmarchicus produces copepodamides both when grazing and when starving [8]. Prior to the experiment, the copepods starved for 24 h so that the gut content was cleared. To ensure a constant supply of predator cues during the experiment, the animals were feeding on Pseudo-nitzschia or Fragilariopsis in the adjacent chamber, respectively.
The experiment was conducted in a temperature controlled room at 4 °C with light intensity of ~ 100 µmol photons m −2 s −1 and a light:dark cycle of 16:8. The experiment was initially carried out in quadruplicates, but the analyses was conducted in triplicates (for experimental workflow see Additional file 1: Figure S4). After 3 days, the experiment was terminated.

Diatom concentration and toxin measurements
After sampling for RNA extractions (below), 2 mL were sampled for diatom cell concentration and fixed in acidic Lugol's solution. The cells were counted in a Sedgewick-Rafter chamber using an inverted light microscope (Olympus CKX31 at a 100× magnification), with a minimum of 400 cells.
The growth rate (µ) was determined using an exponential model: C is the cell concentrations at t 1 start and t 2 end of the experiment.
Division rate per day was calculated as: µ is the growth rate from Eq. 1. For DA measurements 10 mL samples were transferred to glass tubes and centrifuged for 5 min at 4000g, supernatants were discarded and pellets transferred to 2 mL centrifugation tubes (Eppendorf, Hamburg, Germany) and frozen at − 20 °C. DA was measured using liquid chromatography coupled with tandem mass spectrometry following [45]. In brief, samples were measured on a Sciex API 4000QTrap hybrid triple quadrupole-linear ion trap mass spectrometer (Sciex, Darmstadt, Germany) coupled to a LC 1100 liquid chromatograph (Agilent, Waldbronn, Germany). Separation was performed on a reverse-phase chromatography on a C8 phase (50 × 2 mm, 3 μm Hypersil BDS 120 Å (Phenomenex, Aschaffenburg, Germany) at 20 °C. The flow rate was 0.2 mL min −1 and gradient elution was performed with two eluants, wherein eluant A was water and B was acetonitrile/water (95:5 v/v), and both contained 2.0 mM ammonium formate and 50 mM formic acid. DA was detected in the positive mode by using the quantitative mass transition m/z 312 > 266 and the qualitative transition m/z 312 > 161. An external four point calibration curve (10, 50, 100 and 500 pg µL −1 ) was used for calibration and the detection limit (S/N = 3) was determined as 0.7 pg µL −1 .

Harvesting and RNA extraction
After 3 days of incubation with predator cues, the incubators were taken off the wheel. The incubators were carefully turned around minimum 10 times to ensure equal mixing of the culture prior to sampling. From each container, 60 mL were poured into sterile 4 × 15 mL sterile centrifugation tubes (Sigma-Aldrich) and centrifuged at 4 °C and 3300g for 10 min (Algera X-12R, Beckam Coulter, USA). Resulting cell pellets were pooled and immediately mixed with 1 mL 60 °C hot TriReagent (Sigma-Aldrich, Steinheim, Germany) and transferred to a 2 mL sterile centrifugation tube containing acid washed glass beads. The cells were lysed using a tissue lyser (Power lyser 24, Qiagen, Hilden, Germany) at maximum speed for 45 s and immediately afterwards frozen in liquid nitrogen and stored at − 80 °C until further use.

RNA-isolation
The cell lysate was thawed on ice and 200 µL of pure chloroform was added to each vial and vortexed for 15 s. The mixture was incubated for 5 min at room temperature, shaken and incubated again for 5 min, and afterwards centrifuged for 15 min at 4 °C with 12,000g. The upper aqueous phase was transferred to a new vial, 2 µL of 5 M linear acrylamide (10-20 µg mL −1 ) and 1/10 volume of 3 M sodium-acetate (pH 5,2-5,5) was added, the mixture was shaken and 1:1 of 100% isopropanol was added. The mixture was then vortexed and incubated overnight (up to 14 h) at − 20 °C to precipitate RNA. The RNA-pellet was collected by 20 min centrifugation at 4 °C and 12,000g. The liquid was carefully removed and the pellet was washed twice, first with 1 mL 70% aqueous EtOH, centrifuged for 10 min at 4 °C and 12,000g and afterwards with 96% EtOH and centrifuged for 5 min. After removing as much of the EtOH as possible, the pellet was left to air dry for 5-10 min. Finally, the pellet was dissolved in 20 µL RNA/DNA RNase free water (Qiagen, Hilden, Germany). The amount of the RNA was checked on a Qubit 2.0 ® fluorometer (Life Technologies).
A RNA purity check was performed using a NanoDrop ND-100 spectrometer (PeqLab, Erlangen, Germany). The integrity of the RNA was examined using the Nano Chip Assay with the 2100 Bioanalyzer device (Agilent Technologies, Böblingen, Germany). High quality RNAs, OD 260/280 > 2 and OD260/230 > 1.8, with intact ribosomal peaks (obtained from the Bioanalyzer readings) were used for building cDNA libraries.

cDNA libraries
Using the reagents provided in the TruSeq Stranded Total RNA LibraryPrep Kit High Throughput, Illumina ® and following the provided manual, the polyA containing mRNA molecules were purified and fragmented. The RNA was reverse transcribed into double stranded cDNA library by fragmenting and using primers with random hexamers into first strand cDNA by reverse transcriptase. The RNA template was then removed and a replacement strand synthesized, incorporating dUTP in place of dTTP to generate double stranded (ds) cDNA. To prepare the ds cDNA for hybridization onto a flow cell, a single adenylate nucleotide was added to the 3′ ends of the blunt fragments to prevent them from ligating to one another during the adapter ligation reaction. A corresponding single thymine nucleotide on the 3′ end of the adapter provides a complementary overhang for ligating the adapter to the fragment. For ligating multiple indexing adapters to the ends of the ds cDNA RNA Adapter plated provided with the Illumina ® TruSeq ® kit. After ligation, the ds cDNA was amplified with PCR. The cDNA library quality was checked on an Agilent DNA 1000 Chip Assay with the 2100 Bioanalyzer device (Agilent Technologies, Böblingen, Germany).
Each replicate was indexed prior to sequencing, so all replicates were demultiplexed and mapped separately against the reference assembly High-throughput RNA sequencing was performed on a NextSeq Illumina at the AWI.

RNAseq analysis
Using the CLC Genomics Workbench 10 (Qiagen), a de novo assembly was produced from the RNASeq reads obtained for each species. One replicate of each species failed, hence the RNAseq analysis was carried out in triplicates. The annotation of the assembled contigs was done by BLASTx search against the Clusters of Orthologous Groups of proteins, database (KOG) (downloaded 2015.02 from ftp://ftp.ncbi.nih.gov/pub/) [46] and the SwissProt database (release 2016.08) with an e-value cut-off of < 10e −05 . The Kyoto Encyclopedia of Genes and Genomes (KEGG) orthology identifiers (http://www.kegg.jp/) were mapped from SwissProt identifiers at http://www.unipr ot.org. The alignment