Transcriptomic response of Daphnia magna to nitrogen‐ or phosphorus‐limited diet

Abstract Effects of nutrient‐imbalanced diet on the growth and fitness of zooplankton were widely reported as key issues to aquatic ecology. However, little is known about the molecular mechanisms driving the physiological changes of zooplankton under nutrient stress. In this study, we investigated the physiological fitness and transcriptomic response of Daphnia magna when exposed to nitrogen (N)‐limited or phosphorus (P)‐limited algal diet (Chlamydomonas reinhardtii) compared to regular algae (N and P saturated). D. magna showed higher ingestion rates and overexpression of genes encoding digestive enzymes when fed with either N‐limited or P‐limited algae, reflecting the compensatory feeding. Under P‐limitation, both growth rate and reproduction rate of D. magna were greatly reduced, which could be attributed to the downregulated genes within the pathways of cell cycle and DNA replication. Growth rate of D. magna under N‐limitation was similar to normal group, which could be explained by the high methylation level (by degradation of methionine) supporting the body development. Phenotypic changes of D. magna under nutrient stress were explained by gene and pathway regulations from transcriptome data. Generally, D. magna invested more on growth under N‐limitation but kept maintenance (e.g., cell structure and defense to external stress) in priority under P‐limitation. Post‐translational modifications (e.g., methylation and protein folding) were important for D. magna to deal with nutrient constrains. This study reveals the fundamental mechanisms of zooplankton in dealing with elemental imbalanced diet and sheds light on the transfer of energy and nutrient in aquatic ecosystems.


| INTRODUC TI ON
Ecological stoichiometry, a branch of ecology that considers how the balance of energy and elements influences the living systems (Sterner & Elser, 2002), has been widely studied in both land and aquatic systems where stoichiometric imbalances between predator and prey can affect trophic interactions and ecosystem functions.
Elemental compositions (e.g., C:N:P, short for carbon:nitrogen:phosphorus) of phytoplankton can vary easily and frequently following hydrographic conditions of ambient water. Nutrient-limited algae usually are deficient in certain constituents (e.g., certain lipids and amino acids) and are considered as low-quality prey for herbivorous zooplankton (lack of food preference; Weers & Gulati, 1997).
The performances (e.g., ingestion, growth, and reproduction rates) of zooplankton in response to low-quality diet, which have consequences for both phytoplankton and zooplankton population dynamics and energy flow, are important to the understanding of ecological stoichiometry in aquatic systems.
"All life is chemical," and certain chemical elements are critical to organisms for building and maintaining basic biological structures and core functions (Jeyasingh et al., 2011;Sterner & Elser, 2002).
Transcriptomics (expressed mRNA) can be used as a tool to understand the molecular mechanisms underlying many biological processes related closely to ecosystem functions. It can not only explain the physiological responses of organism under environmental stresses, but also provide a global expression pattern which may lead to a more detailed investigation of some "key genes" (Wang et al., 2009). However, transcriptomic studies are much less applied to the ecology of zooplankton (Lenz et al., 2021), compared with other organisms such as bacteria, phytoplankton, land plants, and vertebrates (Hua et al., 2004;Nguyen et al., 2015).
Daphnids, which filter particles in the size range of 0.5-50 μm including algae and bacteria, are key aquatic herbivores and model species to study aquatic ecology. Daphnids may often experience nutrient-limited diet coming from (1) intracellular biochemical difference among ingested food, for example, green eukaryotic algae, cyanobacteria, and fungi, (2) movement of daphnids between heterogeneous food patches, and (3) seasonal succession in the taxonomic composition of phytoplankton (Koussoroplis et al., 2017).
Much work has been done to study the phenotypic response of daphnids, especially on the growth and reproduction rates, to the nutritional constraints induced by diet, which showed that food elemental compositions (particularly N and P) have profound effects on daphnids' population dynamics. For instance, P-deficient food can significantly reduce the growth rate and egg production of daphnids, and the growth rate of daphnids correlates well with algal P content (Sterner et al., 1993;Weers & Gulati, 1997). N-limited diet can also affect the performance of daphnids, but the negative effects were not as severe as P-limited algae in some studies (Sterner et al., 1993;Weers & Gulati, 1997). However, except two studies reporting the transcriptomic response of daphnids to P-limitation using gene microarrays (Jeyasingh et al., 2011;Roy Chowdhury et al., 2015), the metabolic mechanisms (in terms of gene regulation and involved pathways) of daphnids' physiological changes in response to nutrient limitation, especially the different responses between N-and Plimitation, remain unclear.
In this study, we sequenced the mRNA extracted from Daphnia magna (one kind of daphnids species widely distributed in fresh water) fed by N-or P-limited green algae (Chlamydomonas reinhardtii) for 7 days and compared that to those on normal diet (N-and P-sufficient C. reinhardtii). Through gene and pathway enrichment analysis, we found that D. magna upregulated the expression of genes encoding digestive enzymes under both N-limited and Plimited diets, which explains their higher ingestion rate and compensatory feeding behavior. Their growth and reproduction were severely impaired under P-limitation, which can be attributed to the downregulation of genes controlling cell cycle and DNA replication.
Under N-limitation, gene regulation related to methionine metabolism maintained the body growth but reduced reproduction of D. magna. Our results suggest that post-translational modifications (e.g., methylation and protein folding) were important for D. magna to deal with nutrient imbalance and stress.

| Culture of algae
Green algae C. reinhardtii (CC1690) were grown in 1-L vessels with BG11 medium (Stanier et al., 1971). P-limited algae were obtained by reducing K 2 HPO 4 by 10-fold from 40 to 4 mg/L, while N-limited algae were obtained by reducing the concentration of NaNO 3 by 30fold from 1.5 to 0.05 g/L. No modification of medium was made for the normal group. Algae were cultured in a chamber under constant temperature of 23.5℃ and a 14:10 light/dark cycle (20 μmol/m 2 s −1 ).

Concentration of algae was measured daily by a Beckman Coulter Z2
Cell and Particle Counter (Beckman Coulter Inc.). Cellular carbon (C) and nitrogen (N) of algae in three conditions (i.e., normal, N-limited, and P-limited) were quantified by a FlashSmart CHNS Elemental Analyzer (Thermo Fisher). Cellular P was analyzed as orthophosphate after acidic oxidative hydrolysis with 1% HCl. The concentration of PO 4 3− was measured manually according to a previous report (Strickland, 1968).

| Measurement of ingestion, growth, and reproduction rates of D. magna
Before this experiment, D. magna was cultured in the creek water and fed with sufficient C. reinhardtii. During the experiment, newborns of D. magna were collected within 12 hr after birth and transferred to 3 glass beakers (volume = 5 L) with 0.45 μm filtered fresh water. Three kinds of algae (i.e., N-limited, P-limited, and nutrientsufficient C. reinhardtii as described above) were added to each beaker, respectively, as prey for D. magna. Ingestion rate, growth rate, and reproduction rate (numbers of egg in the brood and daily newborns) of D. magna in each treatment were measured on the 1st, 3rd, 5th, and 7th days during the experiment.
For ingestion rate measurement, on each sampling day, 10 individuals of D. magna were taken from the beaker and put in a 600-ml polycarbonate bottle (Nalgene), with three replicates. Algae concentration (as food) was kept constantly to 100,000 cells/ml (i.e.,

~2
.5 mg C L −1 ) which was above the incipient limiting level according to a pre-experiment. Grazing experiment was conducted for 24 hr, and ingestion rate (I, cells daphnids −1 day −1 ) was calculated following the previous studies (Frost, 1972): where Ct ′ and Ct (cells ml −1 ) are the prey (C. To measure growth rate, on each sampling day, 15 individuals of D. magna were randomly taken from the beaker for each treatment. Body length was measured from the anterior margin of the eye to the base of the tail spine of D. magna using a microscope (Olympus IX51). Daily increase in body length was used to represent growth rate.
To measure reproduction rate, 10 individuals of D. magna (with triplicates for each treatment) were taken from the beaker on the first day and cultured in a 600-ml polycarbonate bottle (Nalgene) for 7 days with the same medium in the beaker. On each sampling day (i.e., 1st, 3rd, 5th, and 7th days), a number of eggs and newborns in the brood pouch of D. magna were calculated using an inverted microscope (Olympus CK30). Newborns were removed immediately after calculating. Daily increase in egg production and newborns was used as reproduction rate.

| RNA extraction
On the 7th day, for each type of diet, 50 individuals (rinsed by Milli-Q water for three times. No starvation step was conducted as it will cause extra transcriptomic response) were taken from the beaker and pooled as one sample. Triplicate samples were collected for each type of diet. Body tissues were grinded manually in a PCR tube and kept in RNAlater solution at −80℃. RNA was isolated using TRIzol reagent (Takara) and RNA Mini kit (Thermo Fisher) in accordance with the manufacturer's instructions.
Genomic DNA was removed using the Turbo DNA-free DNase kit (Ambion), and mRNA was isolated with the MicroPoly (A) Purist kit (Ambion). RNA integrity was evaluated using an Agilent 2100 bioanalyzer (Agilent Tech.).

| Sequencing and transcript assembly
Each replicate mRNA was sequenced with Illumina HiSeq2500 (Novogene), generating 150-bp paired-end reads. Raw reads were processed sequentially by quality control, reads filtration, assembly, gene prediction, quantification, and sorting. FastQC (Babraham bioinformatics) was used to assess the quality of raw reads and quality control (length >140 bp, without ambiguous "N," and average base quality >20; Andrews, 2010). SRC_c method (Meng et al., 2018) was employed to sort reads affiliated to D. magna from the whole metatranscriptome (mainly include reads from D. magna, bacteria, fungi, and C. reinhardtii) with indexed k-mers set to 32 and a suggested default similarity value (50%). For this step, library of D. magna was built using RNA-seq datasets from a previous study (Orsini et al., 2016).
Gene showing |log 2 (fold change)| > 1 and p-value < .05 (quasilikelihood F test) at any of the two comparison groups with triplicated gene expression data (N-limitation vs. normal and P-limitation vs. normal) was defined as DEG (Lyu et al., 2019). DEGs were further annotated by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) using DIAMOND (v0.9.21.122) with the following parameters: blastp; k parameter = 1; e-value = 10 -7 (Buchfink et al., 2015). Higher level of GO terms and KEGG pathways were enriched by DEGs using TBtools (Chen et al., 2020), and corrected Benjamini p-value (i.e., q-value) < .05 as cutoff was used to define significantly enriched KEGG pathways.

| Validation of RNA-seq by qPCR
To validate the RNA-seq results, 7 DEGs, mainly within the enriched KEGG pathways, were randomly selected to confirm the expression profiles in our study. Primers (final concentration = 10 mM) were designed by Primer Premier v6.00. For each sample, 500 ng of extracted RNA was used for reverse translation by HiScript III RT SuperMix for qPCR kit (Vazyme Biotech). After the synthesis of cDNA, qPCR was performed with FastStart Universal SYBR Green Master (Roche) on LightCycler 384 (Roche). The qPCR followed an initial hold at 50℃ for 2 min and 95℃ for 10 min followed by 45 cycles of 95℃ for 15 s and 60℃ for 1 min. The relative abundances of genes were calculated by a 2 −ΔΔCt method, with β-actin as an internal control (Lyu et al., 2019).

| Statistical analysis
Statistical analysis was conducted in R software (v 4.0.0). Welch two sample t test was used to test the significant differences in ingestion rate, body length, and numbers of eggs in the brood pouch of D. magna on each sampling day (i.e., day 1, day 3, day 5, and day 7

| C:N:P of algae as food for D. magna
Manipulation of media nutrient contents yielded significantly different elemental composition in C. reinhardtii. Cells grown in N-limited medium had the lowest N content (1.95% ± 0.15%, atomic ratio) and higher C:N ratio (16.7 ± 0.1), while cells grown in the P-limited medium had the lowest P content (0.12% ± 0.01%, atomic ratio) and higher C:P ratio (467.7 ± 89.8), as compared with normal culture (C:N = 7.5 ± 0.2, C:P = 65 ± 6.1, p < .01, by Welch two sample t test; Table A1).

| Physiological characteristics of D. magna under different diet conditions
During the 7-day experiment, D. magna with P-limitation diet had the shortest body length compared to that with N-limitation or nutrientsufficient diet ( Figure A1a). After 3 days, D. magna exhibited higher ingestion rate but lower egg production rate when fed with N-or Plimitation diet than that with nutrient-sufficient diet ( Figure A1b,c).
Significant difference of the physiological performances of D. magna between two treatments at each sampling point is shown in Table A2.
For the entire experiment period, taking ages into consideration, both N-and P-limitation had significant influences on the ingestion and egg production rates of D. magna (Table A3). P-limitation had a significant effect on the growth rate (i.e., slow increase in body length) of D. magna, while this effect was weak for D. magna under N-limitation (Table A3; Figure A1). Here, we focused on the comparison between nutrient-limited group and nutrient-sufficient group on the 7th day when transcriptomic samples were collected. As shown in Figure 1 Figure 1b). An average of 20 ± 9.00 eggs were produced for each D. magna fed with normal algae on the 7th day, which was much higher than D. magna fed N-limited algae

| Transcriptome sequencing and assembly
Transcriptome libraries were constructed using total mRNA and subjected to Illumina deep sequencing (about 18 G per sample).
After quality control, clean reads were assembled to contigs: number of contigs in each sample ranged from 227,920 to 306,664, with N 50 ranged from 1,459 to 1,597 bp. In total, 19,815 unique genes (10,271, on average) of D. magna were predicted after sorting and redundancy filtration, representing 89.2% of the total genes in the transcriptome libraries (Table A4).
PCA using all unique genes showed that the biological replicates of each treatment (N-limitation, P-limitation, and normal) were close to each other and far from other treatments, supporting the reliability of RNA-seq data. Difference between each experimental group and control group was significant by analysis of similarity (ANOSIM) test (p < .01), suggesting the significant effects of low-quality diet on the gene expression of D. magna (Figure 2a,b).

| GO and KEGG enrichment of DEGs
All DEGs in this study were assigned to GO terms with three categories: biological process (BP), molecular function (

| Validation of genes by qPCR
To validate the DEGs identified in the RNA-seq data in this study, the expression of seven selected DEGs (see detailed information in Tables A5 and A6), which were mainly from enriched KEGG pathways, was confirmed by qPCR. As shown in Figure A2

| Effects of low-quality diet on D. magna physiology: why the effect of P-limited diet is more severe?
Our results were in accordance with previous studies reporting the effects of low-quality diet on the growth and reproduction of daphnids, with a conclusion that phosphorus-limited algae is a much poorer food for D. magna (Sterner et al., 1993;Weers & Gulati, 1997). is also required in other important processes, including synthesis of polyamines which are essential for various cellular functions affecting growth and development (Hayashi et al., 2018). Another upregulated gene was adenosylhomocysteinase (achY, EC:3.3.1.1) which hydrolyzes S-adenosylhomocysteine (SAH) to homocysteine (Hcy).
The overexpression of the two genes (metK and achY, from both qPCR and RNA-seq results) would increase the ratio of SAM:SAH which is frequently used as an indicator of cellular methylation capacity (higher ratio means higher capacity; Caudill et al., 2001).
Methylation is an important way for Daphnia to response and adapt to environmental stress (Asselman et al., 2017). It is also closely related to cell growth and tissue differentiation of metazoans and could be induced by nutrient restriction (Kusari et al., 2017). Therefore, in our study, the increased degradation of methionine and SAM: SAH ratio may promote the methylation level and stimulate the development (or growth) of D. magna with N-limited diet.
Methionine is necessary for the increase in fecundity of arthropods (e.g., fruit fly and copepods) and plays an important role in controlling the lifespan of animals (Grandison et al., 2009). Methionine is easy to be degraded during the hydrolysis process. Therefore, decrease of methionine, at the same time, may lead to the reduced reproduction of D. magna under N-limitation in our study.

F I G U R E 4
Framework showing the biological pathways of D. magna influenced by low-quality diet. Incoming food from environment is ingested, digested, and assimilated into two nutrient reserves (N and P) which will support the performance (maintenance, growth, and reproduction) of D. magna. Five significantly enriched KEGG pathways are illustrated to explain the changes in the performance of D. magna under nutrient constrains, with names of pathways in black solid lines. DEGs (some representatives) within these pathways were in green boxes. Detailed information of these DEGs is shown in Table A5. ER: endoplasmic reticulum; OST: oligosaccharyltransferase complex; SAM: S-adenosylmethionine; SAH: S-adenosylhomocysteine; SSB: single-stranded DNA-binding proteins; RPA: replication protein A; and Pol α (or ε): DNA polymerase α (or ε). G1, S, G2, and M are the four stages of cell cycle Secondly, D. magna fed with P-limited algae had the lowest growth rate and fewest number of newborns in our study, supporting that P-limitation causes much severer damage than N-limitation for daphnids (Sterner et al., 1993;Weers & Gulati, 1997). Our results showed that DEGs (annotated by KEGG) involved in DNA replication and cell cycle were all downregulated in D. magna fed with P-limited diet. This was in accordance with a previous study where genes related to nuclear structure, replication, recombination, and repair were downregulated for daphnids fed with low-P diet (Roy Chowdhury et al., 2015). Inhibition of P-limitation on nucleic acid metabolism and cell division could strongly affect the growth rate of phytoplankton, yeast, and bacteria (Brauer et al., 2008;Feng et al., 2015;Vaulot et al., 1996). Thus, our result suggested similar mechanisms may exist in zooplankton, leading to low growth rate of daphnids' population.
Cell cycle in metazoans is controlled by a number of mechanisms on the gene level. Effects and mechanisms of nutritional limitation on cell cycle and growth rate were well demonstrated for phytoplankton and yeast in previous studies (Brauer et al., 2008;Vaulot et al., 1996), while little is known about zooplankton. For instance, the decrease of growth rate was closely related to the number of cells arrested at the G0/G1 phase for yeast (Brauer et al., 2008) and G2 + M phase for marine diatom Thalassiosira pseudonana (Claquin et al., 2002). In our study, for D. magna under P-limitation, downreg-

| How do Daphnia deal with stoichiometric constrains? Lessons from transcriptomic data
We employed a modified framework of the dynamic energy budget model (Sperfeld et al., 2017) to track metabolic pathways of elemental nutrients from the aspects of maintenance, growth, and reproduction. To investigate how D. magna deal with diet-induced elemental constrains, we focused on the investment of nutrients to the three above aspects and related molecular metabolisms. We mainly focused on the genes within the significantly enriched pathways although other central pathways (e.g., carbon metabolism) are important and essential in dealing with nutrient constrains for daphnids.
Increasing in feeding rate, which is also called compensatory feeding, is one strategy for herbivorous zooplankton to deal with specific nutritional constrains in diet. While previous study showed that daphnids increased their ingestion rate when fed low-quality algae (Plath & Boersma, 2001), no significant difference in feeding rates between daphnids fed P-sufficient and P-depleted algae was also reported (DeMott et al., 1998). These contrasting results suggest that compensatory feeding varies with changes in diet including food abundance, digestibility, and the elemental ratio in food.
In our study, ingestion rate of D. magna fed with either N-limited or P-limited algae was higher on the 7th day than that fed with normal diet, production. This result is in agreement with the report that copepods under N-limitation were unable to utilize dietary N efficiently for egg production due to the N demands for maintenance (including growth; Kuijper et al., 2004).
Under P-limitation, although both growth rate and reproduction rate decreased significantly, a remarkable number of upregulated genes of D. magna were detected, which may represent the mechanisms in dealing with phosphorus constrains. Regulation of protein process in ER is one of the post-translational modifications (PTMs, e.g., phosphorylation, monoubiquitination, and glycosylation) which are important ways for organisms to enhance nutrient acquisition and utilize efficiency (Plaxton & Shane, 2018 (Shemi et al., 2016). The upregulation of genes within post-translation process and cell structure (e.g., membrane) biogenesis was also mentioned in a study comparing ancient and modern daphnids under P-limitation (Roy Chowdhury et al., 2015).
Nevertheless, the underlying P-related molecular pathways and gene regulation mechanisms need further investigation.
A recent study showed that low-quality diet (cyanobacteria) could reduce the output of Daphnia gut parasites (Manzi et al., 2020).
Besides, food with high C:P ratios can significantly reduce bacterial infection rates in daphnids (Frost et al., 2008). Similarly, in our study, Overall, in this study, using transcriptomic data, we investigated the metabolic mechanisms regulating the physiological changes of D. magna under nutrient stress caused by low-quality diet. We attributed the increased ingestion rate (i.e., compensatory feeding) of D. magna under nutrient limitation to the overexpression of genes encoding digestive enzymes. We found that severe negative impacts of P-stress on D. magna (e.g., low growth and reproduction rates) can be explained by the downregulation of genes corresponding all phases of cell cycle. By tracking metabolic pathways of elemental nutrients from the aspects of grazing, maintenance (e.g., post-translational modification), growth, and reproduction, we showed how some key genes regulated these biological processes that are important to the demography of both zooplankton and phytoplankton and ecosystem functions. However, the genes (and involved pathways) focused in our analysis are a very small percentage of the total number of the annotated genes. Other enriched metabolic pathways in our study, such as retinol metabolism which is essential for the reproduction of zebrafish and mammals (André et al., 2014), are well studied and proven important for the development and homeostasis maintenance of vertebrates but remain unclear for most invertebrates. Their relationships with the phenotypic responses of D. magna to nutrient stress require further investigation. Zhang for suggestions on nutrient concentration measurement.

CO N FLI C T O F I NTE R E S T
All authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Transcriptome sequencing data were deposited in GenBank (Sequence Read Archive) and are available under the BioProject PRJNA597965.