Genome-wide transcription profiles reveal genotype-dependent responses of biological pathways and gene-families in Daphnia exposed to single and mixed stressors.

The present study investigated the possibilities and limitations of implementing a genome-wide transcription-based approach that takes into account genetic and environmental variation to better understand the response of natural populations to stressors. When exposing two different Daphnia pulex genotypes (a cadmium-sensitive and a cadmium-tolerant one) to cadmium, the toxic cyanobacteria Microcystis aeruginosa, and their mixture, we found that observations at the transcriptomic level do not always explain observations at a higher level (growth, reproduction). For example, although cadmium elicited an adverse effect at the organismal level, almost no genes were differentially expressed after cadmium exposure. In addition, we identified oxidative stress and polyunsaturated fatty acid metabolism-related pathways, as well as trypsin and neurexin IV gene-families as candidates for the underlying causes of genotypic differences in tolerance to Microcystis. Furthermore, the whole-genome transcriptomic data of a stressor mixture allowed a better understanding of mixture responses by evaluating interactions between two stressors at the gene-expression level against the independent action baseline model. This approach has indicated that ubiquinone pathway and the MAPK serine-threonine protein kinase and collagens gene-families were enriched with genes showing an interactive effect in expression response to exposure to the mixture of the stressors, while transcription and translation-related pathways and gene-families were mostly related with genotypic differences in interactive responses to this mixture. Collectively, our results indicate that the methods we employed may improve further characterization of the possibilities and limitations of transcriptomics approaches in the adverse outcome pathway framework and in predictions of multistressor effects on natural populations.

Introduction

Standardized ecotoxicology assays are designed to ensure highly reproducible results. However, this consistency, which is achieved by limiting sources of variability typically encountered in natural environments, can restrict their utility. For example, toxicity tests are typically conducted by exposing inbred or clonally derived laboratory populations to a single stressor, which contrasts with the diversity of natural populations and the complexities of their environments.[1−4] This contrast reflects in part methodological limitations of assessing complex, multivariate systems with a limited set of apical end points.[5] Modern genome-wide approaches may offer a possible solution and are increasingly being used to unravel the molecular mechanisms that define responses to environmental stressors and to predict their adverse effects on critical biological pathways.[6−11] However, the application of these variable-rich (i.e., end point-rich) methods to multiple stressor exposures has been limited to a narrow-range of genetic backgrounds.[12−15] In addition, only a few studies have assessed phenotypic responses of different genotypes to stressor mixtures.[16,17] However, these have not focused on environmental genomic approaches designed to understand how gene function is influenced by environment conditions while accounting for the variation that exists within and among natural populations.

As observed by Altenburger et al. (2012),[18] previous studies that assessed the effect of multiple stressors on transcriptomic patterns, failed to compare these patterns against a theoretical baseline derived from the null expectation of a noninteractive response (e.g., independent action or concentration additivity). Yet, this is required to test if one stressor influences the effects of another stressor. Thus if genomic tools are to improve predictions of stressor effects by regulatory agencies, they need to be developed in parallel with the appropriate bioinformatic and statistical approaches required for their interpretation.[19,20]

The present study aims to address these needs through application of genomic resources developed for the waterflea, Daphnia pulex, which includes a well-annotated reference genome.[21]Daphnia sp. are an integral component of freshwater ecosystems that naturalists have employed for centuries as an animal sentinel to gauge the quality of freshwater lakes and ponds.[22] Here we will focus on how we can begin to disentangle the complexities of stressor mixtures faced by natural populations using two distinct Daphnia pulex genotypes, and we do so by simultaneously investigating the phenotypic and transcriptomic response to cadmium, Microcystis and their combined stress. Specifically, this study is, to the best of our knowledge, the first to provide a method for analyzing the interactive effects of mixtures on the transcriptome using the theoretical framework of the independent action (IA) model of joint stressor action.[23]

The two D. pulex genotypes used in the present study differ in their sensitivity to cadmium stress, because of different histories of metal exposures.[24] The genomes of these isolates are shaped by over a century of differential exposure to the selective forces of cadmium and we test the hypothesis that their transcriptomic responses differ upon exposure to cadmium, as well as, upon exposure to stressors that operate via partly similar mechanisms (i.e., Microcystis).[24] Cadmium is an ubiquitous environmental stressor that still imposes risks to some aquatic ecosystems.[25]Microcystis sp. is one of the most common cyanobacteria found in harmful algae blooms.[26] It produces microcystin, a known neurotoxin,[27] and like other cyanobacteria it is predicted to increase in incidence and bloom intensity, even at more Northern latitudes, as a consequence of global climate change.[28]Microcystis has also been observed in cadmium-contaminated lakes,[29,30] which underlines the ecological relevance of this stressor combination. Cadmium and Microcystis also share known mechanisms of toxicity. Both influence the nuclear factor erythroid 2-related factor (Nrf2) oxidative stress pathway[31−37] and both affect the digestive enzymes in Daphnia.[6,38−40]

The present study highlights the utility of implementing an environmental genomics approach that takes into account environmental variation and genetic background of two distinct genotypes as a step forward to a better understanding of the response of natural populations to stress. In summary, we describe and discuss the data by: (1) characterizing the effects of two single stressors on genome-wide transcription and the response of biological pathways and gene-families, (2) isolating the genotype-dependent response to single-stressors, (3) defining interaction effects of mixtures at the transcription-level, and (4) assessing if the nature of mixture interactions varies across genetic background. Finally, we discuss how the application of these four approaches, which are highlighted by characterizing the molecular mechanisms of divergent Daphnia pulex isolates to mixtures of two co-occurring stressors—cadmium and Microcystis—may improve characterization of Adverse Outcomes Pathways and risk-based predictions of stressor effects on natural populations.

Materials and Methods

Experimental Animals

Two Daphnia pulex genotypes were obtained from isoclonal laboratory cultures of the isolate, K10, originating from Kelly lake, Greater Sudbury, Ontario, Canada and the isolate BH14, originating from Basshaunt lake, Dorset, Ontario, Canada. Previous studies showed high tolerance to cadmium in K10 and low tolerance in BH14[24] and as such these genotypes will be referred to as the tolerant (K10) and sensitive (BH14) genotype. Both D. pulex genotypes and the cyanobacterial microcystin-producing Microcystis aeruginosa strain UTEX LB2385 were cultured as described in Asselman et al.[6] Toxin composition analysis of the Microcystis strain indicated the presence of 0.042 to 0.046 mg microcystin·L–1. (6)

Fitness: Reproduction and Growth

A 16 days chronic life-table test was performed to assess the fitness of the genotypes exposed to cadmium, Microcytis aeruginosa and their combination at EC50 effect concentrations that were defined in previous tests across a background of 24 D. pulex genotypes (Shaw, personal observation). The life-table test followed a full-factorial “cube” design with 3 factors, each having 2 levels: genotype (tolerant vs sensitive), cadmium (control vs 0.5 μg Cd·L–1) and M. aeruginosa (0% vs 50% M. aeruginosa). Tests were performed according to OECD guideline 211.[6,41] Cadmium exposure media were spiked with CdCl2 prior to use to a nominal concentration of 0.5 μg Cd·L–1. Regional concentrations (away from point sources) of cadmium in Europe are reported to be between 0.01 and 0.31 μg Cd·L–1.[42] However, local concentrations (close to point sources) can be much higher.[42] Animals in the control and cadmium treatments were fed daily 1.5 mg DW·L–1Ankidostresmus falcatus. Fifty percent of this diet (DW based) was replaced with M. aeruginosa in the Microcystis and combined stressor treatments. Growth was defined as the difference in body length at the start and the end of the experiment. Data were analyzed by means of a 3-way ANOVA followed by Duncan’s posthoc test. Reproduction data did not meet assumptions of normality (Shapiro-Wilk’s test) and homoscedasticity (Levene’s test) and were square root transformed. All tests were performed at a significance level of 95%.

Chemical Analyses and Internal Microcystin Levels

Dissolved cadmium concentration was determined as described in De Coninck et al.[17] Briefly, samples were collected twice a week from freshly prepared and 48 h-old medium (just after media renewal), filtered through a 0.45 μm Acrodisc filter (Sterlitech, Kent, OH) and acidified with 1% (v/v) 14N HNO3 prior to analysis with graphite furnace atomic absorption spectroscopy (SpectrAA-100, Varian, Mulgrave, Australia). Internal microcystin levels of daphnids were determined with the QuantiPlate microcystin kit (Envirologix) following manufacturer’s protocol (see also Supporting Information for more details).

Genome-Wide Transcription, Pathway and Gene-Family Responses

For the microarray experiment, twenty less than 24h old neonates were exposed during 16 days in one liter of medium using the same full-factorial design and under the same experimental conditions as those used in the life-table test. All exposures were performed simultaneously. Sixteen days of exposure allowed accurate assessment of the impact of the stressors on reproduction (see 2.2), given that the daphnids generally produced three clutches of offspring within that period. Furthermore, we believe that this longer period is a better mirror of ecologically relevant settings where daphnids are exposed during longer periods and to lower stressor concentrations. From the 20 exposed animals (i.e., one biological replicate) RNA was extracted, processed and hybridized on the NimbleGen D. pulex 12-plex long-oligonucleotide microarray (GEO:GPL11278)[21] following a full-factorial design with four biological replicates and as described in Asselman et al.[6] with minor modifications (Supporting Information Figure S1). The microarray contained probes to query the expression of the 30 000 validated gene models within the Daphnia genome.[21] Data of the entire microarray experiment was analyzed using linear models implemented in the LIMMA package for R, Bioconductor[6,21] and was submitted to GEO under accession number GSE25843. All statistics and further analyses were performed by defining different contrasts over the linear model fitted to the entire data set. By defining some of these contrasts in an “ANOVA-like” manner (Supporting Information Table S1) we were able to determine main effects of the factors cadmium, genotype and Microcystis and their first and second degree interactions (Supporting Information Figure S2). Other contrasts for further downstream analyses and reaction norm constructing are given in Supporting Information Table S2. The reaction norms provide a more detailed insight in how individual genes within enriched pathways or gene-families respond to the particular stressor or combination of stressors. In particular, where tables and figures identify which pathways and which genes are enriched or differentially expressed, reaction norms visualize how this enrichment or differentially expression differs under different forms of stress or between different genotypes.[43] All contrasts bare a unique number, used for cross-referencing in the manuscript’s text. Pathways and gene-families enriched with significantly expressed genes (i.e., genes of which the contrast value is significantly different from zero) were determined using Fisher’s exact test according to Asselman et al.[6] All tests were corrected for multiple testing using the Benjamini-Hochberg method at a false discovery ratio of 1%. More information can be found in Supporting Information. Sixty percent of the genes on the array are annotated of which the majority comprised those genes associated with known and conserved biological pathways.[21]

Results

Chemical Analyses and Internal Microcystin Levels

Mean dissolved cadmium concentrations were 0.42 and 0.29 μg·L–1 for fresh and 48 h-old medium, respectively and did not differ significantly between replicates or life-table or microarray experiments (F-test, p > 0.05, Supporting Information Table S3). The sensitive D. pulex genotype accumulated significantly more microcystin (80 ± 22 ng·(g tissue)−1) than the tolerant genotype (27 ± 7 ng·(g tissue)−1) (Mann–Whitney U test, p < 0.05).

Single Stressor Responses and Its Variation among Genotypes

Cadmium adversely affected reproduction in the sensitive genotype only, but negatively impacted growth in both genotypes (Figure 1). In contrast, cadmium elicited a weak response on gene expression in both genotypes. Only nine and eleven differentially expressed genes were detected in the sensitive and tolerant genotype, respectively, which was too few to explore their functional significance within biological pathways or gene-families. Only one of the significantly expressed genes was common between both genotypes (Supporting Information Figure S3).

Figure 1

Life-table test results for total reproduction per female (upper panel) and growth, measured as the difference in body length between the start and the end of the exposures (lower panel). Letters denote homogeneous groups based on Duncan’s posthoc test (p < 0.05). Ctrl, control treatment (i.e., absence of stressors); Cd, cadmium exposure; MC, Microcystis exposure; Cd+MC, combined cadmium and Microcystis exposure.

Exposure to Microcystis resulted in significantly reduced reproduction and growth in the sensitive genotype, but not in the tolerant genotype (Figure 1). However, in contrast to cadmium, exposure to Microcystis did result in several differentially expressed genes (contrasts no. 8–9). The gene and pathway responses to Microcystis in the sensitive genotype have previously been described.[6] Five times more differentially expressed genes were detected in the sensitive genotype compared to the tolerant genotype (Supporting Information Figure S4). In total, 471 differentially expressed genes were identified in the tolerant genotype, of which 107 were also detected in the sensitive genotype (Supporting Information Figure S4). These 471 differentially expressed genes were enriched in the steroid hormone biosynthesis and histidine metabolism pathways, none of which were reported in the sensitive genotype (Supporting Information Table S4). Seven gene-families were found to be enriched in the tolerant genotype. Three of these, the trypsin, neurexin IV, and serine/threonine protein kinase families, were also enriched in the sensitive genotype (Supporting Information Table S4).

The above results show that each genotype responded differently to both single stressors. This is also confirmed by significant cadmium × genotype and Microcystis × genotype interaction terms (Figure 1, Supporting Information Table S5). Indeed, we found that 100 and 4676 genes showed a cadmium × genotype and Microcystis × genotype interaction, respectively (contrasts no. 5 and 6; Table 1). Of the 100 genes showing a cadmium × genotype interaction, 51 genes also showed a genotype main effect (Supporting Information Figure S5). Of the 4676 genes that showed a Microcystis × genotype interaction, 1773 also showed a genotype main effect, 408 genes also showed a main Microcystis effect, and 632 also showed both a genotype and Microcystis main effect (Supporting Information Figure S6).

Table 1

Number of Significant Genes in ANOVA-like Main and Interactive Effects (Contrasts No. 1–7 in Supporting Information Table S1; q ≤ 0.01)a

	genotype	Cd	MC	Cd × genotype	MC × genotype	Cd × MC	Cd × MC × genotype
number of sig. genes	8235	7	1351	100	4676	258	430
	27.87%	0.02%	4.57%	0.34%	15.83%	0.87%	1.46%
number of sig. genes with M > 0	4248	0	435	46	2302	32	209
	14.38%	0.00%	1.47%	0.16%	7.79%	0.11%	0.71%
number of sig. genes with M < 0	3987	7	916	54	2374	226	221
	13.49%	0.02%	3.10%	0.18%	8.04%	0.77%	0.75%

The percentage shows the ratio of affected genes compared to all putative genes on the array (29 546 genes). Sig., significant; Cd, cadmium; MC, Microcystis aeruginosa; M, log2(contrast) value. Biological interpretation is dependent on the contrast (i.e., ANOVA-like effects) as described in Supporting Information Table S1. For instance, an M-value > 0 for the main genotype effect indicates that a gene has a higher expression level in the sensitive genotype compared to the tolerant one and an M-value < 0 for the Microcystis × genotype interaction indicates a higher differential expression under Microcystis exposure (relative to the control) in the sensitive genotype compared to the tolerant genotype under Microcystis exposure (relative to the control).

No pathways enriched with genes showing a cadmium × genotype interaction were detected. In contrast, seven pathways were enriched with genes showing a Microcystis × genotype interaction and these were related to oxidative stress, energy metabolism, and lipid metabolism (Figure 2). Three pathways (glutathione, starch- and sucrose, and linoleic acid metabolism) were also significantly enriched with genes showing a genotype main effect (Figure 2). Gene-expression in the pathways was mostly affected by Microcystis exposure in the sensitive genotype while gene-expression in the tolerant genotype was almost not affected (Table 2, Supporting Information Figure S9).

Figure 2

Overview of median log2(contrast) values of significantly enriched KEGG-defined pathways (panel A) and functionally KOG-annotated gene-families (panel B). Pathway and gene-family names are indicated on the right of the color matrix, ANOVA-like effects (i.e., main effects and interactions as defined by contrasts no. 1–7) are given below the matrix. A gray square in a given row-column combination indicates that the pathway or gene-family in that row was not significantly enriched with genes showing the ANOVA-like effect in that column. Red colors indicate a median log2(contrast) value significantly greater than 0 (or contrast >1) and green colors indicate a median log2(contrast) value significantly lower than 0 (or contrast <1). Biological meaning of the log2(contrast) value depends on the definition of the contrasts (Supporting Information Table S1). For instance, for the Microcystis × genotype interaction, a red color indicates that gene-expression under Microcystis exposure in the cadmium sensitive genotype is upregulated relative to control conditions compared to gene-expression under Microcystis exposure (relative to control conditions) in the cadmium tolerant genotype. A green color indicates the opposite: downregulation of gene-expression under Microcystis stress in the cadmium sensitive genotype compared to gene-expression under Microcystis stress in the cadmium tolerant genotype. Numbers to the left indicate clusters of pathways of gene-families based on their functions: 1, oxidative stress related pathways; 2, energy metabolism related pathways; 3, lipid and poly unsaturated fatty acids metabolism related pathways; 4, other pathways; 5, digestion related gene-families; 6, oxidative stress related gene-families; 7, signal transduction related gene-families; 8, transporter related gene-families; 9, transcription and translation related gene-families; 10, other gene-families. The two upper rows give a visual representation of the ANOVA results of reproduction and growth data: blue squares indicate significant ANOVA terms, gray squares nonsignificant terms.

Table 2

Between-Genotype Comparison of Significant Genes Following Exposure to Microcystis (MC) (i.e., MC versus Control for Both Genotype Separately, Contrasts No. 8–9 in Supporting Information Table S1) in Significantly Enriched KEGG-Defined Pathways and Functionally KOG-Annotated Gene-Families Detected in the Microcystis × Genotype Contrast (Contrast No. 6)a

	MC × genotype		tolerant genotype MC vs control			sensitive genotype MC vs control
	q-value pathway	no. of genes in pathway	q-value pathway	q < 0.01 M > 0	q < 0.01 M < 0	q-value pathway	q < 0.01 M > 0	q < 0.01 M < 0
Pathways
oxidative phosphorylation	<0.0001	148	0.5084	0	0	<0.0001	46	4
mitochondrial dysfunction	<0.0001	107	1	0	1	<0.0001	39	2
glutathione metabolism	<0.0001	172	1	3	0	0.0002	9	31
ribosome	0.0023	351	1	5	0	<0.0001	72	19
starch and sucrose metabolism	<0.0001	127	0.0369	8	0	0.0010	8	28
linoleic acid metabolism	0.0006	24	0.2935	2	0	0.0003	2	10
steroid hormone biosynthesis	0.0013	29	0.0266	3	0	0.0010	3	9
Gene-families
trypsin	<0.0001	255	<0.0001	14	10	0.2676	20	33
zinc carboxypeptidase	0.0005	39	0.0528	6	0	0.0006	9	9
NADH:ubiquinone oxidoreductase	<0.0001	30	1	0	0	<0.0001	21	0
neurexin IV	<0.0001	50	<0.0001	10	0	0.0006	3	18
vacuolar H+-ATPase	0.0029	10	1	0	0	0.7432	4	0
40S ribosomal protein	<0.0001	36	1	0	0	<0.0001	24	0
60S ribosomal protein	0.0029	51	1	0	0	<0.0001	27	1
RNA polymerase II	<0.0001	121	1	6	0	0.0829	5	26
apoptosis-inducing factor	<0.0001	26	1	0	0	<0.0001	0	18

q, Benjamini-Hochberg corrected p-value; M, log2(contrast) value.

Similar to the pathway analysis, no enriched gene-families were detected under cadmium exposure (contrast no. 5). In contrast, nine were detected under Microcystis exposure (contrast no. 6; Figure 2) and these represent different biological processes, such as digestion, oxidative stress response, translation/transcription and signaling (Figure 2). Only one of these nine families, that is, trypsins, was also enriched in both the main genotype and main Microcystis contrast. While most gene-families and pathways contained almost no differentially expressed genes in the tolerant genotype (Table 2), two gene-families, that is, trypsins and neurexins, did have a number of differentially expressed genes in the tolerant genotype. Yet, an almost completely different set of significant trypsin isoforms was detected in both genotypes (Figure 3). This was less the case for neurexins (Figure 3).

Figure 3

Reaction norms for expression of genes following Microcystis exposure (presence vs absence) in the cadmium tolerant (dotted line) and sensitive genotype (solid line) for (A) the trypsins gene-family and (B) the neurexins IV gene-family. Numbers close to each dotted line indicate the number of genes showing the corresponding particular expression pattern in the tolerant genotype. The sum of those numbers in each plot gives the total number of genes showing the specific expression pattern of the sensitive genotype (left: no response; right: upward response; right: downward response). For instance, in panel A, center plot, we can see that while upregulation of gene-expression in 20 genes is detected after exposure to Microcystis in the sensitive genotype, one gene in the tolerant genotype was upregulated after exposure to a level that equaled the level in the sensitive genotype. Two genes in the tolerant genotype were not differentially expressed and had constitutively already the same level of expression as the upregulated genes in the sensitive genotype in presence of Microcystis. Similarly, in panel A, right plot, nine genes in the tolerant genotype had constitutively the same expression level in presence of Microcystis as those genes that were downregulated in the sensitive genotype. In panel A, left plot, finally, the constitutive expression of 33 genes in the tolerant genotype was higher than the constitutive expression of genes in the sensitive genotype. For 19 genes the opposite was true. Contrasts used for constructing the reaction norm plots are given in Supporting Information Table S2.

Cadmium × Microcystis Mixture Interactions

When assessing mixture interaction effects relative to the Independent Action (IA) reference model, no deviations from additivity were detected for reproduction and growth (Supporting Information Table S5). This was also reflected in a limited number of genes showing a cadmium × Microcystis interaction effect (only 0.9% of the genes on the array; contrast no. 4; Table 1). Of these 258 genes, 222 genes uniquely showed the interaction effect, while 35 genes also showed a main Microcystis effect, and only two genes also showed a cadmium main effect (Supporting Information Figure S7). Only one enriched pathway (ubiquinone biosynthesis) and few gene-families (collagens and MAPK related serine-threonine protein kinases) were detected (Figure 2). Similar to pathways enriched with genes showing a Microcystis × genotype interaction, the expression of the genes in these few pathways and gene-families was studied more in depth by means of reaction norm plots (Supporting Information Figure S10).

Genotypic Variation in Mixture Interactions

At the organismal level, no genotypic differences in mixture response were detected (i.e., no genotype × cadmium × Microcystis interaction; Supporting Information Table S1). However, a limited number of genes, that is, 1.5% of the putative genes on the array, showed a 3-way interaction (contrast no. 7; Table 1). Of these 430 genes, 188 and 77 genes also showed a main genotype response or main Microcystis response, respectively. Forty-three of these genes showed both a main genotype and Microcystis effect (Supporting Information Figure S8). None of the genes showing a three-way interaction overlapped with genes showing a cadmium main response (Supporting Information Figure S8). One pathway, the ribosome pathway, was enriched with members of this gene set (Figure 2). Two gene-families were enriched with members of these genes showing a three-way interaction. Both were related to transcription and translation processes and did not overlap with gene families enriched with genes showing one of the main effects (Figure 2). Reaction norm plots showed an influence of the genotype on the significant genes in the cadmium × Microcystis interaction (Supporting Information Figures S11–S13).

Discussion

The present study shows the possibilities and limitations of implementing an environmental genomics approach that takes into account genetic and environmental variation to better understand the response of natural populations to stress. The methods we present allow for characterizing the effects of multiple stressors on genome-wide transcription and the response of biological pathways and gene-families, for isolating the genotype-dependent response to single-stressors, for defining interaction effects of mixtures at the transcription-level, and for assessing if the nature of mixture interactions varies across genetic background. Our results for each of these effects are discussed in terms of application to Adverse Outcome Pathways (AOP) or risk-based predictions of stressors effects in natural populations. AOPs represent a sequence of events that begins with a molecular initiating event, spans multiple levels of biological organization and ends with an adverse or toxic outcome at the whole-organism level.[10,44]

The application of transcriptomics data to an AOP approach has predominantly been used to interrogate acute exposures to single stressors,[9−11,45] but in the current study the transcriptome response to single stressors did not always correlate with higher-level effects (i.e., growth and reproduction). Molecular approaches, such as microarrays, are typically used to study perturbations in molecular toxicity pathways.[10] Use of this information within an AOP framework should ultimately lead to predictions of effects at higher biological levels, which are currently still more useful in ecological risk assessment.[10] However, our study demonstrates that this response relationship is not always clear-cut. For example, cadmium elicited an adverse effect on both reproduction and growth in one genotype, yet a negligible effect of the stressor on gene-expression was observed (only 9 differentially expressed genes). Although our analysis was quite stringent by applying a significance cutoff of p < 0.01, this conclusion was still maintained when a less stringent cutoff of p < 0.05 was applied (Table S6 in Supportive Information). Our observed results differ from the literature where robust gene-expression changes following acute exposure (i.e., 24 to 96 h) to cadmium have been reported,[46−50] and these genes could be coupled with some known mechanisms of action of and response pathways to cadmium.[46−50] However, our studies measured gene transcription after chronic (i.e., several days) exposures to sublethal cadmium concentrations to better reflect the ecological realities of cadmium exposures in natural environments.[25] More complex regulatory pathways such as those resulting in acclimation have been reported following these longer-duration exposures.[50−53] Specific detoxification mechanisms, such as metallothionein (MT) proteins, have an important role in the defense against cadmium toxicity.[47,50] However, these MTs, and other genes, are often expressed in a time-dependent manner.[50,52] Asselman et al.,[50] for instance, showed a significant response of the four MT genes in D. pulex exposed to cadmium after two to eight days, but not after 16 days of exposure. This may be due to the long turnover time of MTs causing a discrepancy between gene expression and protein levels. Such dynamic processes could explain why transcriptomic responses may not always correlate with higher-level responses for all exposure levels or durations.

In contrast to cadmium, transcriptomic patterns in response to Microcystis stress did show a complexity that mirrored results at the organismal level. In this case, the gene-expression patterns revealed that reduction in reproductivity and growth was a result of several different mechanisms, such as the lack of essential nutrients (e.g., essential fatty acids or lipids), or the presence of feeding deterrents and toxins, or oxidative stress. The mechanisms revealed by the transcriptional response agree with previous results at both the transcriptomic[6] and organismal level.[54−56] Thus, combining the results obtained for the cadmium and the Microcystis exposure, we show that the ability of transcriptomic profiles to predict higher-level effects is dependent on the stressor and exposure conditions under consideration.

Another potential benefit of gene-expression data is to disentangle the variability of the environments in which natural populations reside. Indeed, using conventional toxicity assays with limited apical end points, it is difficult to assess the complexities of natural environments that include exposure to stressor mixtures beyond their directional influences on toxicity (i.e., greater or less than additive). For example, gene-expression profiles provide a more detailed biological understanding of the mechanisms by which mixture components interact. These data allow interactions to be assessed with greater knowledge than the phenotypes of, for instance, reproduction or growth by identifying pathways and gene-families that drive interactions.

The transcriptomic response to multiple stressors has been studied before in Daphnia,[12−15] but there is a clear difference between all these studies and our present study in the way these responses are analyzed and interpreted. The previous studies have compared the recorded transcriptomic response to mixtures against responses to individual stressors, but have not formally evaluated interactive effects between the stressors.[18] Yet, interactions between multiple stressors, regardless of whether they occur at the gene-expression, physiological, individual or population level can only be judged against an appropriate baseline of additivity.[18,22] By joining well-established mixture toxicity baseline models with gene-expression data the current limitation of providing mainly anecdotal evidence on mixture effects at the transcriptomic level may be overcome.[18] The present study is to our knowledge the first to investigate mixture interactions with transcriptome data using the Independent Action (IA) reference model, which is typically used to set the baseline for mixture interactions between unrelated stressors.[17] The interaction responses of expressed genes identified in our study provided insight into the biological meaning of mixture interactions. For example, the ubiquinone biosynthesis pathway was the only pathway significantly enriched with expressed genes showing a cadmium × Microcystis interaction effect. Ubiquinone functions as an electron-carrier in the electron transport chain in mitochondria,[57] and because its electrons are only loosely bound when ubiquinone is in a reduced state, it has been suggested to play an important role as an antioxidant in the oxidative stress response in Daphnia.[7] Gene-families enriched with expressed genes showing a cadmium × Microcystis interaction effect were MAPK related serine-threonine protein kinases and collagens of type IV and XIII. Mitogen-activated protein kinases (MAPK) are serine-threonine protein kinases that are expressed in all eukaryotic cells and that are activated by diverse stimuli ranging from cytokines, growth factors, neurotransmitters, hormones, cellular stress, and cell adherence.[58] MAPK pathways have been previously reported to be involved in nutrient deprivation and responses to stress stimuli such as osmolarity changes.[59] Nutrient deprivation is most likely, as discussed above, a consequence of the lack of essential fatty acids in Microcystis, while osmolarity changes may be a consequence of oxidative stress.[59,60] Collagens generally support tissue around the organs and under the epidermis in insects.[61] Several endogenous immunostimulatory peptides in insects that may serve as danger and alarm signals are derived from collagen type IV.[62] Although no interaction between cadmium and Microcystis was observed with apical end points, these studies demonstrate the potential use of gene-expression patterns to better understand interactions occurring in stressor mixtures.

In addition to mixture interactions, transcriptome data can also dissect variation introduced by the diversity residing within natural populations. In other words, these methods have the potential to disentangle the sources of variation that give rise to the most central tenet in ecotoxicology, the dose–response relationship.[63] In the current study, differences between the genotypes are reflected both at the organismal and the gene expression level. Pathways and gene-families enriched with expressed genes that show a Microcystis × genotype interaction effect are not differentially regulated following Microcystis exposure in the cadmium tolerant genotype, but are in the sensitive genotype. This observation could suggest that genes involved in tolerance mechanisms to Microcystis stress already have the basal expression level needed to cope with the stress in the tolerant genotype. The reaction norm approach provided in the current study provides a method for visualizing this hypothesis. Indeed, this hypothesis is only valid if genes in the sensitive genotype following Microcystis exposure reach the same expression level as that of the constitutively expressed genes in the tolerant genotype. Based on the reaction norms, this was clearly not case for the majority of the genes (Figure 3; Supporting Information Figure S9). Other mechanisms must thus be at play, for example those preventing that Microcystis can exert its toxicity in the cell, to explain why the tolerant genotype shows little to no effects of Microcystis exposure both at the organismal and the transcriptomic level. First, this tolerance could be acquired by limiting internal, cellular exposure to Microcystis. Measurements of internal microcystin concentrations, which were significantly lower in the tolerant genotype, provide added support for such a mechanism. Second, the observation of a Microcystis × genotype interaction for genes encoding trypsins and neurexins, suggests a role of these gene-families in the genotypic difference in tolerance to Microcystis stress. A previous study has reported both up- and downregulated genes encoding for trypsins in D. pulex following Microcystis exposure.[6] These authors argued that this up- and downregulation is in accordance with the observations by Agrawal et al.[64] and Schwarzenberger et al.[40] who noted differential sensitivity of trypsins to Microcystis, both at the RNA and the protein level. Asselman et al.[6] speculated that trypsin isoforms, which are more sensitive to Microcystis were downregulated while less sensitive isoforms were upregulated. Our data show a very different set of differentially expressed trypsin genes between the sensitive and the tolerant genotype and therefore identify those isoforms that can be considered candidate isoforms that may correlate with tolerance differences and that could be an interesting subject of more detailed future studies (Figure 3). Similar observations were made for neurexins. However, the function of neurexin in Daphnia and its role in response to Microcystis stress remains unclear. Interestingly, the trypsin gene-family is also enriched with genes showing a main Microcystis effect indicating that some members of this gene family respond similarly to Microcystis stress across both genotypes, underscoring the potential power of gene-expression data to accurately disentangle genotypic differences in tolerance. Understanding the mechanisms that account for tolerance difference highlight the fact that these transcriptomic methods can provide ecologically important insights.

Clear differences between both genotypes existed in their response to Microcystis stress. First, the cadmium tolerant genotype proved to be also more tolerant to Microcystis at the organismal level, while the cadmium sensitive genotype was more sensitive to Microcystis stress (i.e., decreased reproduction after Microcystis exposure compared to a control). Second, the cadmium tolerant genotype showed, compared to the sensitive genotype, only few differentially expressed genes in the pathways and gene-families that were enriched with genes showing a Microcystis × genotype interaction effect. These observations strongly suggest that the cadmium tolerant genotype, originating from a historically, metal polluted environment, is better armed to deal with Microcystis stress. Moreover, the common toxicity mechanisms between cadmium and Microcystis that we hypothesized to be the armory for interactions between these two stressors, indeed provided the arms that enable the cadmium tolerant genotype to also better cope with Microcystis stress. Pathways and gene-families related to oxidative stress, food quality and the digestive system were enriched in significantly genes that segregated the two genotypes. Such pleiotropy, when a gene or pathway affects more than one trait, has been previously described as a mechanism that can provide “cross-tolerances” between different stressors.[65]

Our experimental design also allows us to interpret part of the complexity of stress responses in a natural environment, that is, the genotypic-dependent response to stressor mixtures. Only a limited number of pathways and gene-families were significantly enriched with genes showing a cadmium × Microcystis × genotype interaction effect. These included the ribosome pathway and the histone 2A and 60S ribosomal protein gene-families. Interestingly, all of these pathways and gene-families are related to transcription and translation processes, which suggest a general stress response likely to protect proteins.[6,66,67] In addition, our reaction norm approach revealed in more detail genotypic differences in mixture interaction response: cadmium altered the gene-expression of genes responding to Microcystis stress more pronouncedly in the sensitive genotype than in the tolerant genotype.

Overall, our results demonstrate the utility and also highlight some limitations of applying environmental genomics approaches to dissect the complexity of natural environments and the diversity of natural populations. We were able to interrogate the effects of mixture components and genotypes both independently and in combination, and identify interactions responses among gene families and pathways that ultimately contributed to tolerance differences between individuals. However, these approaches were less successful at linking gene-expression results under exposures to single compounds to organismal level responses. Future studies should examine more single stressors and also more (and also more complicated) mixtures across a range of genotypes to better define the limits of complexity that these approaches can dissect. This will help to further reveal the possibilities and limitations of these techniques to dissect more complex environmentally relevant exposure scenarios in an AOP framework.