Running title: Inducible responses to multiple predators
Lei Gu, Shanshan Qin, Yunfei Sun, Jing Huang, Siddiq Akbar, Lu Zhang, Zhou Yang*
Jiangsu Province Key Laboratory for Biodiversity and Biotechnology, School of Biological Sciences, Nanjing Normal University, 1 Wenyuan Road, Nanjing 210023, China
* Corresponding author: Zhou Yang
TEL: +86-25-85891671.
E-mail: yangzhou@njnu.edu.cn
Jiangsu Province Key Laboratory for Biodiversity and Biotechnology, School of Biological Sciences, Nanjing Normal University, 1 Wenyuan Road, Nanjing 210023, China

Abstract

Inducible defenses of prey are evolved under diverse and variable predation risks. However, during the co-evolution of prey and multiple predators, the responses of prey to antagonistic predation risks, which may put the prey into a dilemma of responding to predators, remain unclear. Based on antagonistic predation pressure from an invertebrate (Chaoborus larvae) and a vertebrate (Rhodeus ocellatus ) predator, we studied the responses of multiple traits and transcriptomes of the freshwater crustacean Ceriodaphnia cornuta under multiple predation risks. Chaoborus predation risk altered the expression of genes encoding cuticle proteins and modulated the biosynthesis of steroid hormones, cutin, suberine, and wax, leading to the development of horns and increase in size at the late developmental stage. Meanwhile, fish predation risk primarily triggered genes encoding ribosomes and those involved in unsaturated fatty acid biosynthesis and cysteine and methionine metabolism, resulting in smaller individual size and earlier reproduction. Inducible responses of both transcriptome and individual traits revealed that predator-dependent unique responses were dominant and the dilemma of antagonistic responses was relatively limited. However, the unique individual traits in response to invertebrate predation could be significantly impaired by vertebrate predation risk, even though the unique responses to different predators were extremely weakly correlated and could be elicited simultaneously. These results indicate that diverse predator-dependent unique responses are favored by Ceriodaphnia during its co-evolution with multiple predators. Nonetheless, Ceriodaphnia is not a generalist that can fully adopt all predator-dependent unique responses simultaneously under multiple predation risks.
Keywords: Cladoceran; Chaoborus ; Fish; Inducible defense; Predation risk

Introduction

In the co-evolution of predators and prey, defense is critical for prey survival. Under variable predation risks, inducible defenses triggered by predation cues are favored by prey (Tollrian & Harvell, 1999). For successful defense against predators, prey organisms adopt various inducible protective characteristics, including behavioral (De Meester, 1993), morphological (Gu et al., 2021), chemical defense (Selander et al., 2015), and life history (Kvile, Altin, Thommesen, & Titelman, 2021) traits. Since diverse predation risks can prevent the stable expression of an inducible defensive trait (Steiner & Auld, 2012), the present study sought to understand the responses of prey to multiple predation risks, particularly to predators exerting antagonistic selection pressures.
Inducible defenses are common in aquatic organisms, such as phytoplankton (Lürling, 2020), zooplankton (Diel, Kiene, Martin-Creuzburg, & Laforsch, 2020), amphibians (Mitchell, Bairos-Novak, & Ferrari, 2017), and fish (Brönmark & Miner, 1992). Through the integration of inducible defense research, we classified responses to multiple predators into two major types. The first type is the general response, which evolves through diffusion co-evolution and represents reciprocal adaptation in response to similar predators; for instance, mayflies adopt the same avoidance behavior under predation risks from different fish (Alvarez, Landeira-Dabarca, & Peckarsky, 2014). The second type is the specific response, which evolves through pairwise co-evolution between specific predators and prey, such as the immune responses to pathogens (Westra et al., 2015) and the inducible crests of Daphnia in response to Notonecta predation (Grant & Bayly, 1981). Moreover, the specific responses of different traits can be further subdivided into antagonistic responses to the same trait and unique responses to separate traits. Under antagonistic selection pressures, if unique responses to separate traits are dominant, the complex defense responses may incur maintenance cost, that is, energetic costs of the sensory and regulatory mechanisms. Moreover, if the prey primarily exhibits antagonistic responses of the same traits, inducible responses to a predator may incur an environmental cost, that is, vulnerability to other predators (Auld, Agrawal, & Relyea, 2010); this environmental cost is further linked to selection and plays a dual role in the evolution of inducible defenses (Decaestecker, De Meester, & Ebert, 2002). Therefore, we hypothesized that predator-dependent unique responses are dominant, helping, at least in part, avoid the dilemma of the prey for responding to predators under multiple predation risks.
In aquatic ecosystems, cladocerans are at the middle of the food chain, acting as a food resource for insects and fishes (Miner, De Meester, Pfrender, Lampert, & Hairston, 2012). These invertebrate and vertebrate predators constitute antagonistic selection pressures on the size or habitat selection of waterfleas; for instance, larger plankton are vulnerable to large visual predators (e.g., fish) but less vulnerable to small ambush predators (e.g., Chaoborus larvae) (Swift, 1992). Presumably, inducible responses of Daphnia to visual predators primarily include life history changes (Effertz & Von Elert, 2014), which are rather different from the defensive traits triggered by invertebrate predators, such as the development of “twist” (Herzog, Rabus, Ribeiro, & Laforsch, 2016), neck teeth (Tollrian, 1993), and horns (Gu, Qin, Zhu, et al., 2020). These reports support our hypothesis.
Antagonistic inducible traits are commonly expressed by Daphnia . For instance, Daphnia hyalina shows completely opposite responses of size and reproduction under predation pressures by vertebrates and invertebrates (Stibor & Lüning, 1994). Daphnia galeata prefers deeper habitats under fish predation risk, while inhabits upper water layers under Chaoborus predation risk (Dodson, 1988). In addition, general responses, such as the development of an elongated tail spine, are observed in Daphnia in response to fish,Triops , and Notonecta (Gu, Qin, Lu, et al., 2020; Ritschar, Rabus, & Laforsch, 2020). Consistent with inducible defensive traits, both general and specific responses appear at the molecular level; for instance, in Daphnia magna , actin and tubulin expression is decreased under Chaoborus larvae or fish predation risks (Pijanowska & Kloc, 2004), ribosomal protein and vitellogenin expression is increased under fish predation risk (Effertz, Mueller, & Von Elert, 2015), and cuticle protein expression is increased but vitellogenin expression is decreased under Triops predation risk (Otte, Fröhlich, Arnold, & Laforsch, 2014). Overall, in Daphnia , a given species exhibits various types of responses under antagonistic predation risks. However, from these sporadic studies on differentDaphnia species and clones, we cannot determine the precise type of response preferred by the prey.
Ceriodaphnia cornuta is a widely distributed species with sensitive inducible defensive traits (Gu et al., 2021; Qin et al., 2021), providing a suitable model for testing our hypothesis. Since some inducible traits are hidden (Laforsch, Ngwa, Grill, & Tollrian, 2004), research on a few traits is insufficient. In recent years, omics technologies have furthered our understanding of the mechanisms of inducible defense (Hales et al., 2017; Zhang et al., 2021). Therefore, to test our hypothesis that predator-dependent unique responses are dominant, we assessed multiple traits and transcriptomes of C. cornuta in response to predation pressures from Chaoborus larvae and fish, respectively. Additionally, to test the hypothesis that predator-dependent unique responses help prey avoid the dilemma of responding to predators, we examined the expression of individual traits under joint predation risks and explored the correlations between various inducible responses.

Materials and methods

Predation risks

Predation risks were simulated using different predator-conditioned medium prepared following the methodology described by Gu, Qin, Zhu, et al. (2020). We cultured four Rhodeus ocellatus or 100Chaoborus sp. larvae in aged tap water, fed them C. cornuta for 6 h, and then transferred them into 2 L of COMBO medium (Kilham, Kreeger, Lynn, Goulden, & Herrera, 1998) for 18 h. The predator-conditioned medium stock containing different predator kairomones (Hahn, Effertz, Bigler, & von Elert, 2019; Weiss et al., 2018) was passed through a 0.22 μm glass fiber filter (Millipore), and the filtrate was stored in a refrigerator until experiments. The control (C) included COMBO medium. The fish (F) and Chaoborus (CH) predation risk treatments included respectively 20- (i.e., 1 fish per 10 L) and 2.5-fold (i.e., 20 Chaoborus larvae per liter) diluted filtered medium stocks. Finally, the combination treatment (CH + F), representing joint predation risks, included both diluted medium stocks.

Life history experiment

The C. cornuta clone used in this experiment was sampled from Lake Taihu (31°22′13.548″N, 120°0′16″E), China. C. cornuta was cultured in COMBO medium at 25°C and 500 lx fluorescent light intensity under a 14:10 h light/dark cycle and fed Chlorella pyrenoidosa(1.5 mg C·L-1). Synchronous C. cornuta at a density of 1 individual per 10 mL, was adapted to the above conditions for at least two generations. To test the type of inducible responses of individual traits as well as the response strategy under joint predation risks, we set up a full factor experiment containing C, F, CH, and CH + F. We randomly divided newborn individuals into different treatments within 12 h. Each individual was cultured in 10 mL of medium, with 10 replicates per treatment. The media in different treatments were refreshed daily. Experimental workflow is presented in Figure S1.
Body size and horns were measured at maturity and at the late developmental stage (i.e., at 16th day). Horns were scored following the method described by Gu et al. (2021): absent (score 0), small (score 5), and large (score 10). Individual scores were normalized by a maximum score to define the induction level (0%–100%). In addition, time to the first brood, neonate size, brood number, total offspring number, and average brood size were recorded.

RNA extraction and sequencing

To accurately evaluate the type of response at the transcriptional level, we sequenced the transcriptomes of C. cornuta under C, F, and CH treatments. Groups of 250 newborn individuals were cultured in 2.5 L of medium with three replicates per treatment. During cultivation, responses triggered by different predation risks were verified through inducible traits, that is, horns and body size at maturity. The medium was refreshed daily, and samples were collected within 12 h after the first brood. The samples were frozen in liquid nitrogen and homogenized in TransZol Up. Total RNA was extracted using the TransZol Up Plus RNA Kit (ER501, TRANS, China), following the manufacturer’s instructions. RNA quality was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies, USA) and verified using agarose gel electrophoresis. The RNA integrity number of all samples was > 7.0.
The total mRNA of Ceriodaphnia was enriched by Oligo (dT) beads. The enriched mRNA was fragmented into short fragments using fragmentation buffer and reverse-transcribed into cDNA using random primers. The cDNA fragments were purified, end-repaired, poly(A)-added, and ligated to Illumina sequencing adapters. The ligation products were size-selected and PCR-amplified to develop a cDNA library. The cDNA library was sequenced using Illumina HiSeqTM 4000 by Gene Denovo Biotechnology Co. (Guangzhou, China).

Transcript assembly and annotation

Since genomic sequencing of Ceriodaphnia has not been conducted to date, we adopted de novo RNA-Seq to analyze the C. cornuta transcriptome. To obtain high-quality clean reads, sequenced reads were cleaned by removing reads containing adapters, > 10% unknown nucleotides (N), and low-quality reads (Q-value ≤ 20). Clean reads were assembled into unigenes using Trinity (Grabherr et al., 2011). To annotate the unigenes, we used the BLASTx program with an E-value < 10-5 in the NCBI non-redundant protein (Nr), SWISS-PROT protein, Kyoto Encyclopedia of Genes and Genomes (KEGG), and COG/KOG databases. Functional protein annotations were obtained according to the best alignment results.

Gene expression analysis and RT-qPCR

Unigene expression was quantified and normalized to reads per kb per million (RPKM) (Mortazavi, Williams, McCue, Schaeffer, & Wold, 2008). Differential transcript expression analysis between the control and different predation risk treatments was performed using DESeq2. Genes with a false discovery rate (FDR) < 0.05 and absolute fold change ≥ 1 were considered differentially expressed genes (DEGs). We categorized the DEGs into different types according to our classification and then used KEGG pathways for functional analysis. Pathways with P ≤ 0.05 were considered significantly enriched. The representative DEGs of various significantly enriched pathways were selected according to the following precedence conditions: stable expression, large fold change, and pathways annotated in a closely related species, that is, D. magna and Daphnia pulex .
To validate the RNA-Seq data in the C. cornuta transcriptome, we quantified the expression of 18 random DEGs using RT-qPCR. ddH2O was used as the negative control. cDNA was synthesized from mRNA using Reverse Transcriptase SuperMix (R233, Vazyme, China), and RT-qPCR was performed using ChamQ Universal SYBR qPCR Master Mix (Q711, Vazyme, China). All primer sequences are listed in Table S1. We obtained expression data for four alternative reference genes (Scoville & Pfrender, 2010) and calculated their average gene expression stability using geNorm. Glyceraldehyde-3-phosphate dehydrogenase (G3PD), RNA polymerase II (RNAP II), and elongation factor 1-alpha genes (EF) were determined to be stably expressed and geometrically averaged to calculate the gene expression normalization factor for each sample. Gene expression was quantified using the 2-ΔΔt method. The concordance between RNA-Seq and RT-qPCR data was assessed using regression analysis.

Statistical analysis

To test the effects of different predation risks, a MANOVA was performed on individual traits. When the residuals were normally distributed (Shapiro–Wilk test) and variances were homogeneous (Levene’s test), the data for each trait were evaluated by two-way ANOVA, followed by Bonferroni comparisons between different groups. When the normality test failed, the Scheirer–Ray–Hare test was used, followed by the Wilcoxon rank-sum test to assess significant differences among the different treatments. In the present study, statistical significance was set at P < 0.05. To simplify the effects of joint predation risks on different types of responses, we considered our results based on the conceptual diagram of possible prey responses in the combined predator treatment (Fig. S2). To test the correlations between different types of responses, we determined Pearson’s correlation coefficients between traits as well as between representative DEGs in different categories. The above statistical tests were performed in R (version 3.6.2). In addition, correlation, principal component, DEGs, and pathway enrichment analyses were conducted using OmicShare Tools (https://www.omicshare.com/tools/).

Results

Morphology and life history traits

The different predation risks significantly triggered various responses of morphology and life history traits (Table 1). Compared with the control, the responses induced by fish and Chaoborus predation risks could be classified into the following four categories (Fig. 1). (1) Unique responses to Chaoborus larvae [horn expression (at maturity and 16th day) and total offspring number]: C. cornutadeveloped horns (maturity: P < 0.001; 16th day:P < 0.001) and produced more offspring (P = 0.002) under Chaoborus predation risk, although these responses were not significant under fish predation risk. (2) Unique responses to fish (time to first brood and neonate size): the neonate size (P= 0.006) and time to first brood (P < 0.001) ofC. cornuta were significantly decreased under fish predation risk, although these responses were not significant underChaoborus predation risk. (3) General responses (size at maturity and brood number): the size (CH vs. C: P = 0.007; F vs. C:P < 0.001) and brood number (CH vs. C: P = 0.006; F vs. C: P = 0.017) were significantly decreased under both fish and Chaoborus predation risks. (4) Antagonistic responses (size at 16th day): C. cornuta size increased underChaoborus predation risk (P = 0.034) but decreased under fish predation risk (P < 0.001). Additionally, no significant differences were observed in average brood size. Thus, the unique responses in individual traits of C. cornuta to different predators were dominant: unique responses (five traits) > general responses (two traits) > antagonistic responses (one trait).