Transcriptomic analysis dissects the mechanistic insight into the Daphnia clonal variation in tolerance to toxic Microcystis

Cyanobacteria have become more prevalent than other phytoplankton in freshwater assemblages during summer. In such conditions, cyanobacterial traits may reduce zooplankton fitness and the energy flow efficiency from primary producers to aquatic herbivores. Cladocerans, as the dominant zooplankton grazers in freshwater ecosystems, exhibit clonal variation in their tolerance to cyanobacteria with an increasing gradient in eutrophication history. Hitherto, research on the full modes of action (MoAs) of Daphnia clonal differences in tolerance to toxic Microcystis still remains in its infancy. We conducted fitness and transcriptome analyses on two Daphnia clones, clone TH09 and TH14. A significant decline in body growth rate was detected in the sensitive clone TH09 compared with the tolerant clone TH14 at the presence of toxic Microcystis. Furthermore, transcriptional analysis indicated that major MoAs such as glutathione metabolism, protein processing in endoplasmic reticulum, amino sugar/nucleotide sugar metabolism, and arachidonic acid metabolism were linked to the tolerance fitness in Daphnia similoides. These results provided mechanistic insights into the pathways of genetic and biological processes involved in cyanobacteria tolerance in the Daphnia clonal variation, and propose that the genetic architecture of this fitness‐related trait would be helpful to clarify how zooplankton clones adapt to harmful algal blooms.

Phytoplankton composition in assemblages has shifted to increased cyanobacterial dominance, and therefore disrupt food-web processes during summer because of global warming and remarkable increases in nutrient fluxes (Kosten et al. 2012). As the most common bloom-forming cyanobacterium, Microcystis spp. can perform the biosynthesis of the hepatotoxin, microcystin (MC) (Carmichael 1997;Codd 2000). MC inhibits protein phosphatases (Toivola et al. 1994;Runnegar et al. 1995) and damages DNA through the promotion of oxidative stress (Žegura et al. 2003). In addition, Microcystis is known for containing protease inhibitors that act against the key digestive serine proteases in the alimentary canal of arthropod herbivores (Agrawal et al. 2001). Most herbivorous zooplanktons, such as Daphnia, are negatively affected by toxic cyanobacteria, as cyanobacteria could contribute to enhanced mortality, abnormal development, and lower reproduction (Lürling and van der Grinten 2003;Dao et al. 2010).
Of note, the uniform understanding of this negative relation between cyanobacterial and zooplankton abundance has been questioned in some lab-and field-scale investigations (Gustafsson et al. 2005;Sarnelle and Wilson 2005;Jiang et al. 2016;Lyu et al. 2016b), revealing that Daphnia had the ability to gain tolerance to dietary toxic Microcysits. The reproductive success of Daphnia is characterized by cyclical parthenogenesis, which is involved in alternation between sexual reproduction and parthenogenesis. Therefore, it is reasonable to assume that Daphnia populations consist of coexisting clones. Different Daphnia clones, that is, genotypes, exhibit differences in their individual fitness after exposure to cyanobacteria (reviewed by Ger et al. 2014). The initial evidence for clonal variation in tolerance to cyanobacteria was reported by Gilbert (1990), where an Anabaena strain inhibited the population growth rate in one Daphnia clone without affecting another. Another previous work studied individuals raised from diapaused eggs deposited in Lake Constance sediments and found that their tolerance to cyanobacteria was changed by the historical gradient of eutrophication (Hairston et al. 1999). Subsequent studies then indicated the improving *Correspondence: yangzhou@njnu.edu.cn This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Additional Supporting Information may be found in the online version of this article.
tolerance of Daphnia to toxic cyanobacteria was considerably dependent on clonal variation (Sarnelle and Wilson 2005).
Although we recognize the importance of gene function and annotation, we have marginal knowledge on their changes in response to ecological or environmental fluctuations, and subsequently their influence on phenotypic variation as well as macroscopic population dynamics (Miner et al. 2012). Therefore, it is essential to evaluate the transcriptional responses to changes in environment together with tractable and well-characterized ecologies (Colbourne et al. 2011). Under the scenario of toxic Microcystis, the discovery of the full modes of action (MoAs) of Daphnia clonal differences in tolerance are currently in its infancy. The previous studies of molecular cyanobacterial responses have mainly focused on the identification of candidate genes (Schwarzenberger et al. 2010(Schwarzenberger et al. , 2012(Schwarzenberger et al. , 2014Lyu et al. 2015). For example, Schwarzenberger et al. (2012) investigated that the inter-clonal variation of Daphnia tolerance to protease inhibitors in cyanobacteria strongly depended on the residual activity of digestive proteases. They subsequently revealed that ATP-binding cassette (a gene family for substrate transportation) RNA expression was linked with the microcystin tolerance in four Daphnia clones by feeding dietary supplementation with pure microcystin (Schwarzenberger et al. 2014). However, it remains unclear whether other potential novel MoAs exit in Daphnia tolerance to live cyanobacteria cells.
Since Daphnia, a global grazer of cyanobacteria, inhabits a central position in pelagic food webs (Ger et al. 2014), it is essential to directly reveal the fundamental MoAs of the potential of Daphnia to resist live cyanobacteria cells and therefore to possibly suppress cyanobacterial blooms (Sarnelle 2007). The present study aimed to demonstrate the important genes and the biological pathways related to Microcystistolerance in the Daphnia clonal variation. To compare the expression and observe the parallel body growth rates, which have been proposed as an appropriate proxy for Daphnia population fitness (Lampert and Trubetskova 1996), RNA-seq was executed separately to sample experimental animals from two Daphnia clones which were both fed with toxic Microcystis. Of note, the exposure to live Microcystis aeruginosa cells, rather than the waterborne cyanobacteria extract or microcystin, was chose in our study, because estimating the effects of dietary live Microcystis is helpful to obtain a more authentic knowledge of aquatic environmental influence (Sarnelle and Wilson 2005;Sarnelle et al. 2010;Asselman et al. 2012).

Organisms' origin and maintenance
The two clones of Daphnia similoides in the present experiment were sampled at Lake Taihu (120 09 0 E, 31 10 0 N; East China) and were named clone TH09 and clone TH14, respectively. This Daphnia species was widely distributed in fresh water in low latitudes (Makino et al. 2017;Rizo et al. 2017;Xu et al. 2018). The lake (i.e., Lake Taihu) experiences annual heavy cyanobacterial blooms, during which Microcystis constitutes over 50% of the summer phytoplankton biomass (Liu et al. 2011;Deng et al. 2014). Genetic distinctness was examined in the two clones by using a PCR method (Schwenk 1993) prior to the optimized modification. The two clonal cultures were established from one parthenogenetic mother respectively and were fed by high-quality food Scenedesmus obliquus (FACHB-416; 1.5 mg C L −1 per day). Only offspring from the third brood were used in culture and exposure. These clones grew at 25 C, under a 14 : 10 h light/dark cycle, and in M4 medium (Elendt and Bias 1990) and in the laboratory acclimated for more than 5 months in order to reduce possible disturb caused by maternal effects and environmental variance before the clone exposure. Meanwhile, the Daphnia medium was mildly aerated with clean oxygen for 1 d prior to use and renewed completely two times per week.
The high-quality alga S. obliquus and the toxic M. aeruginosa strain PCC7806 were purchased from the Freshwater Algae Culture Collection of the Institute of Hydrobiology at the Chinese Academy of Sciences. Cyanobacterium M. aeruginosa grew in the form of single cell (5-6 μm diameter spheres) or paired cells, so that it prevented mechanical interference from large-sized colonies that was probably to impact Daphnia ingestion. Furthermore, we executed a pilot experiment under proposal by the previous study (Gan et al. 2010), wherein we evidenced that the M. aeruginosa strain biosynthesize two kinds of MCs, namely, MC-LR as well as MC-RR. Each cell produced 3.6 pg of each MC variant. S. obliquus and M. aeruginosa were stored in 400 mL of BG-11 medium (Stanier et al. 1971). The algae were grown semi-continuously at 25 C under fluorescent light (40 μmol photons m −2 s −1 ) and in a 14 : 10 light/dark cycle. The strains were centrifuged in the exponential stage and fed to the zooplankton. The dietary algal densities were determined with a blood cell counting chamber under an inverted Nikon microscope at a magnification of ×400 (Lyu et al. 2015).

Experimental design
We conducted a body growth experiment to investigate the D. similoides tolerance performance, according to the previous study (Lampert and Trubetskova 1996). Three replicates were conducted and these animals were released within 12 h. Subsequently, 20 neonates from each beaker were transferred to 1 L of M4 medium. These neonates were exposed to the following food treatments for 7 d at a total food concentration of 1.5 mg C L −1 : (1) 100% S. obliquus; (2) a mixture of 50% M. aeruginosa + 50% S. obliquus. We monitored the medium parameter (Mean value: e.g., pH 7.8, dissolved oxygen = 7.0 mg L −1 , temperature 25 C and water hardness [as CaCO 3 ] of 276 mg L −1 ) every 4 h and renewed exposure medium daily, in order to keep parameters steady. Three biological replicates were analyzed for each food treatment, which was administered to Lyu et al. Tolerance difference in Daphnia clones each clone in a steady photoperiod of 14 : 10 h (light/dark) cycle and at 25 C. The body growth rates (BGRs) of the D. similoides clones were determined from the body dry weight of a subsample of the animals at the beginning and end of the experiment (Day 7), calculated based on the study (Schwarzenberger et al. 2014). Number of molting was examined by counting shed carapaces. Living D. similoides individuals were snap-frozen after the BGR measurements. In the present study, we compared the transcriptomic change between the D. similoides clone TH09 and TH14 fed with M. aeruginosa. Subsequently, the living individuals of the two clones exposed to MC-producing M. aeruginosa were subjected to RNA extraction and RNA-seq analyses.
RNA sampling and cDNA library construction D. similoides RNA was isolated by using TRIzol Reagent (Takara, Japan) in accordance with the manufacturer's instructions. The total RNA quality and quantity of the two groups of samples were evaluated by electrophoresis on 1.2% agarose gels, and the absorbance at 260/280 nm was determined using a spectrophotometer (Nanodrop, Thermos Scientifics). Only samples with OD 260 /OD 280 = 1.8-2.0, OD 260 /OD 230 ≥ 2.0 and RIN ≥ 8.0 were used for sequencing. Further details of RNA sampling and cDNA generation are provided in the Supporting Information.

Illumina sequencing, assembly, and annotation
The cDNA library was sequenced on the Illumina HiSeq 2500 sequencing platform that generated approximately 100 bp paired end (PE) raw reads by LC Sciences (Houston, U.S.A.). The raw sequence data were uploaded to the NCBI Sequence Read Archive (SRA). Further details of sequence analysis are provided in the Supporting Information.

Analysis of differentially expressed genes
To analyze the differentially expressed genes in D. similoides clones in the presence of toxic M. aeruginosa, the number of reads for each of the contigs was converted to reads per kilo base per million (RPKM) by using RSEM 1.2.3 (Wang et al. 2017). We used DESeq to determine the FDR (false discovery rate) threshold. If FDR was smaller than 0.05 and FPKM values showing at least twofold (log2FC = 1) difference among samples, this unigene was considered to indicate significant expression abundance (Wang et al. 2017). Differentially expressed genes (DEGs) were further annotated by Gene Ontology (GO) functional enrichment and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway analysis, as demonstrated in the Supporting Information.
Validation of the RNA-seq profiles by real-time quantitative PCR (qPCR) To validate the RNA-Seq results, 10 genes (see Table S4 in Supporting Information) were randomly selected to confirm the expression profiles of D. similoides clones in response to M. aeruginosa by qPCR. The D. similoides clones were treated as described above. qPCR was performed with an SYBR Green Master kit (Takara, Japan) according to the manufacturer's protocol. The primers for qPCR are listed in Supporting Information Table S4. qPCR was carried out in a total volume of 20 μL (10 μL of 2 × SYBR Premix Ex Taq, 1 μL cDNA mix, 0.5 μL of each primer [10 μM], and 8 μL of sterile distilled H 2 O). The PCR program was 95 C for 30 s, followed by 40 cycles of 95 C for 5 s and 60 C for 30 s. qPCR method with β-actin (see Lyu et al. 2016a for the primer sequences) as an internal control was used to explore mRNA expression levels.

Data analysis
Two-way (Daphnia clone difference and food treatment) ANOVA followed by Holm-Sidak test was employed to detect significant differences in BGR changes and number of molting. Data of BGR and number of molting satisfied assumptions of normality and equal variances before the ANOVA test. A significance level of p < 0.05 was applied. All statistical analyses were performed using SigmaPlot v11.0 (Systat Software).

Changes in BGR and number of molting
Significant BGR (p = 0.019) and number of molting (p < 0.001) differences between clones TH09 and TH14 exposed to M. aeruginosa were detected ( Fig. 1; Supporting Information Table S1), which indicated that clone TH14 was more tolerant to toxic M. aeruginosa than clone TH09. In contrast, the BGR of clone TH09 (0.63 AE 0.075) was significantly higher than that of clone TH14 (0.42 AE 0.013) ( Fig. 1A; p < 0.001) under good food conditions. For the clone TH09, the BGR in toxic M. aeruginosa treatment was significantly lower than that in Scenedesmus ( Fig. 1A; Supporting Information Table S1; p < 0.001). At the presence of toxic Microcystis, number of molting in the sensitive clone TH09 was significantly lower than that in the tolerant clone TH14 ( Fig. 1B; p < 0.001).

Transcriptome sequencing and assembly
To screen for the genes related to Microcystis tolerance, we isolated mRNA from clones TH09 and TH14 after 7 d exposure to the mixture containing 50% MC-producing M. aeruginosa and 50% S. obliquus. Two transcriptome libraries were constructed by using total RNA and subjected to Illumina deep sequencing. After quality control and clipping adapter, an amount of 86,978,248 clean reads remained, comprising 46,438,962 in clone TH09 and 40,539,286 in clone TH14 (Supporting Information Tables S2). These reads were transferred in the de novo assembly by the ABySS program. The assembly created 40,653 contigs with an average length of 3265 bp (N 50 = 1364 bp). After splicing and redundancy filtration, the contigs were intergraded into 17,738 unigenes, each of which had an average length of 1617 bp (N50 = 2658 bp). The unigenes length distribution was shown in Supporting Information Fig. S1. These assembly sequences represent transcription and can be used for further analysis.
After assembly, all the unigenes were analyzed. A cut-off Evalue of 10 −5 was used for BLAST analysis in the NCBI nucleotide sequences (NT), NCBI non-redundant protein sequences (NR), and Swiss-Prot databases (Supporting Information Table S3). The results showed that 12,598 (71.02%) unigenes attained significant BLAST hits in the NR database and 9552 (53.79%) hits in the Swiss-Prot database. Furthermore, 9446 unique protein hits were identified after searching the NR and Swiss-Prot databases. Only 4871 unigenes were not identified in the three databases.

Differential expression analysis
To identify whether the DEGs involved in D. similoides tolerant to toxic M. aeruginosa, pairwise comparisons for differential expression analysis were conducted between clones TH09 and TH14 fed by M. aeruginosa. In this pairwise comparison, 1864 genes were differentially expressed, comprising 689 upregulated genes and 1175 downregulated genes (Fig. 2).

GO and KEGG enrichment analysis of DEGs
In the present study, 1864 DEGs were assigned GO terms. The significant GO terms divided into the three fields of biological process, cellular component, and molecular function were shown in Fig. 3. Among the biological processes, 169 (25%) of the DEGs were distributed to apoptotic process (GO: 0006915), and 107 DEGs were distributed to immune response (GO: 0006955). In the field of molecular function, approximate 20% of the DEGs were belonged to protein homodimerization activity (GO: 0042803) function, with another 19% involved in calcium ion binding (GO: 0005509). In cellular component, the high-enriched GO terms were plasma membranes (GO: 0005886, 23.8%) and integral to membrane (GO: 0016021, 23.7%).
The analysis of pathway enrichment was helped to deeply recognize biological functions of the DEGs. In the present study, five significant (p < 0.05; Table 1) pathways were identified through the comparison of clones TH09 and TH14, including glutathione metabolism, protein processing in endoplasmic reticulum, amino sugar/nucleotide sugar metabolism and arachidonic acid metabolism.

Validation of RNA-seq by qPCR
To validate the DEGs identified by RNA-seq, the expression of 10 target genes randomly selected from those with differing expression patterns and from genes of interest based on functional enrichment and pathway results, was confirmed by qPCR (Table 1). The melting-curve analyses revealed a single product for all tested genes. As shown in Fig. 4, all qPCR results were in good agreement with the Illuminia RNA-seq analysis. For example, the expression of genes (glutathione Stransferase D2, gamma-glutamyltransferase 1, serine/threonine-protein kinase PDIK1L, chitinase 10, chitin synthase 2, hematopoietic prostaglandin D synthase, and group XIIB secretory phospholipase A2) showed an immediate significant decline in clone TH09 compared with TH14 in the presence of toxic Microcystis.

Discussion
According to the hypothesis of the arms race, the degree of outbreak and spread of cyanobacterial blooms relies extensively on the dynamic interplay between the capacity of the toxic algae to fend off the herbivore consumers and the resistance ability of the herbivore consumers to tolerate the toxicity. On the one hand, the toxicity and morphology are considered as key traits in the reduction of grazing losses of cyanobacteria by zooplankton (Yang et al. 2006;Jang et al. 2007;Pineda-Mendoza et al. 2014), although cyanotoxins have other various functions, including iron chelator (Utkilen and Gjølme 1995), photosynthesis or other light related processes (Young et al. 2005), and intercellular intraspecies communication (Schatz et al. 2007;Zilliges et al. 2008). On the other hand, zooplankton often exist mutually with harmful algal blooms food treatment containing 100% S. obliquus or food treatment containing 50% S. obliquus and 50% M. aeruginosa. The symbol (★) denoted significant difference (p < 0.05) between the two clones (clone TH09 and clone TH14) in each food treatment, according to two-way ANOVA.
(HABs), suggesting that they have multiple abilities to behaviorally reduce the ingestion of hazardous doses and to be physiologically tolerant to the ingested toxin. In line with the previous studies (Hietala et al. 1997;Chislock et al. 2013;Jiang et al. 2013Jiang et al. , 2015, the present population-level results clearly showed clonal variation in Daphnia fitness in response to toxic Microcystis. Of note, the clones were cultured individually under identical conditions. Dissolved oxygen, pH, and temperature in the medium were steadily maintained as exposure medium were renewed daily. Therefore, it is conceivable that significant BGR differences between the two clones would be arise from genetic factors, rather than their inducible tolerance. The present clones were originated from individual D. similoides adults isolated from the cyanobacterial-polluted lake, which were rational considered as DNA-level distinct, despite DNA traits in each clone were not monitored by genomic analyses.
Live cyanobacterial cells have been observed to decline zooplankton fitness through the following three basic routines: (1) large colony formation and filamentous morphologies can block filtering appendages and thereby inhibit grazing (DeMott et al. 2001); (2) nutritional deficiency in the cyanobacterial cells (Von Elert and Wolffrom 2001); and (3) cyanobacteria produce toxic secondary metabolites (e.g., MC and digestive serine proteases) (Rohrlack et al. 2001;Yang et al. 2012). The frequently observed mechanical clog on filtering apparatus (DeMott et al. 2001) in Daphnia did not appear in the present study as Microcystis we used was mostly consisted of single and paired cells. Hence, the clonal variation in BGR against toxic Microcystis results from toxic metabolites and nutritional deficiency. We then attempted to acquire valuable global-scale transcriptomic evidence by comparing the clones TH09 and TH14 at the presence of dietary Microcystis, in order to support the clonal variation of Daphnia in response to cyanobacterial effects. Given the advantages in GO annotation and KEGG pathway classification, we identified and validated certain sensitivity-related genes, and the most important pathways involved clonal variation in response to toxic Microcystis were further discussed below. Due to the limited scale of the present study, we encourage future research incorporating more treatments, temporal endpoints across zooplankton taxa, and improved gene annotation to enrich the dynamic processes of MoAs of zooplankton tolerance to toxic cyanobacteria.

Glutathione metabolism
Intracellular MCs suppress protein phosphatases (An and Carmichael 1994), which, along with generation of oxidative stress, is treated as the major MoAs of their toxicity (Rohrlack et al. 2004;Lyu et al. 2014). A possible mechanism for buffering the toxic effects of MCs is linked with glutathione Stransferases (GSTs), which catalyzes the biotransformation of MCs with glutathione; this process is widely proposed as the first step in detoxification found in many aquatic livings, including plants, invertebrates, and fish (Pflugmacher et al. 1998;Ortiz-Rodríguez and Wiegand 2010;Lyu et al. 2016a).
Our study further revealed that glutathione metabolism was significantly (p = 0.0049) enhanced in the tolerant clone TH14 as GST D2, aminopeptidase N, thyrotropin-releasing hormone-degrading ectoenzyme (TRH-DE) and gammaglutamyltransferase 1 (GGT1) expression were significantly (p < 0.001) increased compared to the sensitive clone TH09 Fig. 2. Hierarchical cluster analysis shows toxic Microcystis induced gene expression signatures in clone TH09 and clone TH14 (|log2FC| ≥ 1, p < 0.05). Green, red, and black areas indicate, respectively, decreased, increased, and no significant change in clone TH09 compared with clone TH14.
( Fig. 5; Table 1). Based on the finding that antioxidation as well as biotransformation are the two main processes of GST in MCs detoxification (Pflugmacher et al. 1998), it is conceivable that the enhancement in GST D2 gene expression in the tolerant clone compared with the sensitive was observed when the animals are fed with the toxic Microcystis. Of note, another important element in the detoxification process is glutathione (GSH). Given the GSH depletion caused by conjugation with MC, the cells must continue with GSH re-synthesis to meet the detoxification requirements (Ding et al. 2000). In our study, aminopeptidase N and GGT1, which are involved in GSH synthesis (Lee et al. 2004), were up-regulated in the tolerant clone TH09; thus, GSH synthesis was enhanced. This suggested that these effects inhibited the MC toxicity in the cell.

Protein processing in endoplasmic reticulum
In the pairwise comparison of clone TH09/TH14, the expression levels of serine/threonine-protein kinase-PDIK1L (STPK-PDIK1L) and tyrosine-protein kinase (TPK) were significantly (p < 0.001) reduced, whereas that of heat shock protein 70 Bbb (HSP70B) was significantly increased (p < 0.001). In protein processing in endoplasmic reticulum, the STPK-PDIK1L serves as a sensor to detect the unfolded protein response (Fig. 5). Furthermore, TPK can phosphorylate cytosol or nucleus proteins at tyrosine residues during the process (Radha et al. 1996). These enzymes play an essential role in the homeostatic maintenance of stress response and protein folding, which are crucial when the number of unfolded proteins increases (Walter and Ron 2011). As reported previously (Asselman et al. 2012;Lyu et al. 2016c), exposure to toxic Microcystis inhibited STPK-PDIK1L and TPK expression in the sensitive clone TH09 and interfered with normal protein folding and repair. In addition, since HSP70B serves as a molecular chaperone, minimizes protein aggregation, and repairs and protects cellular proteins from stress damage (Mayer and Bukau 2005), it was to be expected that increased HSP70B expression was induced by an abundance of unfolded proteins. Accordingly, the protein processing in endoplasmic reticulum pathway was identified as statistically significant (Table 1), indicating that the repair of dysfunctional protein folding was enhanced in the tolerant clone.
Amino sugar and nucleotide sugar metabolism While the inhibited chitinase 10, beta-N-acetylglucosaminidase and chitin synthase 2 in clone TH09 compared to clone TH14, the increased amino-sugar and nucleotide-sugar metabolism pathway was found in the tolerant clone TH14 (Table 1). Microcrustaceans growth and immunity require molting in some species (Duneau and Ebert 2012).The present study found the number of molting in the sensitive clone TH09 was significantly lower than that in the tolerant clone TH14, at the presence of toxic Microcystis (Fig. 1B), suggesting that Microcystis disordered molting in sensitive clone worse than the tolerant clone. The prerequisite events during molting are the aged cuticle separation from the underlying epidermis and an increased levels in several enzymes and hormones, collectively named as molting fluid (Stevenson 1972). During molting, chitinase digests the chitin in the old cuticle and chitin synthase catalyzes the generation of new cuticles (Fig. 5). The opposite functions of chitin synthase 2 and chitinase 10 suggest that their activities must be highly coordinated in integument during the moltingcycle (Bade and Stinson 1978). In addition, beta-Nacetylglucosaminidase acts synergistically with endochitinases, which cleave crystalline chitin into chitooligosaccharides. They may also alleviate the inhibition of chitinases through the accumulation of chitotriose and chitotetraose in the molting fluid (Fukamizo and Kramer 1985). In our study, these enzymes were up-regulated in the tolerant clone TH14, which directly explaining that the tolerant clone had greater superiority in individual fitness than the sensitive clone ( Fig. 1).

Arachidonic acid metabolism (AAM)
Living cyanobacteria contain lipopolysaccharide (LPS) in their outer cell layers. Previous studies have revealed that bacterial LPS, a potentially inflammatory toxin (Mayer et al. 2011), increased the release of arachidonic acid and its metabolites (Rosenberger et al. 2004). Here, several AAMresponsive genes (i.e., hematopoietic prostaglandin D synthase (H-PDS), group XIIB secretory phospholipase A2 (GXSPA2), and gamma-glutamyltransferase 5) were significantly (p < 0.001) induced in the tolerant clone TH14 compared with the sensitive clone ( Fig. 5; Table 1). H-PGDS requires GSH for the reaction; it is inactive with other thiol compounds. As the inhibition of GSH production in the sensitive clone was witnessed in the present study, we cannot exclude the possibility that the down-regulated H-PDS gene expression resulted from feedback regulation by adequate H-PDS protein. That is, the increase in H-PGDS expression in the tolerant clone TH14 resulted in enhanced prostaglandin formation (Dawson 1998), and potentially mediated inflammation triggered by microcystin. In addition, GXSPA2, as an atypical member of the secreted phospholipase A 2 family, catalyze the hydrolysis of glycerophospholipids at the S n 2 position to release arachidonic acid and other free fatty acids (Schaloske and Dennis 2006). Based on these reports, we speculated that the enhancement of AAM in the tolerant clone TH14 contributed to the development of resistance against cyanobacterial LPS inflammation.

Conclusions
Through the application of a transcriptome-wide top-down strategy, we had access to recognize the genes and complicated pathways involved in cyanobacteria tolerance in Daphnia populations. As expected, not only did we discover evidence compatible with previous results (Asselman et al. 2012(Asselman et al. , 2017Schwarzenberger et al. 2014), but we also identify novel candidate pathways that accounted extensively for the tolerance of Daphnia to cyanobacteria. There was a total of 1864 DGEs that highlighted the pathways of glutathione metabolism, protein processing in endoplasmic reticulum,

Lyu et al.
Tolerance difference in Daphnia clones amino sugar/nucleotide sugar metabolism, and arachidonic acid metabolism (Fig. 5). These significantly changed pathways suggested the strong tolerability of clonal variation was involved in cyanotoxin detoxification, the activation of protein fold repair, inflammation mediation, and molting/growth maintenance (Fig. 5). Our results emphasized the important roles of novel candidate pathways, which have often been missed in previous zooplankton-cyanobacteria studies.

Data accessibility statement
Sequence data is uploaded to the Sequence Read Archive (SRA) with accession PRJNA446001.

Author Contributions
KL conceptualized and designed the experiment together with ZY. Experiments as well as data analyses were conducted by KL, LG, HW, XZ, LZ, YS, and HY. The first draft of the manuscript was written by KL and commented by ZY. All authors have read the final version of the article.