Skip to main content

Comparative transcriptome analysis of two Daphnia galeata genotypes displaying contrasting phenotypic variation induced by fish kairomones in the same environment of the Han River, Korea

Abstract

Background

Phenotypic plasticity is a crucial adaptive mechanism that enables organisms to modify their traits in response to changes in their environment. Predator-induced defenses are an example of phenotypic plasticity observed across a wide range of organisms, from single-celled organisms to vertebrates. In addition to morphology and behavior, these responses also affect life-history traits. The crustacean Daphnia galeata is a suitable model organism for studying predator-induced defenses, as it exhibits life-history traits changes under predation risk. To get a better overview of their phenotypic plasticity under predation stress, we conducted RNA sequencing on the transcriptomes of two Korean Daphnia galeata genotypes, KE1, and KB11, collected in the same environment.

Results

When exposed to fish kairomones, the two genotypes exhibited phenotypic variations related to reproduction and growth, with opposite patterns in growth-related phenotypic variation. From both genotypes, a total of 135,611 unigenes were analyzed, of which 194 differentially expressed transcripts (DETs) were shared among the two genotypes under predation stress, which showed consistent, or inconsistent expression patterns in both genotypes. Prominent DETs were related to digestion and reproduction and consistently up-regulated in both genotypes, thus associated with changes in life-history traits. Among the inconsistent DETs, transcripts encode vinculin (VINC) and protein obstructor-E (OBST-E), which are associated with growth; these may explain the differences in life-history traits between the two genotypes. In addition, genotype-specific DETs could explain the variation in growth-related life-history traits between genotypes, and could be associated with the increased body length of genotype KE1.

Conclusions

The current study allows for a better understanding of the adaptation mechanisms related to reproduction and growth of two Korean D. galeata genotypes induced by predation stress. However, further research is necessary to better understand the specific mechanisms by which the uncovered DETs are related with the observed phenotypic variation in each genotype. In the future, we aim to unravel the precise adaptive mechanisms underlying predator-induced responses.

Peer Review reports

Background

Organisms are affected by a wide range of biotic and abiotic environmental factors, interacting, and adapting to environmental changes by phenotypic plasticity, evidenced through behavioral or physiological changes [1, 2]. Phenotypic plasticity refers to the capacity of genotypes to generate a range of phenotypes in response to environmental conditions, facilitating the survival and reproduction of organisms in diverse habitats [1, 2]. Phenotypic variation within a population is critical to its persistence, as less variation increases the risk of extinction [3, 4]. Genetic variation due to environmental change, gene flow, positive selection, mutation, and recombination can lead to loss or gain of phenotypic plasticity [5,6,7,8]. In addition, as phenotypic plasticity allows to locally adapt to a changed environment, if a population produces a phenotype more adapted to its environment, then the genotype of that population will contribute more to the genetic makeup of the entire population [9].

Among all environmental factors, predation is an important biological factor that maintains species diversity, organizes entire communities, and drives natural selection in populations [10, 11]. Vertebrate and invertebrate aquatic predators emit kairomones into the surrounding water, which can be detected by prey, eliciting predator-specific responses to reduce vulnerability [12, 13]. Adaptive strategies in response to predators, a typical example of phenotypic plasticity, include life-history changes, behavioral plasticity, and changes in morphological and physiological traits [14], and have been extensively studied across diverse Daphnia species [15, 16].

Daphnia species are small branchiopod crustaceans widely used as model organisms in ecology, evolution, and ecotoxicology [17, 18]. This genus primarily feeds on algae and is highly vulnerable to predation risk because it is the preferred food of many fish [17, 19]. In the genus Daphnia, broad alterations in morphology, behavior, and life-history traits have been observed in response to predation risk. Daphnia species have developed suitable survival and reproductive strategies in various aquatic environments and vegetation conditions, allowing them to inhabit regionally diverse underwater environments. The ability of Daphnia to locally adapt to different stressors has been demonstrated for fish as a vertebrate predator and Chaoborus as an invertebrate predator, for instance [20,21,22]. Such adaptations are associated with the maintenance of genetic diversity in each region and are related to intraspecific phenotypic variation. The distinct phenotypic responses observed among Daphnia populations may be due to unique colonization patterns, monopolization effects, and limited gene flow between populations, as well as significant genetic diversity [23]. In addition, one population of a Daphnia species can have individuals showing different genotypes; these genetic variations act as adaptive strategies to improve their ability to cope with environmental changes [9, 24]. Different genotypes have their own characteristics and survival strategies, allowing Daphnia populations to spread and adapt within their ecological range. Thus, studying genetic diversity within populations is essential for understanding the biodiversity and adaptive strategies of local ecosystems, as it plays an important role in developing different adaptive strategies and adapting to specific environmental conditions.

Despite previous reports on clonal variation in Daphnia [25, 26], previous studies have primarily drawn conclusions based on a single genotype (clonal line) rather than in ≥ 1 genotypes per population. Although intra-population variation is typically associated with population maintenance in response to predation pressure, the significance of such variation has rarely been studied. Previous studies have shown that when 16 D. magna genotypes from four different habitats were exposed to fish kairomones, response traits could vary between genotypes [20]. Moreover, Tams et al. found differences in life-history traits across 24 D. galeata genotypes within the same population and between populations when exposed to fish kairomones [9].

D. galeata is a relatively large-sized zooplankton found throughout the northern hemisphere and identified as a dominant plankton species in the Han River, Korea [27, 28]. It is also a member of the Daphnia longispina complex, which includes D. longispina, D. cucullata, D. hyalina, and D. galeata, and a widely distributed species, whose genetic diversity has been investigated in several countries in Asia, Europe, and North America [28,29,30]. Although this species does not exhibit substantial morphological changes or diel vertical migration in response to vertebrate predator cues, it exhibits marked phenotypic variation in life-history traits when facing predation risk [31]. Accordingly, D. galeata shows an adaptive strategy to maintain generations by producing more offspring when exposed to fish kairomones [32]. In addition, Daphnia can increase its body size in the presence of predators as a response to predator-induced stress [33]. A larger body size can provide advantages such as increased fecundity and improved competitive ability, which can help Daphnia thrive in predator-rich environments [33]. Conversely, Daphnia can also decrease its body size in response to predators [34, 35]. This adaptive response allows them to better avoid being detected and eaten by the predator. A smaller body size reduces visibility, making it more difficult for the predator to see and capture the prey. Previous studies have also shown that when Daphnia is exposed to fish kairomones, energy resources are preferentially allocated to somatic growth, increasing the rate at which neonates reach adulthood, thereby advancing their reproductive age [36]. These results show that the somatic growth rate is also related to reproduction [32, 37].

Extensive research has been conducted on phenotypic plasticity responses to predation stress in Daphnia, with efforts made to establish connections between these responses and genome-wide expression patterns. For example, Tams et al. performed transcriptional profiling of D. galeata exposed to fish kairomones and identified differentially expressed genes related to digestion and growth [32]. In addition, under vertebrate predation risk, D. ambigua also showed expression changes in genes involved in digestion, reproduction, and exoskeleton structure [38]. Another study identified up-regulated genes such as those encoding for the cuticle, vitellogenin, and others in response to invertebrate predation risk in D. pulex [22, 36]. However, the genetic basis of the appropriate functional responses associated with phenotypic plasticity in Daphnia species to predation risk is not well understood. Furthermore, it is not clear how different genotypes inhabiting the same environment exhibit different responses in their life-history traits under predation stress from a genotype-specific perspective.

The present study aimed to compare the phenotypic variation and associated transcriptional profiles of two D. galeata genotypes exposed to fish kairomones, which show different life-history traits despite inhabiting the same environment in the Han River, Korea. We applied a transcriptomic approach (RNA-seq) to compare the gene expression differences related to phenotypic variations between the two genotypes, followed by differential expression (DE), Gene Ontology (GO), and Kyoto Gene and Genome Encyclopedia (KEGG) pathway analyses. Additionally, network analysis was performed to infer the functions of the detected differentially expressed transcripts (DETs). Network-based research offers a structured approach to comprehensively analyze complex datasets and facilitate a holistic understanding of multiple interacting groups. The present results offer the foundation for a better understanding of the adaptation mechanisms related to the reproduction and growth of two Korean D. galeata genotypes under predation stress.

Results

Life-history traits

To compare the changes in life-history traits under predator stress in two D. galeata genotypes, six life-history traits related to reproduction and growth were measured. Reaction norms of life-history traits show intraspecific variation of life-history traits between the two D. galeata genotypes (Fig. 1). A one-way ANOVA showed that “fish kairomones” had a significant effect (P < 0.01) on the six life-history traits in both genotypes (Fig. 1). Both genotypes exhibited significant changes (P < 0.001) in life-history traits associated with reproduction, specifically decreasing the “Age at First Reproduction” in the presence of fish kairomones (Fig. 1A). In contrast, the “Number of Offspring First Brood,” “Total Number of Broods,” and “Total Number of Offspring” significantly increased (P < 0.05) in both genotypes in the presence of fish kairomones (Fig. 1B–D). Both genotypes exhibited similar patterns in life-history traits related to reproduction under fish kairomone exposure, which may be an adaptive strategy to produce more offspring and maintain generations by advancing reproductive age. However, the two genotypes exhibited contrasting patterns in terms of life-history traits associated with growth. In the life-history traits related to growth, “Somatic Growth Rate” of the genotype KE1 significantly increased (P < 0.001) in the presence of fish kairomones, and as a result, “Body Length” significantly increased (P < 0.01; Fig. 1E and F). Conversely, the “Somatic Growth Rate” of the genotype KB11 did not differ between fish kairomone-exposed and non-exposed conditions (Fig. 1E). However, “Body Length” was significantly decreased (P < 0.01) in the genotype KB11 (Fig. 1F). In summary, genotype KE1 exhibited an adaptive strategy by increasing the “Somatic Growth Rate” and “Body Length”, thereby promoting faster growth and then investing in high reproductive effort. In contrast, genotype KB11 decreased the “Body Length” to better avoid being detected and eaten by the predator. Hence, it can be observed that each genotype has a distinct adaptation strategy to cope with predation by fish.

Fig. 1
figure 1

Box plots of reaction norms for selected life-history traits of two D. galeata genotypes (median +/− SD). A, Age at First Reproduction. B, Number of Offspring First Brood. C, Total Number of Broods. D, Total Number of Offspring. E, Somatic Growth Rate. F, Body Length. Control.KE1 indicates the control group (without fish kairomones) of genotype KE1 and fish.KE1 indicates the experimental group (fish kairomone-exposed) of genotype KE1. Control.KB11 indicates the control group (without fish kairomones) of genotype KB11 and fish.KB11 indicates the experimental group (fish kairomone-exposed) of genotype KB11. Stars indicate significant differences (P.adjust *** <0.001, ** <0.01, * <0.05, N.S., no significance) between fish kairomone-exposed and non-exposed in both genotypes

Datasets and de novo transcriptome assembly

RNA-seq was performed on the control and experimental groups of the two D. galeata genotypes. A total of 704,814,567 raw reads were generated from RNA-seq, of which 97.86% (689,735,869 reads) were retained after removing low-quality reads and trimming. Table S1 provides a summary of the RNA-seq data. A de novo transcriptome assembly generated a total of 186,084 transcripts, after removing redundant sequences, the longest transcripts from the CD-Hit clustering results were selected as unigenes, resulting in a total of 135,611 unigenes. Among them, the longest had a length of 28,597 bp while the shortest was 181 bp long. The mean GC content of unigenes was 40.7% (Table S2). The BUSCO estimation of completeness of the assembled transcripts showed that most orthologs among arthropods were represented in the assembled transcriptomes. In each dataset, 99.1% of complete orthologs belonging to the phylum Arthropoda were detected out of the 1013 BUSCO groups examined (Table S3). This includes both complete single copy and duplicates (putative paralogs or complete genes with multiple copies).

Identification of coding regions and unigene functional annotation

A total of 135,611 unigenes were employed as query sequences for functional annotation, with TransDecoder identifying 89,738 unigenes having the longest open reading frames (ORFs) within the unigene datasets. Only ORFs that were at least 100 amino acids long were considered; sequences shorter than 100 residues were excluded. To identify ORFs with homology to known proteins, these transcripts were scanned against the UniProt and Pfam databases. The annotation result and corresponding accession numbers for all unigenes are provided in Table S4. Among 89,738 unigenes, a BLASTp-based homology search for D. galeata resulted in the annotation of 56,306 (62.74%) unigenes (Table S5). Similarly, about 75% and 66% of unigenes were mapped to the Pfam and eggNOG databases, respectively. In contrast, only 1,405 unigenes had KEGG annotations (Table S5). Overall, a total of 74,082 unigenes were annotated using various databases.

DE and GO enrichment analysis

Pearson correlation coefficient and principal component analysis (PCA) were performed to confirm the deviation of triplicated samples by visualizing differences in transcript expression (Fig. 2). The pearson correlation coefficient showed a strong correlation between samples for each genotype (Fig. 2A), and the PCA clearly clustered samples according to genotypes. The first principal component (PC 1) explained 50% of the variance between genotypes, and PC 2 explained 26% of the variance, an effect possible related to variance between replicates (Fig. 2B). Individually, each genotype revealed obvious clustering between control and experimental group (Figure S1). Therefore, we confirmed the consistency between triplicate samples, so the triplicate data was used for DE analysis.

Fig. 2
figure 2

A, Pairwise correlation of all biological replicated RNA-seq samples in D. galeata. B, PCA plot of gene expression in D. galeata RNA-seq samples. Red, control group (without fish kairomones); Cyan, experimental group (fish kairomone-exposed). Circles, genotype KE1; Triangles, genotype KB11

To understand the association between phenotypic variation and genotype-specific responses to fish kairomones, we selected and described DETs based on three criteria: (a) shared DETs that showed consistent expression profiles across both genotypes; (b) shared DETs that showed inconsistent expression profiles across both genotypes and; (c) DETs uniquely identified in each genotype (Table 1; Fig. 3). The complete list of DETs from both genotypes with detailed annotations is provided in Table S6. First, 194 DETs were shared among the two genotypes, 85 exhibited consistent expression profiles (e.g., up or downregulation) in both genotypes (Fig. 3 and Table S7). Consistently up-regulated DETs were homeobox protein cut (CT), dipeptidase 1 (DPEP1), vitellogenin (VG), and others (Table 1 and Table S7). In contrast, consistently down-regulated DETs were vitellogenin (VG), perlucin, and others (Table 1 and Table S7). Second, the remaining 109 DETs showed inconsistent expression (i.e., these transcripts exhibited contradictory expression patterns in both genotypes; Fig. 3 and Table S7). For example, vinculin (VINC) and protein obstructor-E (OBST-E) were up-regulated in KE1 but down-regulated in KB11 (Table 1 and Table S7). Similarly, the protein flightless-1 (FLIL) and dynamin (SHI) were down-regulated in KE1 but up-regulated in KB11 (Table 1 and Table S7). Finally, to associate the uniquely found DETs in each genotype with observed phenotypic variation, we compiled a list of 896 and 801 unique DETs in KE1 and KB11, respectively (Figure S2). To accomplish this, we initially excluded the DETs shared by both genotypes. Subsequently, we excluded DETs from each genotype that were annotated with the same Uniprot ID but displayed inconsistent expression patterns. In KE1, up-regulated DETs related to growth, including matrix metalloproteinase-2 (MMP2), dopamine D2-like receptor (DOP2R), and others were found; this was not the case in KB11 (Table 1 and Table S6). Conversely, down-regulated DETs related to growth, such as paxillin (PXN), muscle Lim protein Mlp84B (MLP84B), and others were found in KB11 but not in KE1 (Table 1 and Table S6).

Table 1 The list of consistent, inconsistent, and unique DETs associated with fish kairomones in each genotype
Fig. 3
figure 3

Venn diagram showing the number of DETs identified from each D. galeata genotype. a) Shared DETs showing consistent expression profiles across genotypes; b) Shared DETs showing inconsistent expression profiles across genotypes and; c) DETs uniquely identified in each genotype. ↑ and ↓ indicate up- and down-regulation, respectively

Next, the enriched GO categories related to the biological process (BP), molecular function (MF), and cellular component (CC) were identified for the DETs. In KE1, GO terms related to development, such as “system development,” “cellular developmental process,” “cell development,” and “animal organ development” were among the top 10 enriched biological processes (Fig. 4A). In addition, 32 enriched GO BP were related to signaling and included regulation of signaling, cell-cell signaling, anterograde trans-synaptic signaling, trans-synaptic signaling, G protein-coupled receptor signaling pathway, and tachykinin receptor signaling pathway (Table S8). Most enriched GO MF were associated with binding of various molecular compounds (e.g., cation, metal ion, anion, small molecule, carbohydrate, etc.) or hydrolase and transmembrane transporter activities in KE1 (Fig. 4A). The GO CC showing the highest enrichment were mainly related to membranes (e.g., plasma membrane, integral component of plasma, etc.), cell junction, synapse, cell projection, cytosol, and neuron projection in KE1 (Fig. 4A).

Similarly, GO terms related to development were the enriched BP in the KB11 (Fig. 4B). In addition, 27 enriched GO BP related to signaling such as signal transduction, regulation of signaling, cell surface receptor signaling pathway, cell-cell signaling, anterograde trans-synaptic signaling, trans-synaptic signaling, hippo signaling, notch signaling pathway, and regulation of G protein-coupled receptor signaling pathway were detected in KB11 (Table S8). Interestingly, the tachykinin receptor signaling pathway was enriched only in KE1, while hippo signaling and notch signaling pathways were enriched only in KB11 (Table S8). In addition, the GO BP related to growth such as somatic muscle development, developmental growth, and cuticle development were enriched only in KE1, where “Somatic Growth Rate” and “Body Length” were significantly increased under fish kairomone exposure (Table S8). The GO MF of KB11 were found to be similar to those of KE1, with the most enriched categories associated with the binding of various molecular compounds such as cations, metal ions, anions, small molecules, carbohydrates, and more (Fig. 4B). We identified three down-regulated DETs, MLP84B, myosin heavy chain, muscle (MHC), and pupal cuticle protein Edg-78E (EDG78E), only in genotype KB11. These three DETs have been involved in MFs related to structural constituents of muscle or/and structural constituent of pupal chitin-based cuticle, and likely associated with the decrease in “Body Length” observed under fish kairomone exposure. Most enriched GO CC were related to membranes, cytosol, chromosome, nuclear protein-containing complex, nuclear lumen, cell junction, nucleoplasm, and cytoskeleton in KB11 (Fig. 4B).

Fig. 4
figure 4

Gene ontology (GO) enrichment of DETs. The provided results display the top 10 processes within each GO category, encompassing biological processes (BP), molecular functions (MF), and cellular components (CC). A, genotype KE1. B, genotype KB11. Percentage (%) indicates the ratio of assigned genes to each GO term

Enrichment analyses of KEGG pathways of DETs

The DETs were used for searching pathways that contribute to the alterations in life-history traits. Enrichment analysis of the pathway was conducted using the KEGG database. Table 2 provides an overview of the top 10 enriched KEGG pathways in KE1 and KB11, respectively; a complete list of pathways is provided in Table S9. Enriched categories with a corrected p-value < 0.05 were selected. The KEGG enrichment analysis identified 24 significant pathways of the DETs for KE1. Most were related to metabolic pathways, amino sugar and nucleotide sugar metabolism, AGE-RAGE signaling pathway in diabetic complications, and fatty acid metabolism. In addition, various metabolic (e.g., arachidonic acid, fructose and mannose, taurine, hypotaurine, etc.) and biosynthetic (e.g., fatty acid, steroid, and folate) pathways were enriched in KE1. Similarly, KEGG enrichment analysis identified 30 significant pathways for KB11, most related to autophagy - animal, protein processing in endoplasmic reticulum, metabolic pathways, AGE-RAGE signaling pathway, and amino sugar and nucleotide sugar metabolism in diabetic complications. Interestingly, four pathways, Wnt signaling pathway, FoxO signaling pathway, ECM-recepter interaction, and MAPK signaling pathway - fly were identified only in KB11.

Table 2 Top 10 enriched KEGG pathways for DETs of two D. galeata genotypes identified by RNA-seq. Only enriched categories with a p-value < 0.05 were selected

Interaction network of DETs

There has been an increasing interest in utilizing biological networks to explore and understand significant biological phenomena [39]. Network-based investigations provide a comprehensive framework for gaining a holistic understanding of data, encompassing multiple interacting groups. The potential of community-based network analysis enables researchers to explore and compare a diverse range of proteins and their interactions within the identified modules or communities of the network [40]. Therefore, we attempted to identify the modular structure of the interaction network between DETs, thus dividing the network into functional six communities (clusters; Fig. 5).

To predict interactions between DETs, a non-redundant list of 2,119 DETs was compiled from both genotypes and mapped to STRING-DB using the common fruit fly (Drosophila melanogaster) as reference organism. Using the UniProt database to annotate unigenes, Drosophila melanogaster, exhibited the highest number of annotated unigenes among arthropods, leading to a subsequent network analysis using this species. Of these 2,119 DETs, only 338 were mapped to the fruit fly proteome. The resultant interaction network was subjected to community (cluster) detection to identify the functional modules. After removing clusters with < 10 nodes, the interaction network was divided into six main clusters comprising 257 nodes (DETs), and 1,314 edges (interactions between nodes; Fig. 5 and Table S10). The number of nodes in each cluster ranged between 11 (cluster 6) and 73 (cluster 1). The node degree of each node in the network ranged between 1 and 29; nodes with a maximum node degree ≥ 20 were considered as hubs. These hub nodes were identified to encode a tyrosine-protein kinase Src42A (SRC42A: node degree = 29) and DNA-directed RNA polymerase II subunit RPB1 (RPII215: node degree = 22) in clusters 1 and 2, respectively. Enrichment analysis indicated that each of these clusters had specific functional roles in their respective BP (Table S10 and S11). Most enriched GO processes in the largest cluster 1 were related to chemical synaptic transmission, ion transport, nervous system process, or development (Fig. 5 and Table S10). The smallest cluster 6 was enriched in processes related to organization, oogenesis, or regulation of growth, including others (Fig. 5 and Table S10). Meanwhile, numerous signaling pathways such as hippo signaling pathway – fly, AGE-RAGE signaling pathway, and phosphatidylinositol signaling system were found in cluster 1 (Fig. 5 and Table S11). Among the 257 nodes, 60 transcripts were commonly expressed in both genotypes, of which 25 transcripts that showed inconsistent expression profiles within the same genotype were excluded from the analysis. The remaining 35 transcripts were used for further analysis; four and six transcripts displayed consistent up and downregulation patterns in both genotypes, respectively (Fig. 5 and Table S12). On the other hand, the remaining 25 transcripts exhibited inconsistent expression patterns in both genotypes (Table S12). The enrichment analysis of GO terms of DETs up-regulated in KE1 and down-regulated in KB11 were related to membrane, cuticle, and muscle (Table S12). This appears to be associated with increased “Body Length” upon exposure to fish kairomones in KE1.

Fig. 5
figure 5

Interaction network of DETs detected in both genotypes. Red and blue nodes represent up- and down-regulated transcripts, respectively. Nodes in the rectangle and triangle indicate the genotype KE1 and KB11, respectively. Green and cyan squares represent up- and down-regulated transcripts in both genotypes. Yellow diamonds represent transcripts showing conflicting expression patterns in both genotypes. Large nodes in clusters 1 and 2 represent the hubs. The top 3 enriched KEGG pathways (red text) and GO terms (green text) are indicated for each cluster

Discussion

Predator-induced responses in Daphnia have been extensively investigated [36, 38]. However, there is still a lack of understanding regarding the genetic mechanisms that regulate the appropriate functional responses associated with phenotypic plasticity in response to predation risk. Additionally, the variation in life-history traits of different genotypes to predation stress in the same environment remains poorly understood from a genotype-specific perspective. Prey species exhibit various adaptations in response to predation, including alterations in their life history, behavior, and morphological and physiological traits [14]. For example, the life history showed the direction of early or late maturation and increased or decreased body size under vertebrate predation risk [9, 41]. Therefore, we selected two genotypes showing contrasting life-history traits to understand the variation in phenotypic plasticity of different genotypes in response to predation stress in the same environment from a genotype-specific perspective. The two genotypes of D. galeata examined in the study were similar in terms of life-history traits related to reproduction, but with respect to growth, one genotype increased body size, while the other decreased body size under predation risk, indicating contrasting adaptive strategies. To investigate the association between phenotypic variation and differentially expressed transcripts in response to predator-induced cues, the two D. galeata genotypes were subjected to gene expression profiling after exposure to fish kairomones, which mimics predation risk. In addition, selected DETs of both genotypes were functionally annotated, and GO, KEGG pathway, and interaction network of DETs analyses were performed to infer functions associated with phenotypic variation. Our results revealed that phenotypic plasticity is relevant to reproduction and growth under predation stress in D. galeata.

Similar to previous studies [12, 20], our results showed altered life-history traits related to reproduction and growth in both genotypes under fish predator stress. The genotype KB11 exhibited defense strategies to maintain the generation by advancing reproductive age, producing more offspring, and evading predators by reducing body size. Previous studies have reported the occurrence of early maturation and reduced Daphnia size in response to vertebrate predators [34, 35]. The ecological advantage of earlier reproduction before predation is the maintenance of the population numbers [34, 42]. If most individuals within a population produce relatively few offspring in a fish environment, it will result in fewer offspring in the next cohort, which could threaten the persistence of the entire population. Therefore, it is presumed to represent a strategy to maintain the generation by producing more offspring. Additionally, exposure to predatory fish kairomones can increase the reproductive rate in Daphnia, which can further enhance their fitness and ability to cope with predation pressure [34, 42]. In KE1, life-history traits related to reproduction changed similarly to KB11 in the presence of fish kairomones (Fig. 1A–D). However, growth-related life-history traits such as “Somatic Growth Rate” and “Body Length” increased, showing opposite patterns to KB11 (Fig. 1E and F). Exposure to predatory fish kairomones can increase body length in Daphnia as an adaptive response to predator-induced stress [33]. When Daphnia perceive the presence of predatory fish, they activate certain physiological, and behavioral mechanisms to increase their chances of survival. One of these mechanisms involves allocating more energy toward growth and development, which leads to an increase in body length. This response is thought to be mediated by the activation of certain signaling pathways and gene expression changes, which ultimately result in the up-regulation of growth-related processes. In addition, previous studies suggested that Daphnia under the threat of fish might grow faster during juvenile stages and then invest in high reproductive effort and fast clutch release once they reach maturity [33]. Furthermore, previous studies suggested that vulnerability to fish predation in Daphnia could be increased by a larger size at first reproduction, which can be interpreted as an ecologic cost (community trade-off) [43] of the defensive reaction, deviating from the classical prediction of decreased body size as a response to risk of fish predation [44, 45].

DE analysis was performed to understand the genetic mechanisms regulating the appropriate functional response associated with phenotypic plasticity in two D. galeata genotypes that exhibit contrasting phenotypic plasticity under predation stress. In the study, among the 194 DETs commonly identified in both genotypes, 85 DETs (13 up-regulated and 72 down-regulated) were consistent (Fig. 3) and likely commonly expressed in response to fish kairomones, regardless of differences in contrasting life-history traits in the two genotypes; with an increasing number of collected D. galeata genotypes, the number of consistent DETs will probably decrease. Among the 85 DETs, there was an up-regulated DET encoding CT (Table S7), an important regulator of the stress response and plays an important role in protecting cells from stress-induced damage in various species including Drosophila [46]. Moreover, one of the consistently up-regulated transcripts encoded DPEP1 which is implicated in digestive functions (Table S7). Exposure to predator kairomones for one generation in D. ambigua led to an up-regulation of genes related to digestive functions [38]. In addition, it has also been demonstrated in D. magna that digestive enzymes such as peptidases affect juvenile growth rates [47], and that cyanobacterial protease inhibitors impair digestion by inhibiting gut proteases, causing significant damage to D. magna populations [48]. These studies corroborate our findings, suggesting that the activation of digestive functions improves digestion and feeding efficiency, which plays a vital role in resource allocation that comes with shifts in life history, such as producing a greater number of offspring. This is supported by the observed phenotypic variation, where both genotypes exhibited a higher number of offspring under predator stress. In addition, a transcript encoding vitellogenin (VG), a precursor to a major yolk protein, annotated with the Q868N5 in the UniProt database was up-regulated in both genotypes (Table S7). A previous study unveiled a set of DETs in D. pulex, where the up-regulated genes included vitellogenin genes, known for their significance in the response to invertebrate predation risk [36]. Increased vitellogenin expression seems to indicate early onset of vitellogenesis and increased fertility, and in the case of D. galeata, fish kairomone induction increases fertility, leading to the production of more offspring, which presumably requires a larger yolk pool. On the other hand, a transcript encoding vitellogenin-3 (VIT-3) annotated with the Q9N4J2 was down-regulated in both genotypes. This could be as it represents a trade-off between reproduction and survival. In the presence of a predator, it may be advantageous for Daphnia to prioritize survival over reproduction by allocating resources toward physiological processes that enhance their ability to escape or defend against predators, rather than investing in reproductive processes that require energy and resources. However, further studies are needed to understand the defense strategies mediated by the transcripts encoding the two vitellogenins, which exhibit opposite expression patterns.

The remaining 109 DETs showed inconsistent expression patterns depending on genotypes. Among them, VINC was up-regulated in KE1 but down-regulated in KB11. Previous studies have shown that D. magna exposed to cyclophosphamide, one of the most widely used drugs in cancer chemotherapy, showed upregulation of vinculin, a cytoskeletal protein involved in the regulation of focal adhesions and embryonic development [49, 50]. There is currently limited information on the specific function of VINC in the predator stress response of Daphnia, but it may play a role in increasing growth rate and fertility in response to predator cues, mediating changes in cell morphology, adhesion, and migration. However, further research is needed to determine the specific role of VINC in Daphnia’s predator stress response. The OBST-E, which belongs to the chitin binding peritrophin-A domain, was up-regulated in KE1 but down-regulated in KB11 (Table S7). Usually, growth-related chitin and cuticle genes were up-regulated by fish kairomone exposure in Daphnia [22, 32, 36]. In addition, this protein plays a role in the cuticle formation in Drosophila [51]. It has been reported that Drosophila evolved chitin-containing structures, such as cuticle or the peritrophic membranes, that serve to protect the body from hostile environments, and expression of OBST-E is involved in this [51]. OBST-E may increase the ability to survive by evolving the cuticle structure in Daphnia, which plays a protective role in the predator stress response. In summary, as shown by the phenotypic variation observations, genotype KE1 increased “Somatic Growth Rate” and “Body Length” under predator stress, while genotype KB11 decreased, which may be related to VINC and OBST-E protein. Moreover, 44 transcripts were down-regulated in KE1, but up-regulated in KB11 (Table S7). For example, FLIL, a conserved cytoplasmic protein found in a wide range of organisms, including fruit flies [52, 53], was down-regulated in KE1 but up-regulated in KB11. FLIL is known to play a role in cytoskeletal dynamics and actin remodeling, which are important processes in cellular morphology and motility [54]. However, given the limited research on FLIL expression in Daphnia and Drosophila under predation stress, it is difficult to explain the association of these DETs with the measured phenotypic variation. Further research is needed to fully understand the specific mechanisms regulating FLIL expression levels in response to predatory-induced signaling in Daphnia.

The two D. galeata genotypes found in the same environment exhibit a larger number of genotype-specific DETs compared to the common ones. Based on the six life-history traits measured, “Body Length” is the only trait showing the opposite phenotype in both genotypes (Fig. 1F). Growth-related muscle protein synthesis, which may affect body length, plays an important role in various metabolic processes, growth, and reproductive activities in crustaceans. The significance of muscle-specific genes and proteins in crustaceans has been highlighted by several research on invertebrates [55, 56]. In this study, up-regulated MMP2 and DOP2R, related to growth, were only found in KE1. MMP2 plays a role in extracellular matrix remodeling and is involved in several physiological processes including growth and development [57, 58]. MMP2 expression is up-regulated in D. pulex under predator stress, indicating its potential involvement in predator-induced phenotypic plasticity [59]. Overall, MMP2 appears to be an important mediator of predator-induced phenotypic plasticity in Daphnia, and may play a role in the regulation of growth and development under predator stress. DOP2R is a neurotransmitter known to be involved in various physiological processes, including growth and development [60]. Under predator stress, dopamine can have a role in modulating the response to stress and the allocation of resources toward to growth in Daphnia. A previous study has reported that dopamine application results in larger individuals at sexual maturity, where dopamine may affect cell proliferation and/or increase cell growth [60]. In contrast, two DETs belonging to the Lim domain, such as PXN, and MLP84B, were down-regulated in KB11, but not found in KE1 (Table S6). Lim domain proteins play crucial roles in various BP, including cytoskeleton organization, cell fate determination, and organ development [61]. Specifically, MLP84B, a cytoskeletal protein primarily expressed in Drosophila muscles, has been implicated in the regulation of actin filament organization and stability in muscle cells, which can affect muscle growth and function [62, 63]. This suggests that genotype KB11 represents an adaptive strategy to reduce body size to avoid being detected and eaten by predators by inhibiting growth by down-regulating Lim domain proteins involved in growth. However, further studies are needed to fully comprehend the specific mechanisms through which the described genotype-specific DETs are linked with the observed phenotypic variation in each genotype.

There are reports suggesting that signaling pathways regulate developmental changes during exposure to predators, which in turn affect reproduction and growth. The tachykinin (TK) receptor signaling pathway was only identified in the GO enrichment analysis of genotype KE1 (Table S8), which matures early, and has an increased body size. This pathway is involved in various physiological processes such as neurotransmission, muscle contraction, regulation of feeding behavior, and immune response in various organisms including Drosophila [64,65,66]. Especially, the TK in Drosophila brain is one of several factors that regulate the insulin-producing cells, which play important hormonal roles in the regulation of metabolic carbohydrates and lipids, as well as reproduction, growth, stress resistance, and aging [64,65,66]. Activation of the TK receptor in various environmental stress responses may lead to the release of neuropeptides, which modulate various physiological and behavioral responses such as increased heart rate, changes in locomotor activity, and altered feeding behavior. In addition, activation of the TK receptor signaling pathway has been shown to increase growth rates in Drosophila [64,65,66]. These responses are thought to be adaptive, allowing the flies to quickly respond to danger and increase their chances of survival. Overall, the TK receptor signaling pathway appears to be involved in promoting growth, and since the TK receptor signaling pathway was only activated in the genotype KE1, which showed increased growth rates under fish kairomone exposure, this is likely to have influenced growth. On the other hand, the cell differentiation, proliferation, and apoptosis-related pathways (hippo, FoxO, ECM-receptor interaction, and MAPK signaling pathways) and related to neuronal development such as Wnt signaling pathways were inhibited under fish kairomone exposure in only genotype KB11, which mature late and have decreased body size (Table S8). All these signaling pathways are involved in the regulation of embryonic development [67,68,69]. Hence, the disrupted regulation of these pathways could have detrimental impacts on embryonic development, potentially linking it to the impaired growth of genotype KB11 when exposed to fish kairomones. Among these pathways, a total of 28 DETs were mapped to the Wnt pathway, which was identified only in KB11 (Table S9). The Wnt pathway is a key signaling pathway that affects the late developmental stage of Daphnia; transcriptional profiling of D. mitsukuri exposed to fish kairomones revealed that it responds to predation risk by regulating the activity of the Wnt signaling pathway [70]. Thus, down-regulation of DETs within the Wnt pathway, which plays a pivotal role in animal development and growth [71, 72], may contribute to the smaller body size observed in genotype KB11. Further research is needed to fully elucidate the precise mechanisms and interactions within the Wnt pathway that contribute to body size determination in Daphnia.

While this study describes common and genotype-specific DETs in response to fish kairomone exposure, it is likely that DETs will change with additional populations and genotypes. Thus, to elucidate the association between the DETs of both genotypes and the observed phenotypic variations, it is crucial to generate RNA-seq data for additional D. galeata genotypes from the same and other populations. This approach would offer insights into shared or divergent life-history traits, allowing for a comprehensive understanding of the genetic basis underlying these phenotypic variations. Our findings could also be explained by the influence of epigenetic modifications, such as cytosine methylation [73]. In addition, in future studies, quantitative RT-PCR (qRT-PCR) should be performed on differentially expressed genes to find candidate genes and quantify the expression. Epigenetic factors can have a significant impact on phenotypic plasticity by modulating gene expression and influencing the response of an organism to environmental cues. In addition, intraspecific variation can have a significant impact on alternative splicing, and the diversity and abundance of isoforms can explain the intraspecific variation found under predation risk [74]. Alternative splicing is a key regulatory mechanism in eukaryotes that plays a pivotal role in transcriptional diversity and activity, significantly contributing to gene expression regulation. Long-read sequencing can serve to understand how different isoforms contribute to phenotypic variation in response to fish kairomone exposure. Our future research aims to elucidate the precise adaptive mechanisms underlying predator-induced responses.

Conclusion

This study aimed to compare the predator-induced phenotypic variation and transcriptional profiles of in two contrasting D. galeata genotypes found in the Han River, Korea. The two D. galeata genotypes investigated in the study exhibited similarities in terms of life-history traits associated with reproduction; however, in terms of growth, one genotype increased body size, whereas the other genotype decreased it under fish kairomone exposure, indicating divergent adaptive strategies. To better understand the transcripts involved in phenotypic plasticity under predation stress, we conducted RNA-seq on their transcriptomes. Our results showed that the consistently up-regulated DETs in both genotypes encoded proteins involved in digestion and reproduction, which can be interpreted as life-history traits related to reproduction. Furthermore, among the DETs exhibiting inconsistent expression patterns in both genotypes, genes associated with growth were up-regulated in KE1 but down-regulated in KB11, suggesting a potential link with life-history traits associated with growth. The present findings are in line with the results obtained from GO, KEGG, and network analyses, which demonstrated the involvement of various processes associated with reproduction and growth. Nevertheless, additional research is required to gain a comprehensive understanding of the precise mechanisms through which the described DETs result in the observed phenotypic variation in each genotype. Our analysis allows for a better understanding of the adaptation mechanisms related to reproduction and growth of two Korean D. galeata genotypes under predation stress.

Methods

Sample collection, culture, and molecular identification

Two D. galeata genotype samples were collected in the same environment (37°30’46.0” N, 126°59’53.7” E) from the Han River, Korea, on June 2020, using a 55 μm mesh sized Kitahara Zooplankton net (Samjee-tech, Korea). After transferring the collected samples to the laboratory, each egg-bearing female individual was transferred to a separate beaker for parthenogenetic reproduction. D. galeata individuals were cultured under laboratory conditions (20℃, 16 h light/8 h dark cycle, ISO medium); 1.0 mg C L− 1Chlorella vulgaris was used as food source once daily. The two genotypes examined were identified by mitochondrial cytochrome oxidase I (cox1) and NADH dehydrogenase subunit 2 (nd2) gene sequence analyses [28, 30]. As a result, we identified only two different genotypes with sequence differences in the cox1 and nd2 genes. In addition, the haplotype network based on nd2 gene revealed genotype KE1 belonged to clade D and KB11 to clade B1 (Figure S3).

Experimental design and life-history experiment

This study was conducted using two D. galeata genotypes (KE1 and KB11) cultured for three years prior to this experiment and had established clonal lines. To mimic fish predation, D. galeata was exposed to a kairomone-enriched medium produced by growing approximately 20 freshwater mandarin fish (Siniperca scherzeri) in 100 L of water. This medium was obtained from the aquaria where freshwater mandarin fish were raised. The medium in the fish aquaria was changed three times a week, and the medium containing kairomones was prepared by filtering it through a 0.45-µm-pore-sized, cellulose acetate membrane (Hyundai-Micro, Korea) and mixing it the ISO medium at a 1:2 ratio. The ratio of fish kairomone medium was set as the ratio at which the changes occur in life-history traits while ensuring the survival rate of D. galeata after several trials. The experimental design and procedures based on the previous research [32]. To minimize inter-individual variances, both genotypes were bred for two consecutive generations in a kairomone-free medium (control group) and fish kairomone medium (experimental group) prior to the experiment. After culturing until the second generation, 20 egg-bearing females of each genotype from the medium were considered as grandmothers (F0) to the experimental animals (F2). According to Tams et al., D. galeata showed strong changes in life-history traits after three generations of fish kairomone exposure [9]. To investigate the long-term effects of fish kairomones on transcript expression levels, the experiment lasted 14 days for each experimental individual. During this time, life-history traits related to reproduction such as “Age at First Reproduction,” “Number of Offspring First Brood,” “Total Number of Broods,” and “Total Number of Offspring” were recorded. In addition, life history variables related to growth such as “Somatic Growth Rate” and “Body Length” were estimated. The “Age at First Reproduction” refers to the day when the first neonates were released from the brood pouch. The “Total Number of Offspring” was calculated by the number of summing offspring produced by each individual during the induction experiment. Additionally, the “Somatic Growth Rate” was measured by subtracting the body length of the neonate at the start of the experiment (t0) from that of the adult at the end of the experiment (t14) and dividing it by the total duration of the experiment (14 days). The measurement of body length excluded the spine itself and was taken from the top of the head, passing through the middle of the eye, ending at the ventral basis of the spine [9]. The individuals were observed using CX22LED microscope (Olympus, Japan), and photographed using eXcope T500 microscope camera (DIXI Science, Korea). The statistical analyses for life-history traits were conducted using the aov() function in R version 4.2.2 [75, 76]. The egg-bearing females (F2) were pooled (n = 20) after 14 days; three biological replicates per experimental condition (control and fish kairomone exposed) were used for the two genotypes, resulting in a total of 240 individuals (two genotypes × two conditions × 20 individuals × three biological replicates; Figure S4). The samples were preserved in Qiagen RNAprotect tissue reagent (Qiagen, Valencia, CA, USA) at 4ºC overnight and subsequently stored at -80ºC until RNA extraction.

Data collection and analysis

RNA preparation and high-throughput sequencing

Since we could not obtain a sufficient amount of RNA from one individual, we pooled the experimental individuals. As embryos contain high amounts of RNA and can perceive predator cues, we accounted for the stage of egg development within the brood pouch [77]. Because sampling females in the inter-molt stage can ensure stable gene expression, only egg-bearing experimental females were pooled [77].

Total RNA was isolated from pools of 20 egg-bearing adults using the Qiagen RNeasy Plus Universal Mini Kit (Qiagen, Valencia, CA, USA) following manufacturer’s instructions after homogenizing with a disposable pestle and homogenizer. Samples were stored at − 80ºC until RNA isolation. RNA quality was assessed using an Agilent 2100 bioanalyzer (Agilent Technologies, Amstelveen, The Netherlands), and RNA quantification was conducted using an ND-2000 Spectrophotometer (Thermo Inc., DE, USA) for 12 samples (two genotypes × two conditions × three biological replicates). Libraries were constructed from the total RNA using the NEBNext Ultra II Directional RNA-seq Kit (New England BioLabs, Inc., UK). The Poly(A) RNA Selection Kit was employed to isolate mRNA (Lexogen, Inc., Austria). The isolated mRNAs were utilized for cDNA synthesis and shearing following manufacturer’s instructions. Indexing was conducted utilizing Illumina indexes 1–12, followed by enrichment by PCR. Subsequently, the libraries were assessed using the TapeStation HS D1000 Screen Tape (Agilent Technologies, Amstelveen, The Netherlands) to determine the mean fragment size. Library quantification was conducted using a library quantification kit on a StepOne Real-Time PCR System (Life Technologies, Inc., USA). High-throughput sequencing was conducted using the NovaSeq 6000 (Illumina, Inc., USA) platform, generating paired-end reads of 101 base pairs from each end.

RNA-seq quality control andde novotranscriptome assembly.

A quality control of raw sequencing data was performed using FastQC v.0.11.9 [78]. Adapter and low-quality reads (< Q20) were removed using Trimmomatic v.0.36 with the following parameters: ILLUMINACLIP: TruSeq3-PE.fa:2:30:10 TRAILING: 20 SLIDINGWINDOW: 4:15 MINLEN: 36 [79]. After trimming, the read quality was checked again with FastQC for successful adapter removal. Since the annotation of the reference genome of D. galeata is incomplete and insufficient for further analysis, a de novo approach was employed to assemble the RNA-seq reads. In addition, a de novo assembly was performed to use the transcripts found in both genotypes after fish kairomone exposure. The de novo transcriptome assembly was performed by combining the resulting quality-filtered reads of RNA-seq data of two genotypes using Trinity v.2.14.0 default parameters [80].

Assembly assessment and unigene calculation

To calculate the assembly statistics, the TrinityStats.pl script from the Trinity pipeline was utilized. The transcriptome completeness of all datasets was assessed using BUSCO v5.3.2 against the Arthropoda dataset (Arthropoda Odb10), composed of 1013 single-copy ortholog genes with default parameters [81]. To eliminate sequence redundancy, unigenes with 95% similarity were calculated using CD-Hit [82]. The reads were mapped back to the unigenes using RNA-seq by Expectation Maximization [83].

Functional annotation of unigenes

The TransDecoder pipeline with default parameters was employed to predict potential coding regions and ORFs [84]. All ORFs analyzed in this study were required to be ≥ 100 amino acids in length. BLASTp analysis was conducted on the unigenes using an e-value cutoff of 1e-5 against the UniProt database [85, 86]. The top hit for each query sequence was selected for transcriptome annotation and further characterization. HMMER v.3.1b2 was utilized against the Pfam database to identify protein domains [87, 88]. The online functional annotation tool, eggNOG-mapper, was used to determine orthologous groups [89]. GO terms enrichment of the DETs was analyzed using the ShinyGO v.0.61 [90]. Enrichment of KEGG pathways of the DETs was analyzed using KEGG mapper and KOBAS v.3.0 web-based tools [91, 92]. Only enriched categories showing a corrected p-value < 0.05 were selected.

Differential transcript expression analysis

DE analysis was performed at the isoform level using the DESeq2 package for R v4.1.3 within the Trinity pipeline [75, 93]. Results were subjected to post hoc filtering to reduce the false discovery rate by applying an adjusted p-value (padj) < 0.05 [94]. Additionally, a fold change ≥ 2 (log2FC ≤ -1 or log2FC ≥ 1) was used to filter the results. To confirm the quality of the triplicated samples, PCA was performed with DEseq2 and ggplot2 [95]. The number of shared transcripts between genotypes was visualized using the web-based tool Venn (https://bioinformatics.psb.ugent.be/webtools/Venn/).

Network analysis and community detection

To predict the interactions between these DETs, the annotated UniProt IDs of the DETs were mapped using Cytoscape stringApp v.1.5.1 using Drosophila melanogaster (fruit fly) as reference organism [96]. The interaction network of the DETs was visualized and further investigated using Cytoscape v.3.7.0 [97]. Clusters with < 10 nodes were excluded from the analysis. Node degree distribution and community structure were assessed using the NetworkAnalyzer and GLay plugins of Cytoscape, respectively [98]. The modular structure of the interaction network was identified by using the GLay implementation of the fast-greedy algorithm within the clusterMaker plugin of Cytoscape [99]. This algorithm identifies clusters by iteratively removing edges from the network and then reevaluating the connected nodes [100].

Data Availability

Sequence reads generated and/or analyzed from this study are available on the Sequence Read Archive (SRA) SRR24475174 ~ SRR24475185 under BioProject PRJNA970532 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA970532).

Abbreviations

DET:

Differentially Expressed Transcript

DE:

Differential Expression

GO:

Gene Ontology

KEGG:

Kyoto Gene and Genome Encyclopedia

ORF:

Open Reading Frames

PCA:

Principal Component Analysis

BP:

Biological Process

MF:

Molecular Function

CC:

Cellular Component

coxI:

cytochrome oxidase I

nd2:

NADH dehydrogenase subunit 2

References

  1. Stearns SC. The evolutionary significance of phenotypic plasticity. Bioscience. 1989;39:436–45.

    Article  Google Scholar 

  2. Agrawal AA. Ecology – phenotypic plasticity in the interactions and evolution of species. Sci. 2001;294:321–6.

    Article  CAS  Google Scholar 

  3. Scheiner SM, Holt RD. The genetics of phenotypic plasticity. X. variation versus uncertainty. Ecol Evol. 2012;2(4):751–67.

    Article  PubMed  PubMed Central  Google Scholar 

  4. Forsman A. Effects of genotypic and phenotypic variation on establishment are important for conservation, invasion, and infection biology. Proc Natl Acad Sci. 2014;111(1):302–7.

    Article  CAS  PubMed  Google Scholar 

  5. Biswas S, Akey JM. Genomic insights into positive selection. Trends Genet. 2006;22(8):437–46.

    Article  CAS  PubMed  Google Scholar 

  6. Vanoverbeke J, De Meester L. Clonal erosion and genetic drift in cyclical parthenogens—the interplay between neutral and selective processes. J Evol Biol. 2010;23(5):997–1012.

    Article  CAS  PubMed  Google Scholar 

  7. Swillen I, Vanoverbeke J, De Meester L. Inbreeding and adaptive plasticity: an experimental analysis on predator-induced responses in the water flea Daphnia. Ecol Evol. 2015;5(13):2712–21.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Griffiths AJF, Miller JH, Suzuki DT, Lewontin RC, Gelbart WM. An introduction to genetic analysis. New York: W. H. Freeman; 2000.

    Google Scholar 

  9. Tams V, Lüneburg J, Seddar L, Detampel JP, Cordellier M. Intraspecific phenotypic variation in life history traits of Daphnia galeata populations in response to fish kairomones. PeerJ. 2018;6:e5746.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Fine PV. Ecological and evolutionary drivers of geographic variation in species diversity. Annu Rev Ecol Evol. 2015;46:369–92.

    Article  Google Scholar 

  11. Aldana M, Maturana D, Pulgar J, García-Huidobro MR. Predation and anthropogenic impact on community structure of boulder beaches. Sci Mar. 2016;80:543–51.

    Article  Google Scholar 

  12. Stibor H, Lüning J. Predator-Induced phenotypic variation in the Pattern of Growth and Reproduction in Daphnia hyalina (Crustacea: Cladocera). Funct Ecol. 1994;8:97–101.

    Article  Google Scholar 

  13. Schoeppner NM, Relyea RA. Interpreting the smells of predation: how alarm cues and kairomones induce different prey defences. Funct Ecol. 2009;23:1114–21.

    Article  Google Scholar 

  14. Lima SL. Nonlethal effects in the ecology of predator-prey interactions. What are the ecological effects of anti-predator decision-making? Bioscience. 1998;48:25–34.

    Article  Google Scholar 

  15. Yin MB, Laforsch C, Lohr JN, Wolinska J. Predator-induced defense makes Daphnia more vulnerable to parasites. Evolution. 2011;65:1482–8.

    Article  PubMed  Google Scholar 

  16. Herzog Q, Rabus M, Ribeiro BW, Laforsch C. Inducible defenses with a twist: Daphnia barbata abandons bilateral symmetry in response to an ancient predator. PLoS ONE. 2016;11:e0148556.

    Article  PubMed  PubMed Central  Google Scholar 

  17. Lampert W. Daphnia: development of a model organism in ecology and evolution. International Ecology Institute; 2011.

  18. Miner BE, de Meester L, Pfrender ME, Lampert W, Hairston NG. Linking genes to communities and ecosystems: Daphnia as an ecogenomic model. Proc Royal Soc. 2012;279:1873–82.

    Google Scholar 

  19. Shapiro J. The importance of trophic-level interactions to the abundance and species composition of algae in lakes. Hyper Eco. 1980;2:105–16.

    CAS  Google Scholar 

  20. Boersma M, Spaak P, De Meester L. Predator-mediated plasticity in morphology, life history, and behavior of Daphnia: the uncoupling of responses. Am Nat. 1998;152(2):237–48.

    Article  CAS  PubMed  Google Scholar 

  21. Jansen M, Coors A, Stoks R, De Meester L. Evolutionary ecotoxicology of pesticide resistance: a case study in Daphnia. Ecotoxicol. 2011;20(3):543–51.

    Article  CAS  Google Scholar 

  22. An H, Do TD, Jung G, Karagozlu MZ, Kim CB. Comparative transcriptome analysis for understanding Predator-Induced Polyphenism in the water flea Daphnia pulex. Int J Mol Sci. 2018;19:2110.

    Article  PubMed  PubMed Central  Google Scholar 

  23. De Meester L, Gómez A, Okamura B, Schwenk K. The monopolization hypothesis in dispersal-gene flow paradox in aquatic organisms. Acta Oncol. 2002;23(3):121–35.

    Google Scholar 

  24. Ravindran SP, Herrmann M, Cordellier M. Contrasting patterns of divergence at the regulatory and sequence level in european Daphnia galeata natural populations. Ecol Evol. 2019;9(5):2487–504.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Castro BB, Consciência S, Gonçalves F. Life history responses of Daphnia longispina to mosquitofish (Gambusia holbrooki) and pumpkinseed (Lepomis gibbosus) kairomones. Hydrobiologia. 2007;594(1):165–74.

    Article  Google Scholar 

  26. Beckerman AP, Rodgers GM, Dennis SR. The reaction norm of size and age at maturity under multiple predator risk. J Anim Ecol. 2010;79(5):1069–76.

    Article  PubMed  Google Scholar 

  27. Brooks JL. The systematics of north american Daphnia. Mem Conn Acad Arts Sci. 1957;13:1–180.

    Google Scholar 

  28. Choi TJ, Do TD, Kim JI, Kim CB. Analysis of the complete mitogenome of Daphnia galeata from the Han River, South Korea: structure comparison and control region evolution. Funct Integr Genom. 2023;23(1):65.

    Article  CAS  Google Scholar 

  29. Bernatowicz P, Dawidowicz P, Pijanowska J. Phenotypic plasticity and developmental noise in hybrid and parental clones of Daphnia longispina complex. Aquat Ecol. 2021;55:1179–88.

    Article  Google Scholar 

  30. Ishida S, Taylor DJ. Quaternary diversification in a sexual holarctic zooplankter, Daphnia galeata. Mol Ecol. 2007;16:569–82.

    Article  PubMed  Google Scholar 

  31. Stich HB, Lampert W. Predator evasion as an explanation of diurnal Vertical Migration by Zooplankton. Nature. 1981;293:396–8.

    Article  Google Scholar 

  32. Tams V, Nickel JH, Ehring A, Cordellier M. Insights into the genetic basis of predator-induced response in Daphnia galeata. Ecol Evol. 2020;10(23):13095–108.

    Article  PubMed  PubMed Central  Google Scholar 

  33. Carter MJ, Silva-Flores P, Oyanedel JP, Ramos-Jiliberto R. Morphological and life-history shifts of the exotic cladoceran Daphnia exilis in response to predation risk and food availability. Limnologica. 2013;43(3):203–9.

    Article  Google Scholar 

  34. Lampert W. Phenotypic plasticity of the size at first reproduction in Daphnia: the importance of maternal size. Ecology. 1993;74(5):1455–66.

    Article  Google Scholar 

  35. Gliwicz ZM, Boavida MJ. Clutch size and body size at first reproduction in Daphnia pulicaria at different levels of food and predation. J Plankton Res. 1996;18(6):863–80.

    Article  Google Scholar 

  36. Rozenberg A, Parida M, Leese F, Weiss LC, Tollrian R, Manak JR. Transcriptional profiling of predator-induced phenotypic plasticity in Daphnia pulex. Front Zool. 2015;12(1):18.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Stibor H. The role of yolk protein dynamics and predator kairomones for the life history of Daphnia magna. Ecology. 2002;83:362–9.

    Article  Google Scholar 

  38. Hales NR, Schield DR, Andrew AL, Card DC, Walsh MR, Castoe TA. Contrasting gene expression programs correspond with predator-induced phenotypic plasticity within and across generations in Daphnia. Mol Ecol. 2017;26:5003–15.

    Article  CAS  PubMed  Google Scholar 

  39. Malik A, Lee EJ, Jan AT, Ahmad S, Cho KH, Kim J, Choi I. Network Analysis for the identification of differentially expressed hub genes using myogenin knock-down muscle Satellite cells. PLoS ONE. 2015;10:e0133597.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Malik A, Lee J, Lee J. Community-based network study of protein-carbohydrate interactions in plant lectins using glycan array data. PLoS ONE. 2014;9:e95480.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Riessen HP. Predator-induced life history shifts in Daphnia: a synthesis of studies using meta-analysis. Can J Fish Aquat Sci. 1999;56:2487–94.

    Article  Google Scholar 

  42. Lynch M. The evolution of cladoceran life hitories. Q Rev Biol. 1980;55(1):23–42.

    Article  Google Scholar 

  43. Garay-Narvaez L, Ramos-Jiliberto R. Induced defenses within food webs: the role of community trade-offs, delayed responses, and defense specificity. Ecol Complex. 2009;6(3):383–91.

    Article  Google Scholar 

  44. Tollrian R, Harvell CD. The ecology and evolution of inducible defenses. Princeton university press; 1999.

  45. Lass S, Spaak P. Hemically induced anti-predator defences in plankton: a review. Hydrobiologia. 2003;491(1–3):221–39.

    Article  Google Scholar 

  46. Storz P. Forkhead homeobox type O transcription factors in the responses to oxidative stress. Antioxid Redox Signal. 2011;14(4):593–605.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Schwarzenberger A, Küster CJ, von Elert E. Molecular mechanisms of tolerance to cyanobacterial protease inhibitors revealed by clonal differences in Daphnia magna. Mol Ecol. 2012;21:4898–911.

    Article  CAS  PubMed  Google Scholar 

  48. Schwarzenberger A, Zitt A, Kroth P, Mueller S, Von Elert E. Gene expression and activity of digestive proteases in Daphnia: effects of cyanobacterial protease inhibitors. BMC Physiol. 2010;10:6.

    Article  PubMed  PubMed Central  Google Scholar 

  49. Laneve P, Gioia U, Ragno R, Altieri F, Franco CD, Santini T, Arceci M, Bozzoni I, Caffarelli E. The tumor marker human placental protein 11 is an endoribonuclease. J Biol Chem. 2008;283:34712–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Dumbauld DW, García AJ. A helping hand: how vinculin contributes to cell-matrix and cell-cell force transfer. Cell Adhes Migr. 2014;8:550–7.

    Article  Google Scholar 

  51. Behr M, Hoch M. Identification of the novel evolutionary conserved obstructor multigene family in invertebrates. FEBS Lett. 2005;579:6827–33.

    Article  CAS  PubMed  Google Scholar 

  52. Campbell HD, Schimansky T, Claudianos C, Ozsarac N, Kasprzak AB, Cotsell JN, Young IG, de Couet HG, Miklos GL. The Drosophila melanogaster flightless-I gene involved in gastrulation and muscle degeneration encodes gelsolin-like and leucine-rich repeat domains and is conserved in Caenorhabditis elegans and humans. Proc Natl Acad Sci. 1993;90:11386–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. de Couet HG, Fong KS, Weeds AG, Mclaughlin PJ, Miklos GL. Molecular and mutational analysis of a gelsolin-family member encoded by the flightless I gene of Drosophila melanogaster. Genetics. 1995;141:1049–59.

    Article  PubMed  Google Scholar 

  54. Strudwick XL, Cowin AJ. Multifunctional roles of the actin-binding protein flightless I in inflammation, cancer and wound healing. Front Cell Dev Biol. 2020;8(1394).

  55. Jung H, Lyons RE, Dinh H, Hurwood DA, McWilliam S, Mather PB. Transcriptomics of a giant freshwater prawn (Macrobrachium rosenbergii): De Novo Assembly, annotation and marker Discovery. PLoS ONE. 2011;6(12):e27938.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Lv J, Liu P, Gao B, Wang Y, Wang Z, Chen P, Li J. Transcriptome analysis of the Portunus trituberculatus: de novo assembly, growth-related gene identification and marker discovery. PLoS ONE. 2014;9(4):e94055.

    Article  PubMed  PubMed Central  Google Scholar 

  57. Page-McCaw A, Serano J, Santé JM, Rubin GM. Drosophila matrix metalloproteinases are required for tissue remodeling, but not embryonic development. Dev Cell. 2003;11:95–106.

    Article  Google Scholar 

  58. Page-McCaw A, Ewald AJ, Werb Z. Matrix metalloproteinases and the regulation of tissue remodelling. Nat Rev Mol. 2007;11:221–33.

    Article  Google Scholar 

  59. Spanier KI, Leese F, Mayer C, Colbourne JK, Gilbert D, Pfrender ME, Tollrian R. Predator-induced defences in Daphnia pulex: selection and evaluation of internal reference genes for gene expression studies with realtime PCR. BMC Mol Biol. 2010;11:50.

    Article  PubMed  PubMed Central  Google Scholar 

  60. Weiss LC, Leese F, Laforsch C, Tollrian R. Dopamine is a key regulator in the signalling pathway underlying predator-induced defences in Daphnia. Proc R Soc B: Biol Sci. 2015;282:20151440.

    Article  Google Scholar 

  61. Bach I. The LIM domain: regulation by association. Mech Dev. 2002;91:5–17.

    Article  Google Scholar 

  62. Mery A, Taghli-Lamallem O, Clark KA, Beckerle MC, Wu XS, Ocorr K, Bodmer R. The Drosophila muscle LIM protein, Mlp84B, is essential for cardiac function. J Exp Biol. 2008;211:15–23.

    Article  CAS  PubMed  Google Scholar 

  63. Stronach BE, Renfranz PJ, Lilly B, Beckerle MC. Muscle LIM proteins are Associated with muscle sarcomeres and require dMEF2 for their expression during Drosophila myogenesis. Mol Biol Cell. 1999;10:2329–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Birse RT, Söderberg JA, Luo J, Winther ÅM, Nässel DR. Regulation of insulin-producing cells in the adult Drosophila brain via the tachykinin peptide receptor DTKR. J Exp Biol. 2011;214:4201–8.

    Article  CAS  PubMed  Google Scholar 

  65. Nässel DR, Zandawala M, Kawada T, Satake H. Tachykinins: neuropeptides that are ancient, diverse, widespread and functionally pleiotropic. Front Neurosci. 2019;13:1262.

    Article  PubMed  PubMed Central  Google Scholar 

  66. Nässel DR, Broeck JV. Insulin/IGF signaling in Drosophila and other insects: factors that regulate production, release and post-release action of the insulin-like peptides. Cell Mol Life Sci. 2016;73:271–90.

    Article  PubMed  Google Scholar 

  67. Ma S, Meng Z, Chen R, Guan KL. The hippo pathway: Biology and pathophysiology. Annu Rev Biochem. 2019;88:577–604.

    Article  CAS  PubMed  Google Scholar 

  68. Martins R, Lithgow GJ, Link W. Long live FOXO: unraveling the role of FOXO proteins in aging and longevity. Aging Cell. 2016;15(2):196–207.

    Article  CAS  PubMed  Google Scholar 

  69. Zhang S, Zhou X, Zhang C, Zhao C, Li W, Wang D, Xu S. Expression of the senescence-related gene FoxO in Daphnia pulex and its role in the regulation of reproductive transformation (Branchiopoda, Cladocera). Crustaceana. 2022;95(8–9):961–83.

    Article  Google Scholar 

  70. Zhang X, Blair D, Wolinska J, Ma X, Yang W, Hu W, Yin M. Genomic regions associated with adaptation to predation in Daphnia often include members of expanded gene families. Proc R Soc B: Biol Sci. 2021;288.

  71. Goto S, Hayashi S. Proximal to distal cell communication in the Drosophila leg provides a basis for an intercalary mechanism of limb patterning. Dev. 1999;126:3407–13.

    Article  CAS  Google Scholar 

  72. Holstein TW. The evolution of the wnt pathway. Cold Spring Harb Perspect biol. 2012;4:e007922.

    Article  Google Scholar 

  73. Asselman J, De Coninck DIM, Vandegehuchte MB, Jansen M, Decaestecker E, De Meester L, De Schamphelaere KAC. Global cytosine methylation in Daphnia magna depends on genotype, environment, and their interaction. Environ Toxicol Chem. 2015;34(5):1056–61.

    Article  CAS  PubMed  Google Scholar 

  74. Suresh S, Crease TJ, Cristescu ME, Chain FJJ. Alternative splicing is highly variable among Daphnia pulex lineages in response to acute copper exposure. BMC Genom. 2020;21:433.

    Article  CAS  Google Scholar 

  75. R Core Team. R: A Language and Environment for Statistical Computing. 2017. https://www.R-project.org/.

  76. MacFarland TW, Yates JM. Oneway Analysis of Variance (ANOVA). Using R for Biostatistics. Springer, Cham. 2021.

    Chapter  Google Scholar 

  77. Weiss LC, Heilgenberg E, Deussen L, Becker SM, Kruppert S, Tollrian R. Onset of kairomone sensitivity and the development of inducible morphological defenses in Daphnia pulex. Hydrobiologia. 2016;779:135–45.

    Article  CAS  Google Scholar 

  78. Andrews S, FastQC. A quality control tool for high throughput sequence data. 2010. Available online at: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/.

  79. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  80. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, Chen Z, Mauceli E, Hacohen N, Gnirke A, Rhind N, di Palma F, Birren BW, Nusbaum C, Lindblad-Toh K, Friedman N, Regev A. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29:644–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  81. Seppey M, Manni M, Zdobnov EM. BUSCO: assessing Genome Assembly and Annotation Completeness. Methods mol biol. 2019;1962:227–45.

    Article  CAS  PubMed  Google Scholar 

  82. Li W, Godzik A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics. 2006;22:1658–9.

    Article  CAS  PubMed  Google Scholar 

  83. Li B, Dewey CN. RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 2011;12:e323.

    Article  Google Scholar 

  84. Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, Couger MB, Eccles D, Li B, Lieber M, MacManes MD, Ott M, Orvis J, Pochet N, Strozzi F, Weeks N, Westerman R, William T, Dewey CN, Henschel R, LeDuc RD, Friedman N, Regev A. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protoc. 2013;8:1494–512.

    Article  CAS  PubMed  Google Scholar 

  85. Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, Madden TL. BLAST+: Architecture and applications. BMC Bioinform. 2009;10:421.

    Article  Google Scholar 

  86. The UniProt Consortium. UniProt: a worldwide hub of protein knowledge. Nucleic Acids Res. 2019;47:D506–15.

    Article  Google Scholar 

  87. Finn RD, Clements J, Arndt W, Miller BL, Wheeler TJ, Schreiber F, Bateman A, Eddy SR. HMMer web server: 2015 update. Nucleic Acids Res. 2015;43:W30–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  88. El-Gebali S, Mistry J, Bateman A, Eddy SR, Luciani A, Potter SC, Qureshi M, Richardson LJ, Salazar GA, Smart A, Sonnhammer ELL, Hirsh L, Paladin L, Piovesan D, Tosatto SCE, Finn RD. The pfam protein families database in 2019. Nucleic Acids Res. 2019;47(D1):D427–32.

    Article  CAS  PubMed  Google Scholar 

  89. Huerta-Cepas J, Forslund K, Coelho LP, Szklarczyk D, Jensen LJ, von Mering C, Bork P. Fast genome-wide functional annotation through orthology assignment by eggNOG-Mapper. Mol Biol Evol. 2017;34(8):2115–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  90. Ge SX, Jung D, Yao R. ShinyGO: a graphical gene-set enrichment tool for animals and plants. Bioinformatics. 2020;36:2628–9.

    Article  CAS  PubMed  Google Scholar 

  91. Kanehisa M, Sato Y. KEGG Mapper for inferring cellular functions from protein sequences. Protein Sci. 2020;29:28–35.

    Article  CAS  PubMed  Google Scholar 

  92. Bu D, Luo H, Huo P, Wang Z, Zhang S, He Z, Wu Y, Zhao L, Liu J, Guo J, Fan S, Cao W, Yi L, Zhao Y, Kong L. KOBAS-i: intelligent prioritization and exploratory visualization of biological functions for gene enrichment analysis. Nucleic Acids Res. 2021;49:W317–25.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  93. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  PubMed  PubMed Central  Google Scholar 

  94. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Series B Stat Methodol. 1995;57:289–300.

    Google Scholar 

  95. Wickham H. ggplot2: elegant graphics for data analysis. Springer; 2009.

  96. Doncheva NT, Morris JH, Gorodkin J, Jensen LJ. Cytoscape StringApp: network analysis and visualization of Proteomics Data. J Proteome Res. 2019;18:623–32.

    Article  CAS  PubMed  Google Scholar 

  97. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  98. Su G, Kuchinsky A, Morris JH, States DJ, Meng F. GLay: Community structure analysis of biological networks. Bioinformatics. 2010;26:3135–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  99. Morris JH, Apeltsin L, Newman AM, Baumbach J, Wittkop T, Su G, Bader GD. Ferrin TE.clusterMaker: a multi-algorithm clustering plugin for Cytoscape. BMC Bioinform. 2011;12:436.

    Article  CAS  Google Scholar 

  100. Koh GC, Porras P, Aranda B, Hermjakob H, Orchard SE. Analyzing protein-protein interaction networks. J Proteome Res. 2012;11:2014–31.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We acknowledged Jong-Won Baek and Min-Ho Mun (Sangmyung University, Korea) sample collection for the assistance of sample collection in field.

Funding

This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) [No. NRF- 2022R1A2C3002750].

Author information

Authors and Affiliations

Authors

Contributions

TJC and CBK contributed to the conception and design of the research. TJC, SMH, and AM contributed to drawing the figures, tables, and analysing the data. TJC, SMH, and AM contributed to the writing and drafting of the manuscript. AM and CBK contributed to the critical revision of the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Chang-Bae Kim.

Ethics declarations

Competing interests

The authors declare no competing interests.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic Supplementary Material

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Choi, TJ., Han, SM., Malik, A. et al. Comparative transcriptome analysis of two Daphnia galeata genotypes displaying contrasting phenotypic variation induced by fish kairomones in the same environment of the Han River, Korea. BMC Genomics 24, 580 (2023). https://doi.org/10.1186/s12864-023-09701-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-023-09701-x

Keywords