Transcriptomic, biochemical and individual markers in transplanted Daphnia magna to characterize impacts in the field.

Daphnia magna individuals were transplanted across 12 sites from three Spanish river basins (Llobregat, Ebro, Jucar) showing different sources of pollution. Gene transcription, feeding and biochemical responses in the field were assessed and compared with those obtained in re-constituted water treatments spiked with organic eluates obtained from water samples collected at the same locations and sampling periods. Up to 166 trace contaminants were detected in water and classified by their mode of action into 45 groups that included metals, pharmaceuticals, pesticides, illicit drugs, and other industrial compounds. Physicochemical water parameters differentiated the three river basins with Llobregat having the highest levels of conductivity, metals and pharmaceuticals, followed by Ebro, whereas the Jucar river had the greatest levels of illicit drugs. D. magna grazing rates and cholinesterase activity responded similarly than the diversity of riparian benthic communities. Transcription patterns of 13 different genes encoding for general stress, metabolism and energy processes, molting and xenobiotic transporters corroborate phenotypic responses differentiated sites within and across river basins. Principal Component Analysis and Partial Least Square Projections to Latent Structures regression analyses indicated that measured in situ responses of most genes and biomarkers and that of benthic macroinvertebrate diversity indexes were affected by distinct environmental factors. Conductivity, suspended solids and fungicides were negatively related with the diversity of macroinvertebrates cholinesterase, and feeding responses. Gene transcripts of heat shock protein and metallothionein were positively related with 11 classes of organic contaminants and 6 metals. Gene transcripts related with signaling paths of molting and reproduction, sugar, protein and xenobiotic metabolism responded similarly in field and lab exposures and were related with high residue concentrations of analgesics, diuretics, psychiatric drugs, β blockers, illicit drugs, trizoles, bisphenol A, caffeine and pesticides. These results indicate that application of omic technologies in the field is a promising subject in water management.


Introduction
Identifying indicators of adverse change in ecological systems which can diagnose causal agents is a major challenge in environmental risk assessment (Baird and Burton, 2001).
Traditionally, biomonitoring of fresh waters has been based on measures of community structure, focussing on biodiversity metrics (Rosenberg and Resh, 1993). Such diagnostic measures have been widely used to establish the ecological quality of water (Munné and Prat, 2009). Nevertheless, biodiversity metrics are based on the occurrence of species and hence they are mostly sensitive to dramatic changes such as reductions of individuals within species or even to species extinction. These features make biodiversity metrics unable to detect subtle changes of individual physiological responses. Furthermore, in many cases biodiversity metrics respond to other environmental factors than to trace contaminants (Baird and Burton, 2001).
The development of new bioassays with caged single species has allowed determining pollutant effects in situ. Key advantages of in situ bioassays over whole effluent toxicity tests and biological surveys of invertebrate communities include: a greater relevance to the natural situation, especially with respect to the contamination scenario; and their ability to detect effects more rapidly (hours to days) than resulting changes in community structure (months to years) measured during macroinvertebrate sampling (Maltby et al., 2002) . More specifically, a set of in situ and cost effective bioassays based on feeding and biochemical responses of invertebrate species have permitted detecting lethal and sublethal responses that are biologically linked with key ecological processes such as detritus processing and algal grazing rates, and of specific toxicological mechanisms (Barata et al., 2007;Burton Jr et al., 2005;Damásio et al., 2010;Maltby et al., 2002;Maltby et al., 2000;Schulz and Liess, 1999).
Nevertheless such developments were still limited to few responses. The use of omic technologies may offer the possibility to extent those responses to many genes within individuals.
The water flea Daphnia magna is possibly the invertebrate species most used in toxicology and experimental ecology and together with its close relative D. pulex is used as a model for environmental genomics research . D. pulex genome has been fully sequenced and about 50% of its genome is annotated (Colbourne et al., 2011), thus that of its close relative D. magna, despite of being incomplete, may benefit from the former. Indeed several studies have identified gene markers within D. magna genome that respond to specific pollutants. These genes include general stress genes such as the heat shock protein 70 (HSP70) and metallothioneins (MT2), whose transcripts are induced by several metals (Ho, 2008;Poynton et al., 2007); genes involved in metabolic pathways, i.e. the metabolism of the sugars (UDP-glucose pyrophosphorylase, UGP), lipids (thiolase, THIO ), amino acids(fumaril acetoacetase, FAA) and the Krebs cycle ( isocytrate dehydrogenase (IDH), aconitase (ACON), (Campos et al., 2013)); specific genes encoding key processes of growth and reproduction (vitellogenein, VTG; ecdysteroid receptor, EcR; retinoid X receptor, RXR and molting inhibition hormone, MIH) (Montagné et al., 2010;Tokishita et al., 2006;Wang and LeBlanc, 2009;Wang et al., 2007)); and transporter genes from the multixenobiotic resistance mechanisms such as the P-glycoprotein (Pgp) and multidrug resistance protein 4 (MRP4) (Campos et al., 2014).
The Mediterranean basin is one of the world's regions most vulnerable to global change (Barceló and Sabater, 2010) and one of the "hot spots" for ongoing problems in water availability (Giorgi and Lionello, 2008). Here we tested the feasibility of using mRNA responses of D. magna genes in detecting and identify environmental stressors in the field having detrimental effects in river biota. To do that, up to 19 responses in D. magna individuals transplanted in the field were used to evaluate the effects of up to 167 trace contaminants and five general physicochemical parameters in 12 sites belonging to three distinct Mediterranean rivers. The biodiversity of benthic macroinvertebrates were also considered to allow comparison of D. magna responses with those of the whole community.
D. magna responses included that of post-exposure feeding rates, five enzymatic biomarkers (cholinesterase, carboxylesterease, lactate dehydrogenase, catalasa, glutathione S transferase) and 13 genes (HSP70, MT2, ACON, IDH, UGP, FAA, THIO, VTG, EcR, MIH, RXR, PGP, MRP4). Secondly, the response of the studied genes together with that of post-exposure feeding rates was evaluated in D. magna individuals exposed in the lab to organic extracts of water samples collected at the studied sites. This allowed to compare the robustness and repeatability of the studied gene responses in detecting effects of organic contaminant residues. This study focused on three river basins Llobregat, Ebro and Jucar rivers (NE Spain).

Study sites
Like most Mediterranean systems, Llobregat, Ebro and Jucar river basins natural resources have been greatly affected by human activities such as agriculture, urbanization, salinization by mining activities and an intensive water use for human consumption., which together have severely deteriorated the ecological status of the main rivers and tributaries since 1970s (Belenguer et al., 2014;Damásio et al., 2011;Damásio et al., 2010;Damásio et al., 2008;De Castro-Català et al., 2013;Fàbrega et al., 2013). The Llobregat River (northeast Spain) is 156.5 km in length, covers a catchment area of approximately 4,948 km 2 , and its watershed is heavily populated (3,089,465 inhabitants in 1999). The Llobregat River is a paradigm of overexploited Mediterranean river with nearly 30% of its annual discharge used for drinking water. Moreover, the Llobregat River receives extensive urban and industrial wastewater discharges (137,000,000 m 3 /year; 92% comes from the wastewater treatment plants) that cannot be diluted by its natural flow (0.68-6.5 m 3 /s basal flow) (Muñoz et al., 2009). The Ebro river is one of the largest river basins in Spain, 928 km in length and with a drainage basin of 85550 km 2 and around 2 800 000 inhabitants living in the area. The most relevant economic activity in the region is basically agriculture (vineyards, cereals, fruit, corn, horticulture and rice production), but there are also some highly industrialized regions, mainly located in the northern-central part. (Silva et al., 2011). The Júcar River is 497.5 km long drainage area of 21 600 km2 and its mean annual flow is 10 m 3 /s; it flows Eastern Spain, under a typical Mediterranean climate (Belenguer et al., 2014) .
Deployment sites comprised four points along the Llobregat, Ebro and Jucar river systems, respectively (Fig 1). The study stations were chosen as being examples of the different characteristics of the basins.The Llobregat stations were L3 (Pont de Vilomara), L4 (Castellbell i el Vilar), L5 Abrera) and L7 (Sant Joan Despí). L3 site was located at the midsection of the river, without much human impact. L4 was located after the junction with the Cardener river, characterized by its high conductivity and pollutants coming from Manresa sewage treatment plant (STP). L5 and L7 were located at the end of the mid and lower section of the river, downstream the STP of Monistrol de Montserrat and Barcelona, respectively.
The Ebro study sites were all located in the upper-mid section of the river, an area with high water usage for agricultural purposes but also with an important human population settlement.
E2 was located on the upper course at Miranda de Ebro. It is the first major city in the main stream. There is a large industrial area and the STP is less efficient than desired so the concentration of organic compounds in the river is high. E3, E4 and E5 sampling site received the influence of the wine field lixiviates from STP of Haro, Logroño and Tudela cities, respectively.
Jucar Sampling stations were distributed along the full course of the river. The J2 site was located inside the city of Cuenca. Water at J4 site presented a strong eutrophic process and presence of filamentous algae along the reach was observed. J5 was just downstream of a small Dam and J6 downstream the bypass of the irrigation channels that are used in the agricultural explorations further down.

Environmental measurements
A set of environmental variables were measured on each deployment occasion. Water physicochemical parameters including temperature (T; o C), pH, conductivity (S/cm), dissolved oxygen (O2, mg/l) and suspended solids (SS, mg/l) were obtained following Damásio et al. (2008) procedures. Briefly, T, pH, conductivity and O2 were measured in situ by using a WTW Multi 340i handheld meter, whereas total suspended solids were measured in the lab following ASTM Standard Methods (APHA- AWWA-WEF, 1995). Residue levels of eight metals in water and that of up to 158 organic contaminant residues were analysed.
Organic contaminant residues were classify into 37 functional groups according to their mode of action or/and chemical structure. Further details of the established groups are in Table S1 (Supplementary information). Surface water samples (12) were taken from the studied sites.
Duplicate water samples were collected in the middle of the current river with 2.5 L amber glass bottles. Within 48 h, water samples were vacuum filtered through 0.45 μm glass fiber filters and aliquoted into 1L. The aliquots were extracted using OASIS HLB SPE cartridge (200 mg sorbent/6 mL cartridge, Waters) and eluted with different solvents according to the family of compounds that were to be analysed. Ilicit drugs were eluted with a gradient acetonitrile:water and analysed by isotope dilution on-line solid phase extraction-liquid chromatography-tandem mass spectrometry (SPE-LC-MS/MS) following the method described in (Postigo et al., 2008). Pharmaceuticals were eluted with methanol and the concentrations of the 73 compounds were determined using a multi-residue analytical method based on LC-MS/MS (Osorio et al., 2012). Most relevant environmental endocrine disruptors compounds (EDCs) and compounds suspected to be EDCs such as natural and synthetic estrogens and their conjugates, antimicrobials, parabens, bisphenol A, alkylphenolic compounds, benzotriazoles, and organophosphorus flame retardants were analysed using a fully automated approach in which water samples were directly injected into the chromatographic system and the target compounds were concentrated into the loading column.

Biological condition
Invertebrate samples from sediment were randomly collected with a corer (24 cm 2 area, 5 replicates per site). Samples were sieved through a 500 µm mesh to separate invertebrates and fixed with 4% formaldehyde. The invertebrates were identified at species level and used to calculate species richness (S) and Shannon diversity index (Shannon and Weaver, 1963).
More information is in

Field exposures
In situ D. magna deployments were conducted as indicated by using the same test chambers and procedures of (Mc William and Baird, 2002) with only minor modifications that included 10 test chambers to allow collection of animals for gene and biomarker determination and to increase the number of replicates for post-exposure feeding rate measurements. Chambers were constructed from clear polyvinyl chloride cylindrical piping (13 cm long, 5 cm external diameter). Each chamber had two rectangular windows (7 x 3.5 cm) cut into either side of the cage, covered with 150-m nylon mesh. Pipe ends were sealed with polypropylene caps.
Groups of 4-5 chambers were placed inside a 13-mm 2 wire-mesh cylinder that was positioned in the stream perpendicular to flow. In each deployment, a lab control treatment with animals maintained in the lab and never exposed to the field was also included as a surrogate control.
Deployments were conducted in 2011, on 15-16 th September in Llobregat, 21-22nd September in Ebro and 5-6th October in Jucar. Within each period deployments were conducted simultaneously in four locations that always included at least a low polluted site. Briefly the procedure for the in situ bioassays was as follows. Juveniles were transported to field sites in groups of 10 in 175 glass jars filled with American Society for Testing Materials (ASTM) hard water (Mc William and Baird, 2002). At each site 10 chambers, each containing 10-20 individuals were placed inside a 13-mm2 wire-mesh cylinder that was positioned in the stream perpendicular to flow.

Lab exposures to reconstituted water
A second exposure experiment was set in the lab to compare feeding and gene responses with those of the field and to test if responses observed in the field were related to organic contaminant pollution only. For that purpose we use the same water samples and extraction procedures as those reported for chemical analyses but without the addition of standards.
Final eluates were then combined, evaporated dryness with N2 and re-suspended in 0.3 mL of acetone. Lab exposures were then conducted by exposing 50 animals in 2L ASTM hard water dosed with 0.2 mL of the obtained acetone extract eluates. Exposures were conducted in 3L glass bottles gently shaken in an orbital incubator at 20 ºC under darkness and lasted 24 h.
Controls that were processed similarly as field samples but using ASTM water only were included in each trial. Each river was assayed simultaneously. After exposures animals were used for post-exposure feeding and gene responses as described below.

Post exposure responses
After 24 h, animals were retrieved from chambers or lab exposure media. Twenty five surviving animals of five chambers or collected from lab exposures were used to determine post exposure feeding rates and the remaining individuals were pooled in Eppendorf's in groups of five, immediately frozen in liquid N2 and kept at -80 o C until further gene analysis.
From each of the remaining 5 chambers of field exposures, surviving animals were also pooled, frozen and used for biomarker determination.
Shortly after exposure (within 1h) five surviving juveniles were placed into 60 ml screwcapped glass jars containing 50 ml of ASTM hard water, with Chlorella vulgaris (Beijerink, strain CCAP C211/12) at a concentration of 5 x 10 5 cells/ml, and allowed to feed for 4h (Mc William and Baird, 2002). Three jars containing no animals were used to establish initial algal densities. Biomarker, gene and post-exposure feeding rates were also measured in animals maintained in the lab during the duration of the deployments and transported to the field sites to include a surrogate lab control. Post-feeding experiments were conducted in darkness to avoid algal growth and under constant temperature conditions (20  2ºC) provided by a thermostatised chamber. Individual feeding rates (cells animal -1 h -1 ) were determined as the change in cell density during 4h according to the method given by (Mc William and Baird, 2002). Cell density was estimated from absorbance measurements at  = 650 nm in a dualbeam spectrophotometer (Uvikon 941) using standard calibration curves based on at least 20 data points, with an r 2 > 0.98.

Enzyme assays
Juveniles were homogenized at 4 o C in 1: 4 wet wt./buffer volume ratio in 100 mM phosphate buffer, pH 7.4 containing 100 mM KCl and 1 mM EDTA. Homogenates were centrifuged at 10 000 g for 10 min and the supernatants were immediately used as enzyme sources.
Biochemical measurements were carried out on Uvikon 941 Plus dual-beam and Spectra-max Plus microplate reader spectrophotometers. Assays were run at least in duplicate. AChE was determined by a modification of the Ellman method adapted to mircroplate (Barata et al., 2004). Acetylcholinesterase activity was measured in the presence of 1 mM acetylthiocholine and 0.1 mM 5,5' dithiobis-2-dinitrobenzoic acid (DTNB), and the increase of absorbance was measured at 405 nm. Catalase activity was measured by the decrease in absorbance at 240 nm due to H2O2 consumption (extinction coefficient 40 M -1 cm -1 ) according to (Aebi, 1974) . The reaction volume was 1 mL and contained 50 mM phosphate buffer, pH 6.5, 50 mM H2O2 (Ni et al., 1990). Glutathione S transferase (GST) activity towards 1chloro-2, 4dinitrobenzene (CDNB) was measured as described by (Habig et al., 1974). The reaction mixture contained 100 mM phosphate buffer (pH 7.5), 1 mM CDNB and 1 mM of reduced glutathione. The formation of S -2, 4 dinitro phenyl glutathione conjugate was evaluated by monitoring the increase in absorbance at 340 nm. CbE activity was measured by the UV method of in the presence of 0.25 mM α-naphtyl acetate, and the formation of naphthol monitored by the increase in absorbance at 235 nm (Barata et al., 2004) . Lactate dehydrogenase (LDH) activity was determined according to (Diamantino et al., 2001) . Proteins were measured by the method of (Bradford, 1976) using serum albumin as standard.

Gene responses
Total RNA was isolated from the samples using Trizol reagent (Invitrogen, Carlsbad, USA) following the manufacturer's protocol and quantified in a NanoDrop D-1000 Spectrophotometer (NanoDrop Technologies, Delaware, DE). RNA quality was checked in an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara CA). Quantities of 1µg were retro-transcribed to cDNA using First Strand cDNA Synthesis Kit Roche®77 (Germany) and stored at -20ºC. Thirteen genes were selected for representation of different pathways/gene families: HSP70, MT2, ACON, IDH, UGP, FAA, THIO, VTG, EcR, MIH, RXR, PGP, MRP4. The gene glyceraldehyde 3-phosphate dehydrogenase-G3PDH was used as an internal control. For each of these genes primers were designed with Primer Quest (IDT Technologies, Coralville, IA) and are listed in the Supporting Information (Table S2). Aliquots of 10 ng were used to quantify specific transcripts in Lightcycler® 78 480 Real Time PCR System (Roche, Germany) using Lightcycler 480 SYBR Green I Master® (Roche, Germany). Relative abundance values of all genes were calculated from the second derivative of their respective amplification curve, Cp values calculated by technical triplicates. Cp values of target genes were compared to the corresponding reference genes.

Data analysis
Within each deployment date, post-exposure feeding rates, gene and enzymatic activity responses at the studied sites were converted to proportional responses relative to the surrogate lab controls. Proportional responses of all the studied 12 sites were then compared using one way ANOVA followed by post-hoc Tukey's multiple comparison test (Zar, 1996). Prior to analysis data was log transformed to meet ANOVA assumptions of normality and variance homoscedasticity.
To explore causal-effect relationships between the studied environmental variables and parameters and the biological responses, Principal Component Analysis (PCA) and Partial Least Square Projections to Latent Structures regression (PLS) methods were used (Damásio et al., 2008). PCA was used to investigate the existing relationships between samples and variables and to deduce how many independent sources (components) were needed to explain the observed (experimental) data variance. These included two PCA performed on biological responses obtained in field (18) and lab exposures (12) and one conducted on environmental variables (44).
PLS is a regression extension of PCA which was used to connect the information between biological responses (Y variables) and environmental (X) variables. PLS models can be expressed in terms of traditional regression coefficients, b, of the multilinear model (i.e. y = X b model), simplifying model interpretation for the general case where multiple latent variables are needed for a satisfactory data modeling and prediction. Information about the correlation structure among variables and responses can be obtained using the VIP parameter. This parameter is a weighted sum of squares of PLS weights taking into account the amount of explained Y variance and it summarizes the information content of all latent and X variables (Wold et al., 1993) . X variables having high VIP scores contributed greater to Y variance.
We performed two PLS considering field and lab data. The latter analysis was limited to gene and post exposure responses and organic contaminant residues.
Since variables were very different and they were not measured using the same scale units, the data was auto-scaled prior to analysis (each element was subtracted by its column mean and divided by the standard deviation of its column). The number of PCA and PLS components was finally selected according to cross validation leaving one out prediction errors criteria (Wold et al., 2001). PCA and PLS analyses were conducted using the Matlab 6.0 software (MathWorks, Natick, Massachusetts). [TABLE1]

Environmental water parameters
From the 158 organic residues measured in the water samples collected from the 12 sites, 22 were not detected and the remaining ones were grouped into 37 groups according to their mode of action or chemical structure. In some cases the groups depicted in Table 1 contained a single compound whereas in others the sum of many. The composition of each group and its full name is depicted in Table S2 (Supplementary information). In general, conductivity and suspended solids increase from upper to downstream reaches being Llobregat the river having the highest levels. Water temperature was lower in Jucar sites. Measured oxygen levels in water and pH varied little across sites and were within ecological optimal values (Damásio et al., 2008). Selected contaminant groups varied dramatically across rivers and sites (Table 1).
For most contaminant groups the downstream site of Llobregat (L7) had the greatest values that on average were around 10 times higher than those of the rest of sites. As expected the PCA analysis performed for physicochemical water parameters was dramatically influenced by site L7, which alone explained 52.4 % of data variance (Fig. 2 inlet graph). Without considering L7 and excluding four classes of contaminants that varied little across the remaining sites (warfarin, tricyclic antihistamins, tamsulosin, clopidogrel), the PCA resolved in five interpretable components that explained up to 79 % of data variance (Table 2). Sample scores obtained from each of the five principal components were significantly correlated with 31 out of the 44 environmental factors considered. The first component (PC1) explained up to 34.7% of data variance, separated the three rivers (Fig. 2) and its sample scores were significantly (P<0.05) correlated with measured suspended solids, four metals, eight pharmaceutical groups, two industrial chemical groups (triazoles and alkylphenols), opiod residues of illicit drugs, azole fungicides and dexamethasone. PC2 explained up to 14.8% of data variance, and mainly separated E2 from the rest due to its low levels of conductivity, Ni and bisphenol A. The third, fourth and fifth components explained 30% of data variance and their sample scores were significantly correlated with the values of 8 environmental factors (Table 2). [TABLE2]

Biological responses
The mean percentage of animals recovered (dead and alive) from the chambers after field exposures were always higher than 95% and so mean survival rates. Except for carboxylesterase, the studied proportional responses in field transplanted D. magna varied significantly (P<0.05, ANOVA tests) across sites and rivers (see Table S4 for further information). The PCA analysis performed on the remaining 18 significant (P<0.05) D. magna proportional responses and the Shannon's diversity index resolved into five interpretable components that explained 89.8% of data variance (Table 3). The sample scores of the five obtained components correlated significantly (P<0.05) with the mean proportional responses of the studied biological variables establishing six clusters. Post exposure feeding and ChE proportional responses decreased in Llobregat and Ebro rivers from upstream to downstream reaches but not in Jucar and co-varied with the Shannon's diversity (Fig. 3A). The rest of D.
magna responses were clustered into five additional groups according to their variation across sites ( Fig. 3B-E). Interestingly proportional responses of gene markers but those of the gene

RXR grouped differently from biomarkers. [FIGURE 3][TABLE3]
From the 14 proportional biological responses determined for lab exposures only 12 varied significantly (P<0.05, based on ANOVA) across sites (see Table S5 for further information).
PCA conducted on those 12 variables identified three interpretable components explaining 84.1% of data variance (Table 3). PC1 explained most variance of data (58%) and its sample scores correlated significantly with mean proportional values of 9 out of the 11 genes considered. Sample scores of PC2 correlated with mean proportional responses of MRP4 and MT2 and proportional post exposure feeding rates with PC3's sample scores. The resulting five interpretable clusters of laboratory D. magna responses are depicted in Fig. 4.

[FIGUE 4]
Bi-variate Pearson correlations between mean site proportional responses of D. magna individuals transplanted in the field or exposed in the lab to organic eluates of water samples collected at the studied sites showed significant (P<0.05) correlations for EcR, MIH, FAA, UGP, MRP4 (Table 4).

[TABLE4]
3.4 Relationships between environmental parameters and biological effects PLS models performed to relate biological responses with environmental variables explained from 84 to 97% of the total variance of the former. In all analyses two PLS components were selected. PLS results for field variables are depicted in Table 5. For clarity we only depicted regression coefficients associated to VIP scores higher than 1. For interpretation purposes a cluster analysis was performed on regression coefficients, which differentiated eight clusters for biological and environmental variables, which are depicted in Table 5. Proportional responses of gene markers encoding for specific metabolic and molting processes (PGP, EcR, FAA, MIH, UGP, MRP4) were negatively related with up to 18 classes of contaminants.
Proportional feeding rates, and those of GST, LDH, CAT and RXR were negatively and positively related with measured environmental factors. Up to 9 classes of organic contaminants, 6 metals and water temperature were positively related with responses of HSP70 and/or MT2. Nine environmental variables including (estrogenic compounds, antibiotics, β2-adrenergic agonists) were related positively with species diversity and proportional cholinesterase activity whereas fungicides, pyrine, analgesics, conductivity and suspended solids were negatively related with HSP70 and MT2.
Environmental variables were grouped into three major and five minor clusters. One of the major clusters included pharmaceuticals (analgesics, lipid regulators, anticoagulants, sulfonylurea diuretics), opiod's illicit drugs, triazole and alkylphenol compounds, and covaried negatively with the response of 9 genes, and positively with that of three genes, feeding and three biomarkers. The cluster that included the metals As, Mn, Ni, diuretic, psychiatric drugs and anti-asthma compounds were positively related with the response of MT2 and HSP70 and negatively with that of MRP4, MIH, UGP. The third major cluster included antibiotics, illicit drugs and caffeine that were negatively related with the response of most traits except with those of CHE and H that covaried positively with them. From the remaining minor clusters and environmental variables those having greater coefficients were those of estrogenic compounds that were positively related with the response of VTG, CHE, H; the cluster of fungicides and antipyrine analgesics that covaried negatively with the response of CHE and H; opiate analgesics that covaried negatively with responses of CAT and H: and temperature, which was positively related with feeding and RXR responses.  (Table 4). This means that the relative response of eight genes to the studied organic contaminant classes was similar in field and lab exposures.

Discussion
This study aimed to identify biological active pollutants using in situ D. magna responses. It was implemented in three Mediterranean rivers differing in anthropogenic pressures and hence on pollution impacts. Indeed physicochemical water parameters such as water temperature, conductivity, suspended solids, five metals and 23 organic chemical functional classes differentiated the three rivers and within each river most and less polluted sites. Eighteen D. magna responses measured in individuals deployed in the field also differentiated the three rivers and sites. Biological responses clustered into six distinct groups. Post exposure feeding rates and cholinesterase activity co-varied similarly with the diversity of macroinvertebrate communities from the studied sites. Previous studies performed in the Llobregat and Besós rivers using D. magna and field collected caddisfly larvae also found a good correlation between the ecological quality of riparian invertebrate communities, D. magna feeding rates and cholinesterase activity in both D. magna and caddisfly larvae (Damásio et al., 2011;Damásio et al., 2008). Thus in situ D. magna responses of feeding rates and cholinesterase activity can be considered good markers of ecological quality. Indeed post-exposure feeding rates have already been used in several studies to assess ecological quality of river biota (Barata et al., 2007;Damásio et al., 2011;Mc William and Baird, 2002;Puértolas et al., 2010).
A second cluster included three biomarkers and the retinoic X receptor gene whose levels were inhibited towards downstream locations in Llobregat and Ebro river but increased in Jucar. Wang et al (2017) found that mRNA levels of the retinoic X receptor were high in reproductive females and were inhibited by insecticide terpenoids like pyriproxyfen. Thus, higher levels of RXR may indicate optimal conditions for Daphnia growth and reproduction and low levels increasing concentrations of insecticides. Enzymatic activities of CAT, GST and LDH followed similar response patterns than RXR thought less apparent. In the Llobregat river these enzymatic activities hardly varied across sites (values approached 1) but in downstream sites of the Ebro river increased, whereas in most sites of the Jucar river decreased. These results are consistent with previous studies that also reported no changes of CAT and GST activities across the same stations of Llobregat (Damásio et al., 2011) .
Inhibition or enhanced enzyme activities of CAT and GST in D. magna deployed in the field have been related to the presence of organophosporous pesticides, herbicides, alkylphenols, fungicides, metals and polycyclic aromatic hydrocarbons (Barata et al., 2007;Damásio et al., 2008). Another cluster was composed of genes encoding for specific responses such as sugar metabolism (UGP), molting (MIH) and xenobiotic transporter activity (MRP4). The responses of these genes varied across rivers and sites being down-regulated in Llobregat, up-regulated in Ebro and down and up-regulated in Jucar. There is reported evidence that UGP and MRP4 in D. magna are deregulated by psychiatric drugs and metals (Campos et al., 2014;Campos et al., 2013). A fourth cluster grouped genes codifying two stress proteins MT2 and HSP70 (Ho, 2008;Poynton et al., 2007). Metallothioneins like MT2 are involved in metal detoxification, binging to them and hence facilitating metal metabolism, whereas HSP70 protect other proteins under stress (Asselman et al., 2013;Bond and Bradley, 1997;Haap and Köhler, 2009;Haap et al., 2008). Levels of mRNA of these two genes and specially those of MT2 were upregulated in downstream sites of Llobregat that also had the highest levels of metals and of most other pollutants. A fifth cluster included genes having distinct functions such as those encoding for xenobiotic transporter proteins (PGP), those involve in the metabolism of amino acids (FAA) and lipids (THIO) and the gene encoding for the ecdysone receptor (EcR) (Campos et al., 2014;Campos et al., 2013;Kato et al., 2007). The response of these genes was up-regulated at the downstream sites of Llobregat and Jucar and down-regulated at E4 in Ebro river. Previous studies have reported that pharmaceuticals and pesticides deregulated those genes (Campos et al., 2014;Campos et al., 2013;Mu and Leblanc, 2004;Wang et al., 2011).
A sixth cluster included vitellogenin and two genes from the Krebs cycle (ACON, IDH). mRNA levels of those genes varied between rivers and across sites as follows: within river basins mRNA gene levels decreased towards downstream sites in Llobregat, remained unchanged in Ebro and increased towards downstream reaches in Jucar. Note also that levels of VTG were always down-regulated relative to the surrogate lab control having the greatest levels of deregulation in Ebro. Hannas et al. (2011) recently reported that the putative VTG gene in D. magna acts like a general stress gene that could be either up or down-regulated by many contaminants. There is also evidence that the transcripts of the two Krebs cycle genes (ACON, IDH) are deregulated by pollutants (Campos et al., 2013).
A further characterization of the relationships between measured responses and environmental factors was performed with the aid of the PLS excluding out L7. The most important relationships allowed to group biological responses into eight groups that were affected similarly by environmental parameters. These include cholinesterase and diversity responses that were affected negatively by conductivity, suspended solids in water, fungicides and antipyrine analgesics and positively by -blockers, β2-adrenergic agonists and antibiotics.
Residue levels of antipyrine analgesics and β2-adrenergic agonist were probably too low ≤2 ng/l to have any effect on the measured responses and those of -blockers hardly varied across most sites. In this study conductivity measured salinization that in Mediterranean rivers is an important problem that deteriorates water quality (Damásio et al., 2011). High levels of suspended solids associated to anthropogenic impacts are also known to impair the ecological quality of rivers (Damásio et al., 2011). Fungicides are also known to reduce diversity of riparian communities directly impairing the physiology of aquatic organisms and indirectly by reducing fungus biomass and hence food for shredders (Maltby et al., 2009). There is evidence that low levels of antibiotics reduce the microbial load of invertebrates, promoting its growth and reproduction (Zalewski et al., 2011). Thus low levels of antibiotics could be beneficial for invertebrates and thus may increase diversity. Previous studies have reported that organophosphorous or carbamate insecticides inhibit cholinesterase acticities of D. magna at about 100 ng/L ( Barata et al., 2007;Barata et al., 2004). Therefore, measured organophosphorours and carbamate pesticides residue levels in water were probably too low (≤2 ng/l) to have impaired cholinesterase activity in Daphnia. The cluster of the remaining biomarkers (CAT, GST, LDH) was related positively with a cluster of eight organic chemicals (i.e.analgesics, lipid regulators, diuretics, illicit drugs, trizoles) and negatively with up to 11.
Excluding out pesticides, sublethal effects of organic chemical substances to D. magna rarely occur below µg/l (Constantine and Huggett, 2010;Yang et al., 2013). This means that from the above mentioned pollutant groups, residues of trizole compounds, NSAID analgesics and lipid regulators were the most likely to affect biological responses of deployed D. magna individuals. Nevertheless for these chemical groups effects on D. magna have been reported at the mg/l range (Heckmann et al., 2007;Seeland et al., 2012;Zurita et al., 2007).
The cluster of metallothionein and heat shock protein 70 was related positively with most studied environmental factors. This is expected to occur since these two proteins are known to be induced under stress acting as a detoxification mechanisms to metals (MT2) or protecting proteins (HSP70). (Asselman et al., 2013;Bond and Bradley, 1997;Haap and Köhler, 2009;Haap et al., 2008). Interestingly MT2 responses obtained in the field were not related with those observed in lab exposures to organic eluates neither their PLS regression coefficients.
This means that responses of MT2 in the field were related to other factors than organic pollutants. Metals like As, Mn, Ni and conductivity were positively related with MT2, which support recent findings reporting that mRNA levels of MT2 are inducible by several metals (Asselman et al., 2013). Genes encoding for key enzymes of the Krebs cycle such as IHD and ACON, VTG and those involved in the lipid (THIO) or xenobiotic metabolism (PGP) responded differently in the field and lab exposures. According to our experimental design biological responses to lab exposures should mimic those of the field for organic chemicals.
This means that for those responses, the obtained PLS-relationships with measured organic residues have to be considered with caution since they can be associated to other sources of contamination than those measured in this study. This was the case for mRNA levels of vitellogenin in D. magna that was positively related with estrogenic steroids but not with other known estrogenic compounds such as alkylphenols, bisphenol A and triazoles (De Castro-Català et al., 2013). Furthermore, obtained PLS associations of vitellogenin in field exposures were not correlated with those obtained in lab exposures. In a previous study De Castro-Català et al. (2013) reported a positive relationship between estrogenic compounds and the number of eggs per clutch in freshwater snails deployed to the same rivers and sites. Contrary to crustaceans like Daphnia (Hannas et al., 2011), gastropods respond to estrogenic compounds (Castro et al., 2007;Stange et al., 2012). Differences in the exposure scenario and/or variations in gene responses to environmental stressors may also explain the observed lack of relationships between field and lab exposures. During field assays animals were directly exposed to the river water flow and suspended matter and hence to a broader number of contaminants and other environmental factors than in lab exposures, which were conducted under static conditions and using organic extracts of filtered water samples. Genes encoding for general metabolic processes (IHD, ACON, THIO) are likely to be altered by several stressors apart from organic contaminants and that encoding for PGP is known to respond to many different organic and inorganic chemicals (Campos et al., 2014;Faria et al., 2011).
Therefore, the higher complexity of field exposures may explain the observed discrepancies on gene responses between lab and field. This was clearly illustrated for feeding responses, which were more affected during field assays (Fig 3, 4 A).
For genes related with molting and reproductive processes (MIH, EcR), sugar, protein and xenobiotic metabolism (UGP, FAA, MRP4), lab and field responses were similar and were negatively related with several pharmaceutical groups (diuretics, illicit drugs, lipid regulatorsfibrates, psychiatric drugs, β2-adrenergic agonists), bisphenol A, estrogenic and anti-parasitic compounds, As, Mn, Ni, Co, conductivity and suspended solids. Few of the measured contaminants (trizoles, pesticides) were positively related with MRP4, MIH, FAA or/and UGP. Responses of the retinoic receptor RXR followed a distinct patern being positively and negatively related with the above mentioned and other pollutants like caffeine, analgesic opiates. An excess of suspended solids, salt and of metals like Ni, inhibit growth or/and reproduction in D. magna (Baillieul et al., 1996;Burton Jr et al., 2005;Diamond et al., 1992;Hoang et al., 2007;Pane et al., 2004;Robinson et al., 2010). Psychiatric drugs such as selective serotonine re-uptake inhibitors are known to disrupt signaling gene pathways of molting, reproduction, sugar and aminoacid metabolism at environmental relevant concentrations close to 1 µ/l (Campos et al., 2013). Despite that sublethal effects of fibrates, estrogenic steroids and bisphenol A occur in D. magna at mg/l (Brennan et al., 2006;Jeong et al., 2013;Zurita et al., 2007), in mixtures pharmaceuticals may trigger physiological responses at lower doses (Cleuvers, 2003). Many organic contaminants and metals can interact and act additively or synergically inhibiting xenobiotic transporter activity mediated by multidrug resistance proteins like MRP4 (Campos et al., 2014;Faria et al., 2011;Kurelec, 1997). In the present study many of the above mentioned pharmaceutical and industrial chemical groups occurred together and hence it should be feasible to establish a threshold biological effect at 100 ng/l. Accordingly analgesics, diuretics, psychiatric drugs,  blockers, illicit drugs, trizoles, caffeine and measured pesticide levels should be considered of environmental concern.

Conclusions
In summary, our results led positive support to the use of sublethal specific gene responses in combination with in situ Daphnia feeding and biochemical responses to assess effects and identify in a better way environmentally detrimental factors within a complex multi-stressed river systems in the field, thus contributing to a more realistic assessment of ecological risks.
D. magna responses of feeding and cholinesterase activity were able to assess the general ecological quality of the selected river sites and those of genes assessed specific effects of particular contaminant groups. Interestingly the contaminant groups that differentiated the studied sites (Table 2) did not follow those that were associated with specific biological response (Table 5). This means that risk assessment estimates based on chemical analyses have to be taken with caution (Fàbrega et al., 2013;Ginebreda et al., 2010). Furthermore, subdivision of chemical groups according to known mode of actions and the inclusion of lab exposures allowed to judge false positive relationships. The experimental procedures developed in this study indicate that in multi-stressed rivers biota is often chronically exposed to sublethal levels of contaminants and hence a large set of biological and chemical markers are needed to identify detrimental pollutants. Nevertheless, it is important to consider that, while the approach used in this study drives us closer towards the "in situ environmental hazard identification evaluation", issues arising from other confounding factors influencing in situ Daphnia responses still should be considered with caution in the interpretation of such findings as conclusive diagnostic proofs of individual factors causing effects. Table 1. Measured (Mean ± SD) physco-chemical water parameters across the studied river sites. Contaminant residues have been grouped in classes. See Table S2 for further explanation. Metal and organic contaminant levels are µg/l and ng/l, respectively. The range (min-max) is also depicted.        Table S4. Proportional biological responses (Mean SE; N =5) measured in D. magna individuals deployed at the studied sites. Different letters 1 means significant (P<0.05) differences across sites following ANOVA and Tukey's multiple comparison tests. The number of species (S) and