Insights into the genetic basis of predator‐induced response in Daphnia galeata

Abstract Phenotypic plastic responses allow organisms to rapidly adjust when facing environmental challenges—these responses comprise morphological, behavioral but also life‐history changes. Alteration of life‐history traits when exposed to predation risk have been reported often in the ecological and genomic model organism Daphnia. However, the molecular basis of this response is not well understood, especially in the context of fish predation. Here, we characterized the transcriptional profiles of two Daphnia galeata clonal lines with opposed life histories when exposed to fish kairomones. First, we conducted a differential gene expression, identifying a total of 125 candidate transcripts involved in the predator‐induced response, uncovering substantial intraspecific variation. Second, we applied a gene coexpression network analysis to find clusters of tightly linked transcripts revealing the functional relations of transcripts underlying the predator‐induced response. Our results showed that transcripts involved in remodeling of the cuticle, growth, and digestion correlated with the response to environmental change in D. galeata. Furthermore, we used an orthology‐based approach to gain functional information for transcripts lacking gene ontology (GO) information, as well as insights into the evolutionary conservation of transcripts. We could show that our candidate transcripts have orthologs in other Daphnia species but almost none in other arthropods. The unique combination of methods allowed us to identify candidate transcripts, their putative functions, and evolutionary history associated with predator‐induced responses in Daphnia. Our study opens up to the question as to whether the same molecular signature is associated with fish kairomones‐mediated life‐history changes in other Daphnia species.


| INTRODUC TI ON
Organisms are challenged throughout their lives by environmental changes that have an impact on the health and fitness of each individual. A given phenotype that is advantageous in one environmental setup might become disadvantageous in another. In general, organisms have two possibilities to cope with environmental changes: return to the ecological niche by behavioral (i.e., migration) or physiological changes, or change the boundaries of their ecological niche by genetic adaptation (Van Straalen, 2003). The former is achieved at the phenotypic level and described as a phenotypic plastic response, while the latter is a genetic adaptation process, where genotypes with a higher fitness pass on their alleles to the next generation.
In some instances, kairomones can be detected by their prey, inducing highly variable as well as predator-specific responses that reduce their vulnerability. These predator-induced responses are a textbook example of phenotypic plasticity (Tollrian & Harvell, 1999) and have been reported in detail for a variety of Daphnia species (e.g., Weider & Pijanowska, 1993;Boeing et al., 2006;Yin et al., 2011;Herzog et al., 2016).
Daphnia are small branchiopod crustaceans and are a model organism widely used in ecology, evolution, and ecotoxicology (e.g., Lampert, 2011;Miner et al., 2012;Picado et al., 2007). Members of this genus link trophic levels from primary producers to consumers in freshwater ecosystems and are, therefore, vulnerable to high predation risk (Lampert, 2011). Extensive changes in behavior, morphology, and life-history traits have been observed in response to predation and predation risk. The responses induced by invertebrate predators include morphological changes such as the formation of helmets in D. cucullata (Agrawal et al., 1999) and D. longispina (Brett, 1992) and the formation of neck teeth in D. pulex (Tollrian, 1995). Vertebrate predator cues have been shown to induce behavioral changes linked to diel vertical migration (Cousyn et al., 2001;Effertz & von Elert, 2017;Hahn et al., 2019) as well as changes in life-history traits in D. magna. The specificity of such predator-induced responses by vertebrate and invertebrate kairomones has been shown, for example, for the D. longispina species complex from the Swiss lake Greifensee (Wolinska et al., 2007). The documented changes in life-history traits included a decrease in size at maturity when exposed to fish kairomones and an increase when exposed to kairomones of the phantom midge larvae, a predatory invertebrate of the genus Chaoborus.
Although phenotypic plastic responses to predation risk have been extensively studied in the ecological and genomic model organism Daphnia, their genetic basis is not well understood (Mitchell et al., 2017;Weiss, 2019). Linking predator-induced responses to the underlying genome-wide expression patterns has been attempted from different perspectives (length of exposure time, species, and experimental conditions) in Daphnia. Orsini et al. (2018) investigated the effect of short-term exposure to fish kairomones (several hours) in D. magna, revealing no change in gene expression. Yet another study identified over 200 differentially expressed genes in response to invertebrate predation risk in D. pulex, of which the most prominent classes of upregulated genes included cuticle genes, zinc-metalloproteinases, and vitellogenin genes (Rozenberg et al., 2015). Finally, a study on D. ambigua under vertebrate predation risk revealed ~50 responsive genes involved in reproduction, digestion, and exoskeleton structure (Hales et al., 2017).
Our goal is to investigate the genetic basis of life-history shifts in response to vertebrate predation risk. Daphnia galeata is the ideal candidate species to address this question, since this species does not show diel vertical migration behavior (Stich & Lampert, 1981) or severe morphological changes in the presence of vertebrate predator cues even after long exposure, but strong phenotypic variation of life-history traits under vertebrate predation risk (Tams et al., 2018).
With a combined approach, we aim to understand the complexity of responses to environmental changes such as those induced by predators, which are known to vary across Daphnia species. We applied a transcriptomic approach (RNA-sequencing), followed by differential gene expression, gene coexpression, and orthology analysis. Gene coexpression network analysis allows to infer gene functions because of the modular structure of coexpressed genes and their functional relations; often coexpressed genes share conserved biological functions (Bergmann et al., 2003;Subramanian et al., 2005). A further benefit of the coexpression network analysis lies in the opportunity to correlate gene expression and external information (Langfelder & Horvath, 2008) thus simplifying the process of candidate genes identification. Additionally, orthology analysis allows revealing functional roles as well as the evolutionary history of transcripts. The degree of conservation of the predator-induced response can be estimated by finding orthologous genes in species having diverged million years ago (Cornetti et al., 2019

| Experimental organisms
This study was conducted on two D. galeata genotypes originally hatched from resting eggs collected from Müggelsee (northeast Germany). This lake is inhabited permanently by a diverse fish population (Fischereiamt Berlin, 2013). A previous study involving 24 genotypes (clonal lines) from four different lakes revealed that the variation for some life-history traits increased for genotypes exposed to fish kairomones within the Müggelsee population (Tams et al., 2018), meaning a broader range of phenotypes were displayed for that life-history trait. We chose the genotypes M6 and M9 which differed in all of their life-history traits and were at the contrasting ends of the phenotypic range exhibited by D. galeata exposed to fish kairomones. Genotype M6 displayed a phenotype which matured later, produced less offspring, and stayed smaller under predation risk. Genotype M9 displayed the opposite phenotype, that is, matured earlier, produced more offspring, and became larger under predation risk (Figure 1). Hence, the specific clonal lines showed opposing responses for each life-history trait.

| Media preparation
ADaM (Klüttgen et al., 1994) was used as the basic medium for fish and Daphnia cultures. Two types of media, fish kairomone and control, were prepared and used for breeding and experimental conditions as detailed in Tams et al. (2018). Briefly, fish kairomone medium was obtained by maintaining five ide (Leuciscus idus) in a 20L tank for 24 hours prior to medium use. All media were filtered (Whatman, membrane filters, ME28, Mixed cellulose-ester, 1.2µm) prior to use and supplemented with 1.0 mg C L-1, P rich Acutodesmus obliquus.
Media were exchanged daily (1:2) to ensure a nutrient-rich environment and a constant fish kairomone concentration. The algae concentration was calculated from the photometric measurement of the absorbance rate at 800 nm. Cetyl alcohol was used to break the surface tension during breeding and the experiment to reduce juvenile mortality (Desmarais, 1997). Breeding and experimental phases were conducted at a temperature of 20°C and a 16-h light/ 8-h dark cycle in a brood chamber with 30% of its maximum light intensity (Rumed, Type 3201D).

| Experimental design and procedures
Each genotype was bred in kairomone-free medium (control) and in fish kairomone medium (predation risk) for two subsequent generations before the start of the experiment to minimize interindividual variances. To this end, 20 egg-bearing females per genotype were randomly selected from mass cultures. From these females of unknown age, neonates (<24 hr) were collected and raised under experimental conditions in 750 ml beakers at densities of < 40 neonates per beaker. They served as grandmothers (F0) for the experimental animals (F2). Based upon previous work (Tams et al., 2018), we started the second (F1) generation after 16-20 days to ensure that offspring from the 3rd to 5th brood were used to start the next generation. The third generation of experimental individuals (F2) was started after 18 days. At the start of the experiment, a pair of neonates was introduced in the experimental vessels (50 mL glass tube) to compensate for juvenile mortality. Before the release of the first brood, on day 6, one of the individuals was randomly discarded whenever necessary so that only one individual remained in each vessel. During the 14 days of the experiment, neonates were removed every 24 hr and the number of broods of each experimental female was documented before media renewal. The adult females (F2) were pooled (n = 20) and homogenized in RNAmagic (Bio-Budget technologies, Krefeld, Germany). Five biological replicates were produced per experimental condition (environment) and per genotype, resulting in a total of 400 individuals (two genotypes x two environments x 20 individuals x five biological replicates). Two of these replicates were backup in case of a downstream failure, F I G U R E 1 Reaction norms of selected life-history traits of experimental genotypes (mean ± SE) from a previous experiment (Tams et al. 2018) emphasizing the opposing environmental effect on life-history traits for the two genotypes. (a) Age at first reproduction in days. (b) Total number of broods per female (c) total number of offspring per female (d) somatic growth rate in µm per day. "Blue": genotype M6. "Orange": genotype M9. "Control": environment without fish kairomone exposure. "Fish": environment with fish kairomne exposure (predation risk) and three were processed for sequencing (see below). The experiment lasted for 14 days for each experimental individual to assess the long-term effect of fish kairomones on gene expression level in D. galeata.

| RNA isolation and preparation
Appropriate amounts of RNA were not available from single individuals, and hence, we used pools of experimental individuals. Similar pooling approaches have been used in other Daphnia differential gene expression studies (Rozenberg et al., 2015;Hales et al., 2017;Herrmann et al., 2018;Huylmans et al., 2016;Ravindran et al., 2019;Orsini et al., 2018). Because embryos have high amounts of RNA as well and are able to sense predator cues, we were careful to control for egg developmental stage in the brood pouch (Weiss et al., 2016).
Only experimental females bearing freshly released eggs were pooled, resulting in a minor difference in age and experimental time as some experimental females had been pooled a day later. The advantage of sampling females in their intermolt stage (egg-bearing) is to ensure a stable gene expression (Altshuler et al., 2015). Total RNA was extracted from pools of 20 egg-bearing adults after homogenizing in RNAmagic (Bio-Budget technologies, Krefeld, Germany) for 5 min with a disposable pestle and a battery-operated homogenizer.
Samples were stored at -80°C until RNA isolation. Chloroform

| RNA-seq library construction and sequencing
Library construction and sequencing was identical for all samples and was performed by the company Macrogen (Seoul, South Korea).
RNA-seq libraries were constructed using Illumina TruSeq library kits. Illumina HiSeq4000 (San Diego, CA, USA) instrument was used for paired-end sequencing with 101-bp read length resulting in 48-79 million reads per library.
The option "strata" was used to output only the highest mapping scores for any given read and thus the uniquely mapped reads. The quality of filtering and mapping reads was verified with QualiMap v.2.2.1 (Okonechnikov et al., 2016). Subsequently, the htseq-count python script implemented in HTSeq v.0.9.1 was used to quantify the number of reads mapped to each transcript (Anders et al., 2015).

| Differential gene expression analysis
Differential gene expression analysis was performed in the R environment v.3.4.2 (R Core Team, 2017) with the R package "DESeq2" v.1.18.1 (Love et al., 2014) implemented in Bioconductor v.3.6 (Gentleman et al., 2004). The calculation was based on normalized read counts per environment (control & fish) using negative binomial generalized linear models. Prior to the analysis, all transcripts with a read count lower than 12 across all libraries were excluded.
Results were filtered post hoc by an adjusted p-value (padj < 0.05) (Benjamini & Hochberg, 1995) to reduce the false discovery rate (FDR) and filtered for a fold change ≥ 2. Differentially expressed transcripts (DETs) were binned into four groups: <twofold, two-to fourfold, four-to sixfold and > sixfold difference in expression. The three biological replicates were checked for homogeneity by principal component analysis (PCA). A differential expression analysis of genes between environments, between genotypes, and between environments within each genotype was done. In addition, a twofactor analysis was applied to investigate genotype-environment interactions (GxE). PCA plots were created in R with "ggplot2" v.2.2.1 (Wickham, 2009). The web tool jvenn (Bardou et al., 2014) was used to visualize the number of shared transcripts between groups.

| Gene coexpression network analysis
Variance-stabilized read counts obtained from the previous "DESeq2"-analysis were used in the coexpression analysis. First, an automatic, signed weighted, single gene coexpression network was constructed via an adjacency matrix in the R environment v.3.2.3 with the R package "WGCNA" v.1.61 (Langfelder & Horvath, 2008). Second, gene coexpression modules-clusters of highly interconnected genes-were identified based on the topological overlap matrices (TOM) with a soft cut-off threshold of 14 in "WGCNA". Module eigengenes (ME)-representing the average gene expression of their module-were calculated and used to investigate their relationship with other modules as well as external information (predation risk and genotype). ME-trait correlations were calculated to identify transcripts of interest with correlation values of> 0.5 or < −0.5. Finally, hubgenes-defined as the most interconnected genes per module-were identified to gain insight into the biological role of a gene coexpression module.

| Gene set enrichment analysis (GSEA)
To identify potential function of differentially expressed and coexpressed transcripts, we assigned Gene Ontology (GO) annotations using the reference transcriptome of D. galeata (Huylmans et al., 2016). To shed light on the biological importance of transcripts of interest, we performed a gene set enrichment analysis in R with the package "topGO" v.2.30.0 (Alexa & Rahnenfuhrer, 2016). The default algorithm "weight01" was used taking the hierarchy of GO terms into account, which results in fewer false positive results (Alexa & Rahnenfuhrer, 2016). Given that, a multiple testing correction after the Fisher's exact test was not applied (Timmermans et al., 2009).
GO terms of the three GO categories "Molecular Function" (MF), "Biological Process" (BP), and "Cellular Compounds" (CC) with a pvalue < 0.05 were considered significant.

| Orthology analysis
OrthoMCL cluster information from the reference transcriptome of D. galeata (Huylmans et al., 2016) was used to enhance our understanding of the functional roles of the transcripts of interest and their evolutionary history. "OrthoMCL" is a tool to identify clusters of homologous sequences in multiple species, that is, orthologs.
When assuming that orthologs are functionally conserved, known functions of orthologs in one species can be used to assign putative functions to sequences from other species in the same orthologous group (Li et al., 2003). These OrthoMCL clusters were originally build based on data for three Daphnia species (D. galeata, D. pulex and D. magna) and two insect species (Drosophila melanogaster and Nasonia vitripennis). Further details about the genome versions and annotations are available in the original publication by (Huylmans et al., 2016).
We analyzed the OrthoMCL clusters containing transcripts of interests by counting how many orthologs from other species where comprised in these clusters. Clusters were grouped into different categories: Daphnia galeata specific, Daphnia galeata plus one of the other Daphnia species (D. magna or D. pulex), Daphnia specific, and "all species" for those containing at least one transcript for each of the reference species (three Daphnia and two insect species). This allowed to measure how conserved our signal was and whether the response to predator risk was affecting transcripts specific to this particular species.

| RNA-seq data quality
RNA samples passed all quality steps before RNA sequencing. All 12 samples were successfully sequenced, resulting in 48.2 to 79.2 million reads of 101 bp length. After trimming and quality control, ~90% of trimmed reads were kept for further analysis. An average of 88% of these trimmed reads were uniquely mapped to the D. galeata reference transcriptome. After the filtering process, the full dataset used for further analysis comprised a total of 32,903 transcripts.

| Differential gene expression analysis
Before subsequent analysis, all transcripts with a read count lower than 12 across all libraries were excluded. 23,982 transcripts remained for both genotypes: 21,740 transcripts for genotype M6 and 21,813 for genotype M9.
A principal component analysis (PCA) was performed to visualize the grouping of read counts and to help identify batch effects.
The first principal component (PC 1) explained 83% of the variance between genotypes, revealing no clear clustering of read counts per environment (Figure 2a). PC 2 explained just 10% of the variance, which seems to be related to the variance between replicates.
Separate plots per genotype improved the visualization of replicate and environmental differences (Appendix S1 -Dryad repository) but did not indicate an evident clustering by environment either.
The differential expression analysis considering both genotypes in the two-factor analysis revealed no differentially expressed transcripts (DETs) between environmental groups, but a total of 5,283 DETs between genotypes (up: 2,228 (42%), down: 3,055 (58%); Figure 2b, Figure 3). Because of the strong genotype effect, the genotypes were analyzed separately in a one-factor analysis (Table 1A).

| Gene coexpression network analysis
The single network analysis clustered the expressed transcripts into 16 gene coexpression modules (Appendix S3 -Dryad repository, Table 2). A total of eight modules correlated with the external traits predation risk (PR) or genotype (G) with a correlation coef-ficient> 0.5 or < −0.5 (Figure 4). For readability, these modules of interest will be referred to as "module color_trait" (e.g., "blue_G").
The most highly interconnected gene within a gene coexpression module (hubgene) was identified for each module (Table 2).

| Orthology analysis
GO term and orthoMCL cluster information was available for a total of 9,172 D. galeata transcripts. Out of the 867 transcripts in the predation risk-related set, 600 were assigned to an orthology cluster. Predation risk-related transcripts of interest were distributed among 563 or-thoMCL clusters (2,131 D. galeata transcripts; GO annotation n = 224).
Most of these orthoMCL clusters comprise orthologs for all three Daphnia species (Figure 5), hinting at a common Daphnia response.

From an ecological point of view, predator-induced responses in
Daphnia have been studied extensively. In the past years, few studies addressed the link between such ecological traits and the underlying genetic pathways (Rozenberg et al., 2015;Hales et al., 2017;Orsini et al., 2018). Similar trends in life-history shifts after exposure to predator kairomones have been observed across Daphnia species showing, for example, the predominant trend of early maturation and a decreased body size under vertebrate predation risk (e.g., Riessen, 1999). This is why we chose specifically clonal lines with

| Insights from differential gene expression analysis-similar transcripts, differentially expressed?
In contrast to our expectations, the differential gene expression analysis revealed only a moderate number of differentially expressed with previous studies in our group (Ravindran et al., 2019). Instead, the apparent genotype-specific response in our study might be explained by the phenotypic divergence between the two studied clones.
One explanation for the low number of DETs concerning predation risk (environment) compared to genotype-specific differences is that life-history changes could only marginally correlate with gene expression. The D. galeata genotypes used here only displayed shifts in life history, whereas other Daphnia species show additional adaptations of morphology and behavior that could be caused by or correlated to much stronger differential gene expression, for example, neck-teeth induction that was linked to 230 differentially expressed genes in D. pulex (Rozenberg et al., 2015).
The moderate number of DETs found under predation risk could be explained by other causes than gene expression as well; additional post-translational processes, such as miRNA-mediated regulation or increased degradation (Schwarzenberger et al., 2009), might play a role. Epigenetic modifications, such as cytosine methylation, can be another explanation for our findings. Asselman et al (2015) showed that epigenetic effects might be important in Daphnia in response to environmental changes, such as shifts in predation regimes.
There were three reasons why we expected more pronounced and cumulative changes in differential gene expression in the third experimental generation. First, the chosen D. galeata genotypes displayed strong shifts in life-history traits after three generations of fish kairomone exposure (Tams et al., 2018). Second, the effect of kairomone exposure is expected to be cumulative and to increase over the course of multiple generations; for example, D. cucullata displays the largest helmets when exposed to kairomones from Leptodora kindtii (an invertebrate predator) and Chaoborus for two generations compared to the first generation (Agrawal et al., 1999).
Third, transgenerational plasticity was described in D. ambigua (Hales et al., 2017); genes were significantly differentially expressed after one generation of fish kairomone exposure (n = 48 DEGs) and without kairomone exposure after the second (n = 223 DEGs) and third

| Insights from gene coexpression and gene set enrichment analysis (GSEA)-different transcripts, different functions?
Although different transcripts were identified in the differential expression analysis, the combined approach revealed their similarity of biological functions. In brief, our gene coexpression and gene set enrichment analysis revealed digestion-and growth-related enriched GO terms. These are interesting because predatorinduced responses in Daphnia include changes in body size (e.g., Tams et al., 2018) and morphological modifications (e.g., Laforsch & Tollrian, 2004). Such modifications require energy allocated from nutrients; consequently, digestive enzymes like peptidases have been shown to be important for juvenile growth rate in D. magna (Schwarzenberger et al., 2012).
Transcripts of the "salmon_PR" gene coexpression module in D. galeata were significantly enriched for "serine-type endopeptidase activity," the most important digestive protease in the gut of D. magna (Agrawal et al., 2005). The exposure to predator kairomones for one generation in D. ambigua-a species from the D. pulex-complex more closely related to D. galeata than D. magna (Cornetti et al., 2019;Adamowicz et al., 2009)-led to an upregulation of genes related to digestive functions (Hales et al., 2017).

Cyanobacterial protease inhibitors cause considerable damage
to Daphnia populations by inhibiting the gut proteases, impairing their digestion (Schwarzenberger et al., 2010). These studies concord with our results suggesting that an increase in "serine-type endopeptidase activity" leads to improved digestion and feeding efficiency that is necessary for the resource allocation that comes with shifts in life history, such as producing a greater number of offspring.
The GO term "structural constituent of cuticle" was significantly enriched in both genotypes suggesting that even though there was no overlap in the affected transcripts, similar functions were affected. The "structural constituent of cuticle" was enriched in D. pulex exposed to Chaoborus kairomones Rozenberg et al., 2015) and related to remodeling of the cuticle.
Furthermore, it was also enriched in the proteomic response of D. magna to Triops cancriformis (Otte et al., 2015) and is thought to be related to changes in carapace morphology as well as the formation of ultrastructural defenses of the cuticle (Rabus et al., 2013).
Genes related to body remodeling and activation of cuticle proteins were enriched for D. magna exposed to vertebrate and invertebrate predator kairomones (Orsini et al., 2018). Furthermore, for D. magna, D. pulex, and D. cucullata, not only visible morphology changes but also fortification of the carapace in the presence of invertebrate predator kairomones has been recorded (Laforsch & Tollrian, 2004;Rabus et al., 2013;Kruppert et al., 2017). The investigated D. galeata genotypes did not display pronounced morphological defenses under vertebrate predation risk, but changes in body size and symmetry especially with regard to head shape (Tams et al., 2018).
Hence, our results indicate that ultrastructural defenses could also be present in D. galeata under vertebrate predation risk and could be an interesting field of investigation.
Altogether, cuticle-associated proteins seem to play an essential role in the response to vertebrate or invertebrate predators in Daphnia. DETs found in genotype M6 showed the possible involvement of "metallocarboxypeptidase activity," known to be involved in the stress response to copper in D. pulex (Chain et al., 2019).
Interestingly, "chitin metabolic process", "proteolysis", "structural constituent of cuticle", "chitin binding", "serine-type endopeptidase", and "metallopeptidase activity" were all found to be enriched in a gene expression analysis during the molt cycle in the marine copepod Calanus finmarchicus (Tarrant et al., 2014). Since Daphnia need to shed their rigid carapace in order to grow, molting is directly related to changes in body size. Another analysis of D. magna exposed to Triops cancriformis kairomones revealed the role of proteins related to the cuticle, muscular system, energy metabolism and regulatory proteins that may be involved in morphological carapace defenses and changes in resource allocation (Otte et al., 2014). In conclusion, a number of biological functions hypothesized to be involved in kairomone response could be confirmed, for example transcripts related to body remodeling and growth.
It is worthwhile to mention that some biologically interesting gene functions were only found with the help of the gene coexpression network analysis and would have been overlooked with only a differential expression analysis. For example, the GO term "growth factor activity" occurred in both "red_PR" and "tan_PR" modules, which correlated negatively with fish kairomone exposure and comprising transcripts not identified as DETs. Nevertheless, they could be extremely important for life-history changes and might be directly related to changes in somatic growth rate and body size.

| Insights from orthology analysis-homologous sequences (common ancestor), evolutionary conserved?
Using orthology-orthologs being homologous sequences that differ because of speciation events-is the most popular strategy to derive functional similarity of sequences (Pearson, 2013).
Here, we chose this approach to gain insight into the evolutionary conservation of transcripts of interest. The results obtained with the OrthoMCL approach provide support for the hypothesis that predator-induced plasticity might be evolutionary conserved in Daphnia. However, several strategies might have evolved over time to cope with or adapt to predation risk. This hypothesis seems likely since Daphnia have the ability to rapidly adapt to local predator regimes (Declerck & De Meester, 2003) and our study provides elements supporting it. First, the differential gene expression analysis revealed genotype-specific molecular responses to predation risk for genotypes originating from the same population. Second, the involved transcripts have similar functions relating to life-history changes induced by predation risk, but different transcripts were involved in the predator-induced response for each genotype. This concurs with the suggestion of niche-specific adaptation in D. magna due to the genotype-and condition-specific transcriptional response to environmental changes of biotic and abiotic factors (Orsini et al., 2018). Their gene coexpression analysis revealed that genes of interest were crustacean related, meaning that the conservation of genes did not exceed the level of crustaceans (Orsini et al., 2018). Further insights into evolutionary conservation of differentially and/or coexpressed transcripts linked to phenotypic traits are available for modern and ancient (resurrected) D. pulicaria exposed to different phosphorous regimes (Frisch et al., 2020). With a different, yet similar approach this recent study reveals the importance of a holistic approach to tackle the question: What is the molecular basis of phenotypic responses to environmental changes?
Gene expression analyses to uncover the molecular basis of predator-induced phenotypic changes have focused so far on single species, and invertebrate predators (e.g., Rozenberg et al., 2015;Orsini et al., 2018). Studies conducted in similar conditions on other Daphnia species, that is, long-term exposure to vertebrate predation risk, are lacking and prevent us from drawing conclusions about a general Daphnia response to fish predation risk. In the future, simultaneous exposure of several species to kairomones, and the coupling of phenotyping and gene expression would help to address the question of a conserved response.

| CON CLUS ION
In summary, the aim of this study was to characterize the genetic basis for the predator-induced response of the freshwater grazer D. galeata. Our hypothesis that clonal lines present a common pred-

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R S CO NTR I B UTI O N S
VT, JHN, and MC: Study design and writing. JHN, AE, and VT: Laboratory work. JHN and VT: Gene expression and gene network analysis. All authors: Final approval.