Inducible responses under joint predation risks
Antagonistic responses would put the prey into a dilemma regarding which predator to respond to (Fig. S2). In our experiments, size at 16th day was not significantly expressed under joint predation risks (CH + F vs. C: P = 0.764). Consistently, unique responses to Chaoboruspredation risk were significantly affected by the presence of fish predation risk, that is, horn expression in the combination treatment was greatly impaired compared with that under the Chaoboruspredation risk (maturity: P = 0.012; 16th day: P< 0.001). This effect was highly correlated to individual development, that is, in the combination treatment, C. cornutahorns were significantly expressed at maturity (CH + F vs. C: P = 0.007) but not at 16th day. Conversely, the effects of joint predation risks were observed in neither unique responses to fish nor general responses. As such, the responses of size at maturity (CH + F vs. C:P < 0.001), time to first brood (CH + F vs. C: P= 0.027), and neonate size (CH + F vs. C: P = 0.003) were significantly altered in the combination treatment but were not significantly different from those under the fish predation risk treatment. Thus, various predator-dependent unique responses could be simultaneously expressed under antagonistic predation risks, although the prey could not completely avoid the effect of other predation risks, and the co-expression of unique responses was related to individual development. These results basically support our hypotheses.

Overview of the assembled transcriptome

We obtained respectively 78,341,484, 72,493,072, and 69,930,191 clean reads for C, CH, and F. A total of 37,120 unigenes were assembled using Trinity, and each sample contained > 66.54% assembled unigenes (Table S2). Through comparisons against the Nr, KEGG, COG, and SWISS-PROT databases, 18,343 unigenes (49.4%) were annotated (Fig. S3). Regarding species distribution in the Nr database, C. cornutashowed the highest comparison rate with D. magna (26.33%), followed by D. pulex (2.6%) (Fig. S4A). The unigenes enriched in the KOG database were classified into transcription, ribosomal structure, gene replication, recombination, and repair (Fig. S4B). The annotated GO terms were mainly associated with metabolic processes, cellular processes, cell parts, and binding (Fig. S4C).

DEGs

Paired samples within the same treatment showed high Pearson’s correlation coefficients (≥0.90) and were clustered in principal component analysis (Fig. 2), which conformed to the requirements of biological repetition. Furthermore, the expression patterns of DEGs showed significant concordance between RT-qPCR and RNA-Seq (Fig. S5), indicating that our RNA-Seq analysis of the expression data was reliable.
Compared with the control, Chaoborus and fish predation risks significantly affected the expression of 1,515 and 846 genes, respectively (Fig. 3A). Among these, there were 1,399 unique DEGs in CH, 730 unique DEGs in F, 114 general DEGs, and two antagonistic DEGs (Fig. 3B and Table S3). Considering the DEGs triggered by different predators, we further analyzed the differences in the enriched pathways of C. cornuta (Supplementary Table S4) and identified the major DEGs and pathways related to inducible defensive traits (Table 2). The DEGs ofC. cornuta under Chaoborus predation risk, including cuticle protein, fatty acyl-CoA reductase, and trypsin genes, were mainly enriched in cutin, suberine, and wax biosynthesis; protein digestion and absorption; and steroid hormone biosynthesis. Meanwhile, the DEGs of C. cornuta under fish predation risk, including ribosomal protein, actin, and short-chain type dehydrogenase genes, were mainly enriched in cysteine and methionine metabolism, ribosome, phototransduction, and unsaturated fatty acid biosynthesis. The general DEGs, including cysteine proteinase, HSP70, actin, and alpha-tubulin genes, were enriched in apoptosis pathways and antigen processing and presentation. Therefore, specific unique responses to different predators were dominant at the transcriptional level: unique DEGs (2,129 genes) > general DEGs (114 genes) > antagonistic DEGs (2 genes).

Correlation between different inducible responses

Nine individual traits and 40 representative DEGs (Table S5) revealed strong Pearson’s correlations within unique responses to fish orChaoborus larvae (Fig. 4). For instance, the response of horn expression at maturity was significantly correlated with changes in horns at 16th day (P = 0.05) and total offspring number (P= 0.05). However, the unique responses to fish and Chaoboruslarvae were weakly correlated in both analyses. Only a few significant correlations appeared on the weak correlation area in transcriptional expression. For instance, Unigene0008297 in ribonucleoprotein component protein expression was significantly and positively correlated with Unigene0027903 (P= 0.049) and Unigene0005204 (P= 0.013) in aspartokinase and adenosylhomocysteinase expression, respectively (Fig. 4). Therefore, the unique responses triggered by different predators are weakly coupled.

Discussion

Our experiments proved that C. cornuta responds differently to antagonistic predation risks of fish and Chaoborus larvae. Based on the classification of inducible responses, our results revealed for the first time that specific predator-dependent unique responses are dominant, followed by general responses, whereas the antagonistic responses are rare. In addition, the unique responses triggered by different predators were extremely weakly coupled and could be elicited simultaneously. Overall, the experimental results support our hypothesis that predator-dependent unique responses are dominant, helping the prey avoid the dilemma of which predator to respond to. Nonetheless, the unique responses to Chaoborus larvae cannot completely avoid the effects of fish predation risk, implying the presence of complex costs and limitations of unique responses to multiple predators.
In C. cornuta , horn expression and body enlargement are adaptive inducible traits to Chaoborus larva predation (Gu et al., 2021; Riessen & Trevett-Smith, 2009). A larger individual size at late developmental stages requires rapid growth and high food intake (Gianuca, Pantel, & De Meester, 2016), thereby promoting the brood and total offspring numbers. The horns are formed by the carapace, which comprises two layers of dermal cells and is covered by chitin in combination with cuticle proteins (Charles, 2010). Therefore, the development of morphological defensive traits involves a series of changes in the expression of chitin, hormone, and epidermal formation genes at different times (Christjani, Fink, & Elert, 2016; Miyakawa et al., 2010) as well as the regulation of epidermal cell growth via endocrine hormones (Weiss, Leese, Laforsch, & Tollrian, 2015). In the present study, significant changes in the expression of genes and pathways involved in cuticle protein, cutin, suberine, and wax biosynthesis as well as steroid hormone biosynthesis were noted. The expression of these genes may promote the synthesis of related substances (Fig. 5) and regulate individual growth (Edgar, 2006). However, the growth and development of cladocerans requires continuous molting and formation of a new carapace; thus, the maintenance of horns requires continuous substance synthesis, which may result in constant distribution costs (Auld et al., 2010). Furthermore, in C. cornuta , Chaoborus predation risk altered digestion and absorption by modulating the expression of the trypsin gene, which may affect the digestion and resource allocation strategies (Von Elert et al., 2004).
Smaller size, earlier reproduction, and increased brood number are adaptive responses to fish predation risks, which are similar to the typical responses of other cladocerans under fish predation (Diel et al., 2020). In terms of gene expression, our results showed that genes encoding actin and ribosomal proteins were down-regulated under fish predation risks. Since actin plays an important role in cytoskeletal structure, its inhibition may result in a smaller size. Similar results have been reported in previous studies on inducible defense responses ofD. magna (Effertz et al., 2015; Pijanowska & Kloc, 2004). On the contrary, Schwarzenberger, Courts, and Von Elert (2009) revealed that the actin genes in D. magna were up-regulated under fish predation risks. Since gene expression is jointly regulated by transcriptional regulators and related proteins (Stibor, 2002), differential expression patterns could be observed at different time points (Effertz & Von Elert, 2014). Ribosomal proteins are responsible for protein assembly and translation; thus, the down-regulation of ribosomal protein may inhibit the synthesis of proteins essential for individual growth and development of C. cornuta (Zhou, Liao, Liao, Liao, & Lu, 2015). In enrichment analysis, some DEGs were found to be enriched in multiple pathways. Significantly enriched phototransduction altered the visual perception of Daphnia(Mahato et al., 2014), which may be an adaptation to behavioral responses, such as habitat selection (Loose & Dawidowicz, 1994) and escape behavior (Pietrzak, Pijanowska, & Dawidowicz, 2017). InDaphnia , fish predation reduced unsaturated fatty acid levels in neonates, rendering them vulnerable to starvation (Stibor & Navarra, 2000); therefore, the significantly enriched unsaturated fatty acid pathways may alter the distribution of these components. Furthermore, the longevity regulating pathway was significantly enriched under fish predation risks, which may incur an opportunity cost, that is, a decline in lifespan (Dawidowicz, Predki, & Pietrzak, 2010).
Under different predation risks, C. cornuta showed general responses, such as the expression of cysteine protease, heat shock protein, actin, and tubulin genes. The cDNA sequence of crustacean cysteine is similar to that of insect cathepsin L, which regulates the molting cycle and promotes cell death during development (Agrawal, Bagchi, & Bagchi, 2005). Thus, cysteine protease likely affected molting and increased brood number in the present study. The up-regulation of heat shock proteins is an adaptive response to various environmental stresses, including predation risks (Pijanowska & Kloc, 2004). As this response is rapid and returns to the previous state after long-term treatment (Pauwels, Stoks, & De Meester, 2005; Pauwels, Stoks, Decaestecker, & De Meester, 2007), the down-regulation of heat shock protein genes likely promoted the recovery of heat shock proteins in the present study. Similarly, general responses of the actin and tubulin genes, which are involved in cytoskeleton formation and other life activities have been observed in Daphnia (Pijanowska & Kloc, 2004), although their specific functions warrant further research (Chen et al., 2018).
Based on the classification of inducible responses, we further analyzed the early transcriptional data of Daphnia magna in response to vertebrate (sticklebacks) and invertebrate (Triops ) predation risks (Orsini et al., 2016). The results are consistent with our finding that unique responses to different predators are dominant, while antagonistic responses are rare (Fig. 6 and Table S6). Moreover, previous weighted gene co-expression network analysis revealed that the most significantly co-expressed gene networks under vertebrate and invertebrate predation risks were unique (Orsini et al., 2018). Thus, diverse predator-dependent unique responses are favored by cladocerans during their co-evolution with multiple predators. For successful evolution of diverse predator-dependent unique responses, genotype, selection, and cost are important. First, the genotypes of cladocerans in ponds or lakes are highly diverse and the inducible traits of different clones are uncoupled (Boersma, Spaak, & De Meester, 1998; Decaestecker, De Meester, & Mergeay, 2009; Stoks, Govaert, Pauwels, Jansen, & De Meester, 2016). Second, multiple predators produce variable selection effects, which contribute to predator-dependent unique defensive traits (Herzog & Laforsch, 2013; Heynen, Bunnefeld, & Borcherding, 2017). Finally, environmental costs, such as altered predator regimes, may exceed maintenance costs (Decaestecker et al., 2002; Tollrian, 1995; Yin, Laforsch, Lohr, & Wolinska, 2011). Therefore, diverse predator-dependent unique responses are favored by prey.
In the present study, prey exhibited coupled responses to the same predator but extremely weakly coupled predator-dependent unique responses to antagonistic predation risks. Evidently, prey can alter resource allocation strategies under a single type of predation risk, resulting in an array of adaptive responses (Reede, 1995). In addition, the co-expression of uncoupled unique inducible responses can help prey avoid the dilemma of which predator to respond to, thus improving their survival rate under antagonistic predation risks. For instance, smallerC. cornuta individuals are less likely to be found by fish (O’Brien, 1987). Simultaneously, horn expression renders C. cornuta less vulnerable to predation by Chaoborus larvae (Gu et al., 2021). However, the effects of joint predation risks on the expression of predator-dependent unique responses are prominent and could be influenced by development, indicating a complex trade-off underlying adaptation to multiple predation risks (Riessen & Gilbert, 2019). Therefore, further studies are warranted to elucidate the mechanisms of the phenomenon that diverse inducible responses of prey can be elicited simultaneously under joint predation risks.

Conclusions

Through responses of individual traits and transcriptome, the present study revealed the inducible responses of C. cornuta to predation by Chaoborus larvae and fish. To cope with such antagonistic predation risks, C. cornuta mainly altered cuticle gene expression and developed horns under Chaoborus predation risk, while it altered ribosome gene expression and reduced body size under fish predation risk. Our analysis of these inducible responses revealed for the first time that predator-dependent unique responses are dominant and antagonistic responses are rare, implying that predator-dependent unique responses are favored by cladocerans during their long-term co-evolution with multiple predators, thereby lowering the environmental cost of inducible defenses. However, unique responses to one predator cannot completely avoid their dilemma of responding to other predators, although this dilemma is relatively limited. The present study expands our understanding of the evolution and expression of inducible defenses under multiple predation risks.