# Whole-transcriptome response to wastewater treatment plant and stormwater effluents in the Asian clam, *Corbicula fluminea*

Bertucci Anthony <fnref n="1" />, *, Pierron Fabien <fnref n="1" />, Gourves Pierre-Yves <fnref n="1" />, Klopp Christophe <fnref n="2" />, Lagarde Gauthier <fnref n="1" />, Pereto Clement <fnref n="1" />, Dufour Vincent <fnref n="1" />, Gonzalez Patrice <fnref n="1" />, Coynel Alexandra <fnref n="1" />, Budzinski Helene <fnref n="1" />, Baudrimont Magalie <fnref n="1" />

## Abstract :

The increase in human population and urbanization are resulting in an increase in the volume of wastewater and urban runoff effluents entering natural ecosystems. These effluents may contain multiple pollutants to which the biological response of aquatic organisms is still poorly understood mainly due to mixture toxicity and interactions with other environmental factors. In this context, RNA sequencing was used to assess the impact of a chronic exposure to wastewater treatment plant and stormwater effluents at the whole-transcriptome level and evaluate the potential physiological outcomes in the Asian clam *Corbicula fluminea*.

We de-novo assembled a transcriptome from C. fluminea digestive gland and identified a set of 3,181 transcripts with altered abundance in response to water quality. The largest differences in transcriptomic profiles were observed between C. fluminea from the reference site and those exposed to wastewater treatment plant effluents. On both anthropogenically impacted sites, most differentially expressed transcripts were involved in signaling pathways in relation to energy metabolism such as mTOR and FoxO, suggesting an energy/nutrient deficit and hypoxic conditions. These conditions were likely responsible for damages to proteins and transcripts in response to wastewater treatment effluents whereas exposure to urban runoff might result in immune and endocrine disruptions.

In absence of comprehensive chemical characterization, the RNAseq approach could provide information regarding the mode of action of pollutants and then be useful for the identification of which parameters must be studied at higher integration level in order to diagnose sites where the presence of complex and variable mixtures of chemicals is suspected.

<footnote n="1">¹ Univ. Bordeaux, UMR EPOC CNRS 5805, 33615 Pessac, France</footnote>

<footnote n="2">² Plate-forme bio-informatique Genotoul, Mathématiques et Informatique Appliquées de Toulouse, INRA, 31326 Castanet-Tolosan, France</footnote>

<footnote>* Corresponding author : Anthony Bertucci, email address : anthony.bertucci@u-bordeaux.fr</footnote>
**Graphical abstract**

<figure><img src="image_1.png" /></figure>

**Highlights**

► Clams were caged for 3 months in a river receiving WWTP and stormwater effluents. ► We de-novo assembled a transcriptome from clams' digestive glands. ► 3181 genes were differentially expressed between sites. ► Water pollution induced disorders to energy metabolism and immunity. ► RNAseq allowed the identification of novel early-warning molecular markers.

**Keywords** : Transcriptomics, Wastewater treatment plant, Urban runoff, Water quality, In-situ exposure, Corbicula fluminea
44

# 1. Introduction

45

46

47

48

49

Aquatic ecosystems are threatened by the rapid growing of human population and the parallel increase in urbanization that are driving the increase in the volume of wastewater and urban runoff worldwide (Muñoz et al., 2008; Sato et al., 2013). Direct and indirect effects of chronic exposure to effluents from Wastewater Treatment Plant (WWTP) and urban runoff effluents and the consequences on ecosystem functioning are still poorly understood.

50

51

52

53

54

55

Wastewater treatment plants effluents are characterized by inorganic and organic materials such as metals, biocides, polycyclic aromatic hydrocarbons (PAHs), chlorinated solvents, pharmaceuticals and personal care products (Arlos et al., 2015; Miège et al., 2009; Muñoz et al., 2008). For both pharmaceuticals and personal care products, wastewater is the main route of emission to the environment (Ternes et al., 1999). Negative effects of WWTP effluents were documented on aquatic organisms (including bivalves), from the molecular to the community level (Morris et al., 2017).

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

WWTP effluent exposure were often associated with endocrine disruption and altered reproduction (Fuzzen et al., 2015; Mezghani-Chaari et al., 2015), immune response disturbances (Jasinska et al., 2015) and oxidative stress (Gillis et al., 2014). Urban runoff (also known as stormwater), the precipitation-related discharge of impervious surfaces, is considered another major source of water pollution (Huber et al., 2016). Major pollutants found in urban runoff include both organic and inorganic substances, such as nutrients, organic and particulate matter, PAHs (Krein and Schorer, 2000), metals such as Pb, Zn and Cu and pathogens (Sidhu et al., 2012). *Salmonellae* that is found in humans, pets, farm animals, and wild life feces is the most common pathogen. Municipal sewage, agriculture pollution, and storm water runoff are the main sources in natural waters (Cabral, 2010). A list of 25 selected stormwater priority pollutants was proposed by Eriksson *et al*. (Eriksson et al., 2007). These pollutants may have multiple origin and a high level of qualitative and quantitative variability. The characteristics of runoff waters (physical, chemical or microbial) are dependent on the type of surfaces encountered and their use (Göbel et al., 2007). However, main contamination sources are related to traffic, surrounding land use, atmospheric contamination, meteorology and other environmental factors (Huber et al., 2016). Moreover, the concentration of each single constituent,
71 likely to represent a risk for human, animals or plants, may vary from one precipitation event to the
72 other and from site to site (Kayhanian et al., 2008; Kayhanian Masoud et al., 2003).

73 Transcriptomics (thanks to the development of next generation - high throughput - sequencing), like
74 other "omics", can provide a global investigation of the potential toxicity of contaminated water, as
75 well as the discrimination of clean and polluted sites. Indeed, mRNA sequencing (RNAseq) data may
76 depict the alteration of several genes and molecular pathways simultaneously in non-model organisms
77 in response to stressful conditions (Gonzalez and Pierron, 2015). This method is more and more often
78 used in a multi-stress context (Baillon et al., 2015; Bertucci et al., 2017; Pierron et al., 2011; Poynton
79 et al., 2008; Regier et al., 2013).

80 The Asian freshwater clam *Corbicula fluminea* is an invasive species, widely distributed in rivers,
81 lakes and estuaries throughout Europe (Sousa et al., 2008). This suspension feeder has been used as a
82 bioindicator to assess the impact of different contaminants, such as metals (Arini et al., 2014) or
83 organic molecules (Bebianno et al., 2016). In the present study, we investigated the transcriptome-
84 wide response of field caged specimens to wastewater treatment plant and stormwater runoff effluents
85 in order to get insights on the impacted metabolic pathways and to evaluate the potential use of
86 transcriptomics as early diagnostic tool of water quality *in-situ*.
87

87

## 2. Materials and methods

88

### 2.1. Study organism

89

90 The freshwater clam (Bivalvia: Corbiculoidea) *Corbicula fluminea* (Müller, 1774) is an invasive
91 species original from Asia that has spread through most Europe and America. This benthic species is
92 considered a good model organism due to its wide distribution, abundance, ease of collection and
93 maintenance under laboratory or *in-situ* conditions. Moreover, it filters large volumes of water for
94 respiratory and nutritive purposes. It was used in ecotoxicological studies to evaluate the effects of a
95 variety of pollutants such as metals, flame retardants, pesticides and pharmaceuticals. Adult specimens
96 (approximatively 2 cm long) were originally from a reference site on the Isle River (45°0.878′N,
97 0°9.820′E).

98

### 2.2. Study sites
Three field sites located along the Jalle d'Eysines River (also known as the Jalle de Blanquefort) close to Bordeaux, France (Figure 1) were selected. Site 1 (44°53.652'N, 0°41.534'W), located upstream, was chosen as reference site owing to its use as a source of drinking water. Intermediate Site 2 (44°53.732'N, 0°39.916'W) is located 100 metres downstream of a wastewater treatment plant. The most downstream site, Site 3 (44°54.035'N, 0°36.803'W) receives inlet of urban runoff water from surrounding area and from Bordeaux's ring road. Three plastic racks, each containing 20 clams were transplanted at each site (Salazar and Salazar, 2005) for 3 months (14th December 2015 – 14th March 2016). These dates are outside of the reproduction period of *C. fluminea* (i.e. May to October; Mouthon, 2001). In December 2015 to March 2016, rainfall events were registered on 71 days with a monthly average precipitation of 127 mm (www.infoclimat.fr). Sites 1, 2 and 3 will be referred to as Reference, WWTP exposure and urban runoff exposure, respectively.

110

<figure><img src="image_2.png" /></figure>

111

Figure 1: Overview of the study area. The sites are figured with red paddles numbered 1-3. Green, yellow and red diamonds indicate the reference site, wastewater treatment plant and urban runoff effluent, respectively. Main localities are indicated. Major roads are in red.

115

116

## 2.3. Physico-chemical parameters and pollution analysis

117

Turbidity, temperature, pH, conductivity, RedOx potential and dissolved oxygen (DO) concentration
were measured once a day on each site during the entire exposure period with a multiparameter probe.

The collection of water samples for pollutants analysis was only performed episodically on 5 days
scattered in February and March 2016. For trace metal analysis, water samples were collected in acid
pre-cleaned (72 h in 1.5 N HNO3, and rinsed with ultrapure water) polypropylene bottles. All samples
were immediately filtered on-site through 0.2 µm Sartorius® polycarbonate filters. Filtrates were
collected in pre-cleaned 30 ml polypropylene bottles, acidified (0.5%; HNO3 PlasmaPur Ultrex®) and
stored at 4°C until analysis. Samples were analyzed by ICP-MS (Xserie2, Thermo Scientific) under
standard conditions. The applied analytical methods were continuously quality checked by analysis of
international certified reference materials. Accuracy was within 10% of the certified values and the
analytical error lower than 5 % (r.s.d.) for concentrations ten times higher than detection limits.

Water samples for pesticide analysis were collected with HDPE bottles (Nalgene). Considering the
slight amounts of suspended matter in the river (usually between <2 and 10 ng.L-1) and the polar
nature of the monitored pesticides, samples were filtrated over 0.7 µm glass filters from Wattman
(Fisher Bioblock Scientific, Illkirch, France), and analysis were performed only on dissolved phase.

Waters were then stored at -18°C until extraction. Glyphosate and aminomethylphosphonic acid
(AMPA) were extracted following the protocol described by Fauvelle et al. (2015). Glyphosate-
13C215N and AMPA-13C15N were introduced in samples (3 ng) to play role of internal standard. 40
µL of the final extract were injected in LC-MS/MS (HPLC 1290 and triple quadrupole 6460 from
Agilent Technologies (Santa Clara, CA, USA). Separation was performed at 0.6 mL.min-1 over a
reverse C18 phase Acquity BEH column (50 x 2.1mm; 1.7µm, Waters) at 35°C with ultrapure water
and ammonium acetate 5mM added with ammonium acetate until pH 9, and MeOH as mobile phases.

Acquisition was performed using electrospray in negative mode in dynMRM. Fipronil was
characterized following the method described by Le Coadou et al. (2017), with on-line SPME coupled
with an Agilent 7890 gas chromatography system and a 7000C tandem mass spectrometer from
Agilent-Technologies (Santa Clara, CA, USA). 9 mL of sample were transferred into SPME auto
sampler vial, spiked with 0.3 ng of fipronil-13C15N as internal standard. Extraction was performed for
30 min using a PDMS/DVB fiber (65µm), at 50°C under agitation (250 rpm). Analytes were desorbed
for 10 min in the injector of the chromatograph set at 250°C in pulsed splitless mode (25 psi, 1 min).
Separation was performed using an Agilent J&W GC column HP 5 MS UI (30 m x 0.25 mm x 0.25 µm film thickness; 5% phenylmethylsiloxan). Diuron, hexazinone, isoproturon, metolachlor, metolachlor OA and terbutryn were extracted by solid phase extraction, following the protocol described by Belles et al. (2014) concerning polar pesticides. The extraction was performed on 100 mL of sample, acidified at pH 2 with chloridric acid. 5 ng of diuron-d6, hexazinone-d6, isoproturon-d6, metolachlor-d6 and terbutryn-d5 were added in samples to play role of internal standard. 5µL of the final extracts were injected in LC-MS/MS (HPLC 1290 and triple quadrupole 6460 from Agilent Technologies (Santa Clara, CA, USA). Separation was performed at 0.5 mL.min⁻¹ over a reverse C18 phase Kinetex column (100 x 2.1 mm; 1.7 µm, Phenomenex) at 35°C with ultrapure water (+5mM ammonium acetate and 0.1% acetic acid) and MeOH as mobile phases. Acquisition was performed using electrospray in positive mode in dynMRM. Performance of each protocols were evaluated for each batch with the extraction of artificial samples based on mineral water spiked with known amount of studied molecules. Quantification limits and quantification yields were calculated on these samples to ensure the performances of the different protocols. Blanks were also extracted in parallel to control the potential contamination introduced during extraction procedures and analysis.

## 2.4. Clams collection

Eight individuals were collected randomly and immediately dissected at each site at the end of the exposure period, leading to a total number of 24 individuals. Soft (body) tissues and shells were weighted in order to determine the condition index (CI) with the following formula: Tissue (g, fresh weight) / Shell (g, dry weight) * 100 (Quayle and Newkirk, 1989). The digestive glands were then dissected for RNA extraction. In bivalves, the digestive gland functions both as a site for nutrients / pollutants uptake and as an important reservoir for contaminant storage. Samples were collected in RNA-later solution and kept on ice until return to the lab where they were placed at 4°C overnight before storage at -20°C until RNA extraction.

## 2.5. RNA Extraction, library construction and sequencing

Total RNA was isolated from digestive glands tissue using the SV Total RNA isolation system kit (Promega) according to the manufacturer's instructions. RNA-seq libraries were prepared and sequenced by Genotoul (France) with the Illumina HiSeq 3000 technology. The RNA-seq raw reads
174

used in this study have been deposited in the NCBI Gene Expression Omnibus (Edgar et al., 2002)

175

under GEO Accession GSE104933.

176

## 2.6. Bioinformatics workflow

177

The read quality of the RNA-seq libraries was evaluated using FastQC

178

(http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Reads were cleaned and filtered using the DRAP pipeline (Cabau et al., 2017) and *de novo* assembled with the Trinity assembler version

179

2.0.6 (Grabherr et al., 2011). Assembled contigs were filtered in order to keep only those with at least one fragment per kilobase of transcript per million reads (FPKM). Reads have been realigned back to contigs with bwa (Li and Durbin, 2009) (version 0.7.12). The resulting files were compressed, sorted and indexed with samtools (Li et al., 2009) (version 1.1). The contig transcription counts have been generated with samtools (version 1.1). The alignment files have then been filtered for duplicates with samtools (version 1.1). This Transcriptome Shotgun Assembly project has been deposited at

181

182

183

184

185

186

DDBJ/EMBL/GenBank under the accession GFYR00000000. The version described in this paper is the first version, GFYR01000000. The assembled contigs were aligned with NCBI blast (Camacho et al., 2009) (version 2.2.26) on Refseq protein, Swissprot and Ensembl protein reference files from *Crassostrea gigas*, *Lottia gigantea* and *Lingula anatina* to retrieve the corresponding annotations. The contigs were also processed with rnammer (Lagesen et al., 2007) (version 1.2) to find ribosomal genes, with repeatmasker (Smith et al.) (version open-4-0-3) to list the contained repeats and with interproscan (Quevillon et al., 2005) (version 4.8) for gene ontology (GO) and structural annotation.

187

188

189

190

191

192

193

GO annotations were obtained with an E-value and a homology threshold set to 10⁻⁶ and 50%, respectively. Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic pathways assignments were performed with the the GhostKOALA server (Kanehisa et al., 2016). Detailed statistic information of *de novo* assembly are listed in Table 1. Moreover, all the results have been uploaded in a RNAbrowse instance (Mariette et al., 2014) (Sigenae web address available upon publication).
Transcriptome completeness was assessed with BUSCO v3 (Simão et al., 2015) against the Metazoa orthologs dataset.

194

195

196

197

198

199

200

## 2.7. Real-time PCR
Specific primers amplifying approximately 100 bp were designed for a set of 18 genes using the software primer3 (Rozen and Skaletsky, 2000) (Details are given in Appendix A). Expected lengths of the amplicons were checked by agarose gel electrophoresis after regular PCR amplification using the GoTaq DNA Polymerase (Promega). Primer efficiencies were determined using standard curve analysis (E = 10^(-1/slope)) with a dilution series of pooled cDNA from all conditions and ranged from 90 to 99 %. Transcript level quantification was performed using the GoTaq qPCR mastermix from Promega and a Startagene Mx3000P system. PCR conditions were as follows: 1X GoTaq qPCR mastermix (Promega), 200 nM primers and 10 ng of cDNA in a total volume of 25 µL. PCR parameters were: 95 °C for 2 min, followed by 40-45 cycles of 15 s at 95 °C, 60 s at 60 °C and a dissociation curve step (60-95 °C) to confirm the absence of nonspecific products. The dissociation curves showed a single amplification product and no primer dimers. Three contigs were selected as control genes based on their stability in the RNA-seq data set. The relative quantification for each gene was normalized by the geometric mean of control genes and relative to its expression in the reference site. Fold expression values were calculated as (E^dCT)_target/(E^dCT)_control where E is the amplification efficiency for each pair of primers and dCT = CT-site - CT-reference.

## 2.8. Statistical analyses

The comparison of CIs amongst the three sites was performed by analysis of variance (ANOVA), after checking assumptions of normality and homoscedasticity, followed by a Tukey HSD test. As the assumptions of the parametric ANOVA were not met for physico-chemical parameters we used the Scheirer-Ray-Hare (SHR) non-parametric ANOVA to test the differences amongst sites. Pairwise comparisons within the 3 sites were carried out using Mann-Whitney U-tests. Computations were performed with the R package "Stats". Differences were considered significant at the level p<0.05. Since the number of measures didn't allow proper statistical comparison of micropollutants, minimal and maximal values only were presented.

Differential gene expression was calculated with the R package DESeq2 (Love et al., 2014). P-values were corrected using the Benjamini and Hochberg method, and a false discovery rate (FDR) threshold of 0.05 was used. To avoid large fold changes in low abundance transcripts, only contigs with more than 2 counts per million (cpm) in at least one condition were considered in subsequent analyses
229 (Appendix C). GO-enrichment analyses were carried out with R package topGO (Alexa and
230 Rahnenfuhrer, 2016). In order to reduce redundancy in GO terms and ease the interpretation, we used
231 the REVIGO web server with a value of C=0.5 (Supek et al., 2011). KEGG enrichment was calculated
232 with a Fisher exact test. Both GO and KEGG analyses were performed with a p-value ≤ 0.05.

233

234 **3. Results and discussion**

235 **3.1. Transcriptome assembly**

236 RNA sequencing generated 795 million paired-end 150 bp reads. We assembled a total of 58,291
237 contigs of which 23,758 (40.76 %) were annotated. The assembly statistics are presented in Table 1.
238 The quality of our assembly is evidenced by the high abundance of Metazoa benchmarking universal
239 single-copy orthologues (BUSCO summary: C:98.7%[S:92.5%,D:6.2%],F:0.2%,M:1.1%,n:978). The
240 percentage of reads mapping the transcriptome ranged from 95.8 to 99.2 %. The oyster *Magallana*
241 gigas (previously *Crassostrea gigas*) accounted for 55.84 % (10,793 sequences) of the best BLAST
242 homology. The functional annotation showed that 13,608 (57.3 %) and 9,427 (36.7 %) of these contigs
243 were associated to GO and KEGG terms, respectively (details are given in Appendix B).

244
<table><thead><tr><th colspan="5">Statistics</th></tr><tr><th>Sequencing and mapping</th><th>Total</th><th>Reference</th><th>WWTP</th><th>Runoff</th></tr></thead><tbody><tr><td>Total number of Pair-End reads</td><td>794,791,848</td><td>291,118,523</td><td>245,734,147</td><td>257,939,178</td></tr><tr><td>Number of mapped reads</td><td>777,208,882</td><td>284,107,241</td><td>241,459,784</td><td>251,641,857</td></tr><tr><td>Min. number of reads per sample</td><td></td><td>30,314,069</td><td>20,657,826</td><td>23,463,430</td></tr><tr><td>Max. number of reads per sample</td><td></td><td>43,582,589</td><td>40,380,849</td><td>50,218,485</td></tr><tr><th colspan="5">Assembly</th></tr><tr><td>Number of contigs</td><td>58,291</td><td></td><td></td><td></td></tr><tr><td>Total contigs length (bases)</td><td>98,530,638</td><td></td><td></td><td></td></tr><tr><td>Min. contig length (bases)</td><td>207</td><td></td><td></td><td></td></tr><tr><td>Max. contig length (bases)</td><td>53,892</td><td></td><td></td><td></td></tr><tr><td>Mean contig size (bases)</td><td>1,690</td><td></td><td></td><td></td></tr><tr><td>N50 (bases)</td><td>2,401</td><td></td><td></td><td></td></tr><tr><td>Mean GC %</td><td>41.4</td><td></td><td></td><td></td></tr><tr><th colspan="5">Annotation</th></tr><tr><td>Number of contigs with BLAST hits</td><td>23,758</td><td></td><td></td><td></td></tr><tr><td>Number of contigs with GO terms</td><td>13,608</td><td></td><td></td><td></td></tr><tr><td>Number of contigs with KEGG ortholog</td><td>9,427</td><td></td><td></td><td></td></tr></tbody></table>

Table 1: *Corbicula fluminea*'s digestive gland transcriptome assembly statistics.

## 3.2. Water physicochemistry and pollutants analysis

Significant differences amongst sites were detected for conductivity, pH, Redox potential, DO and temperature (Table 2). However, variations were relatively small and unlikely to produce much metabolic modifications since they correspond to normal range of values and seasonal variations experienced by clams in similar environments (Vidal et al., 2002). For instance, temperature variation was less than 1°C and pH values at each site were in the range of the High Ecological Status class of the European Water Framework Directive (WFD, 2000/60/EC; EC, 2000), i.e. from 6.5 to 8.2. On the contrary, DO appeared as the most variable parameter with an average of 10.13 mg/l (89.39% saturation), 7.47 mg/l (67.18%) and 5.51 mg/l (49.34%) on reference, WWTP and runoff sites, respectively. WWTP and urban runoff could hence be classified as good and moderate, respectively whilst the reference site was of high ecological status (EC, 2000). Trace metal showed no differences amongst sites. Amongst the nine pesticides we measured, the largest differences were observed for glyphosate and its primary degradation product, AMPA (Aminomethylphosphonic acid), and for terbutryn that were higher in the WWTP and urban runoff effluents.
<table><thead><tr><th></th><th></th><th>Reference</th><th>WWTP</th><th>Runoff</th></tr></thead><tbody><tr><td rowspan="6">Physicochemistry</td><td>Turbidity (FNU)</td><td>14.09 ± 3.04</td><td>25.47 ± 6.77</td><td>26.38 ± 11.46</td></tr><tr><td>Conductivity (µS/cm)</td><td>330.50 ± 11.39 <sup>a</sup></td><td>396.22 ± 22.82 <sup>b</sup></td><td>420.88 ± 17.97 <sup>b</sup></td></tr><tr><td>pH</td><td>7.26 ± 0.07 <sup>a</sup></td><td>6.89 ± 0.03 <sup>b</sup></td><td>7.13 ± 0.03 <sup>a</sup></td></tr><tr><td>Redox potential (mV)</td><td>414.82 ± 4.86 <sup>a</sup></td><td>403.85 ± 3.84 <sup>a</sup></td><td>453.99 ± 1.71 <sup>b</sup></td></tr><tr><td>Dissolved O<sub>2</sub> (mg/l)</td><td>10.13 ± 0.08 <sup>a</sup></td><td>7.47 ± 0.47 <sup>b</sup></td><td>5.51 ± 0.33 <sup>c</sup></td></tr><tr><td>Temperature (°C)</td><td>9.77 ± 0.14 <sup>a</sup></td><td>10.57 ± 0.16 <sup>b</sup></td><td>10.38 ± 0.14 <sup>b</sup></td></tr><tr><td rowspan="10">Metal<br>(min-Max µg/L)</td><td>Al</td><td>101.9 - 277.7 (3)</td><td>82.1 - 261.2 (3)</td><td>84.1 - 140.3 (3)</td></tr><tr><td>Cd</td><td>0.02 - 0.03 (3)</td><td>0.02 - 0.03 (3)</td><td>0.01 - 0.03 (3)</td></tr><tr><td>Co</td><td>0.09 - 0.57 (3)</td><td>0.10 - 0.49 (3)</td><td>0.11 - 0.31 (3)</td></tr><tr><td>Cr</td><td>0.54 - 0.64 (3)</td><td>0.49 - 0.64 (3)</td><td>0.50 - 0.58 (3)</td></tr><tr><td>Cu</td><td>0.72 - 1.20 (3)</td><td>1.24 - 1.39 (3)</td><td>1.16 - 1.59 (3)</td></tr><tr><td>Fe</td><td>197.3 - 279.1 (3)</td><td>168.2 - 283.0 (3)</td><td>183.3 - 226.5 (3)</td></tr><tr><td>Ni</td><td>0.56 - 1.09 (3)</td><td>0.57 - 1.20 (3)</td><td>0.55 - 1.15 (3)</td></tr><tr><td>Pb</td><td>0.16 - 1.07 (3)</td><td>0.14 - 0.29 (3)</td><td>0.14 - 0.27 (3)</td></tr><tr><td>V</td><td>0.98 - 1.10 (3)</td><td>1.03 - 1.26 (3)</td><td>1.08 - 1.19 (3)</td></tr><tr><td>Zn</td><td>4.25 - 8.37 (3)</td><td>4.88 - 8.33 (3)</td><td>6.44 - 11.09 (3)</td></tr><tr><td rowspan="10">Pesticides<br>(min-Max ng/L)</td><td>AMPA</td><td>52.3 (1)</td><td>93.2 (1)</td><td>141.2 - 230.3 (2)</td></tr><tr><td>Diuron</td><td>2.6 - 5.9 (4)</td><td>3.8 - 7.3 (4)</td><td>4.3 - 14.1 (5)</td></tr><tr><td>Fipronil</td><td>0.5 - 0.8 (2)</td><td>1.0 - 1.8 (3)</td><td>1.0 - 1.2 (4)</td></tr><tr><td>Glyphosate</td><td>29.1 (1)</td><td>27.1 (1)</td><td>57.9 - 77.2 (2)</td></tr><tr><td>Hexazinone</td><td>2.3 - 17.8 (4)</td><td>2.0 - 16.5 (4)</td><td>1.6 - 14.9 (5)</td></tr><tr><td>Isoproturon</td><td>1.5 (1)</td><td>0.2 - 1.5 (2)</td><td>0.3 - 1.7 (3)</td></tr><tr><td>Metolachlor</td><td>8.9 - 20.1 (4)</td><td>8.3 - 17.7 (4)</td><td>5.7 - 16.3 (5)</td></tr><tr><td>Metolachlor OA</td><td>317.3 - 353.0 (4)</td><td>274.3 - 501.7 (4)</td><td>188.6 - 469.9 (5)</td></tr><tr><td>Terbutryn</td><td>0.3 - 0.8 (4)</td><td>0.6 - 5.2 (4)</td><td>0.6 - 5.4 (5)</td></tr></tbody></table>
Table 2: Physico-chemical parameters and micropollutants concentrations at each site. Physico-chemical parameters are indicated as mean ± SE during the whole exposure period. Letters indicate statistical differences according to the Scheirer-Ray-Hare test followed by Wilcoxon-Mann-Whitney tests when differences were detected (p value ≤ 0.05). Micropollutants concentrations are indicated as minimum and maximum values measured (the number of measures is indicated within brackets).

268

## 3.3. Effects on the Condition Index

At the end of the caging period, means ± SE in reference, WWTP exposure and urban runoff exposure were 51.67 ± 2.09, 53.65 ± 1.14 and 55.37 ± 2.65, respectively, and no significant difference was detected (ANOVA, p-value = 0.46). The significance of other factors such as natural variability, age of the individuals, diet, mixture of pollutants or hydrology is difficult to predict (Petrie et al., 2015) and would require a substantial increase in sampling and analytical effort. Aquatic toxicity is often undertaken under controlled laboratory conditions over short time period / acute exposure to pollutants. Field experiments, on the contrary, allow to determine the impact of mixtures and chronic exposure, hence providing more environmentally realistic, though less reproducible, information.

## 3.4. Gene transcription analyses

A total of 3,181 Differentially Expressed Genes (DEGs, FDR ≤ 0.05) were identified in WWTP exposure and urban runoff exposure compared to the reference site. In details, 1,447 and 412 were up-regulated and 1,040 and 623 were down-regulated in WWTP exposure and urban runoff exposure, respectively (Figure 2 and Appendix C). To validate the sequencing data, the transcription of 15 genes that showed variations in their expression levels was measured by the RT-qPCR method. RNAseq and RT-qPCR gave consistent results with a R² values of 0.91 and 0.82 for WWTP exposure vs. Reference and urban runoff exposure vs. Reference comparisons, respectively (Appendix A).

286
<figure><img src="image_3.png" /><figcaption>Figure 2: Venn diagram and number of differentially expressed genes. Number of up- (red) and down (green) regulated genes in individuals caged downstream of the WWTP and downstream of the urban runoff effluent. The number of contigs with a BLAST hit (1e-5) is indicated between brackets.</figcaption></figure>

There was little overlap amongst the transcripts differentially expressed in the two sites compared to the reference site. Fold change values of the regulated genes were relatively small, ranging from 0.51 to 2.38 in WWTP exposure and from 0.56 to 2.20 in urban runoff exposure (Table 3 and Appendix C).

The number of DEG is sometimes viewed as a proxy of the level of stress experienced by the organisms. However, the nature and the number of affected biological processes, as well as the magnitude of changes certainly provide more relevant information regarding the nature of the stress and the associated biological consequences of water quality. Such detailed information regarding the biological processes in which each DEG was involved were provided by their functional annotation, i.e. GO (Table 4) and KEGG orthology (Table 5).
<table><thead><tr><th>Contig name</th><th>Description</th><th>FC WWTP</th><th>FC runoff</th><th>Putative function</th></tr></thead><tbody><tr><td colspan="5"><strong>Xenobiotic detoxification</strong></td></tr><tr><td>Reg_ADHX.2.2</td><td>Alcohol dehydrogenase class III</td><td>1.51</td><td>1.76</td><td>Formaldehyde detoxification</td></tr><tr><td>Reg LOC108269088.1.3</td><td>Alcohol dehydrogenase class III</td><td>1.77</td><td>1.65</td><td>Formaldehyde detoxification</td></tr><tr><td>Reg LOC108269088.2.3</td><td>Alcohol dehydrogenase class III</td><td>1.85</td><td>1.55</td><td>Formaldehyde detoxification</td></tr><tr><td>Reg LOC108269088.3.3</td><td>Alcohol dehydrogenase class III</td><td>1.93</td><td>1.60</td><td>Formaldehyde detoxification</td></tr><tr><td>Reg_ADHX.1.2</td><td>Alcohol dehydrogenase class III</td><td>---</td><td>1.75</td><td>Formaldehyde detoxification</td></tr><tr><td>Reg LOC106174241.2.2</td><td>Alcohol dehydrogenase class III-like</td><td>---</td><td>1.69</td><td>Formaldehyde detoxification</td></tr><tr><td colspan="5"><strong>Energy metabolism</strong></td></tr><tr><td>Reg LOC105326051.4.4</td><td>Phosphoenolpyruvate carboxykinase, cytosolic [GTP]-like</td><td>1.83</td><td>1.69</td><td>Gluconeogenesis</td></tr><tr><td>Reg LOC105326362</td><td>Pyruvate carboxylase, mitochondrial-like</td><td>1.91</td><td>1.36</td><td>Gluconeogenesis</td></tr><tr><td>Reg LOC105326051.3.4</td><td>Phosphoenolpyruvate carboxykinase, cytosolic [GTP]-like</td><td>1.69</td><td>1.59</td><td>Gluconeogenesis</td></tr><tr><td>Reg LOC105347001</td><td>Neutral cholesterol ester hydrolase 1-like</td><td>---</td><td>0.63</td><td>Lipid catabolism</td></tr><tr><td>Reg LOC105344840.2.2</td><td>Complement C1q tumor necrosis factor-related protein 4-like</td><td>2.38</td><td>---</td><td>Regulator of glucose and lipid metabolism</td></tr><tr><td colspan="5"><strong>Immune response</strong></td></tr><tr><td>Reg LOC102197942</td><td>C-type lectin domain family 4 member E-like</td><td>0.54</td><td>---</td><td>Pathogen recognition</td></tr><tr><td>Reg_Fcer2.1.3</td><td>Low affinity immunoglobulin epsilon Fc receptor</td><td>0.58</td><td>---</td><td>Immune response</td></tr><tr><td>Reg_SAA2.2.2</td><td>Serum amyloid A-2 protein [Homo sapiens]</td><td>1.84</td><td>1.60</td><td>Response to inflammation</td></tr><tr><td>Reg LOC105347518</td><td>Interferon-induced protein 44-like</td><td>1.89</td><td>---</td><td>Response to virus</td></tr><tr><td>Reg ifi44</td><td>Interferon-induced protein 44 [Mus musculus]</td><td>1.89</td><td>---</td><td>Response to virus</td></tr><tr><td colspan="5"><strong>Cell signaling</strong></td></tr><tr><td>Reg LOC106880299.2.2</td><td>Serine/threonine-protein kinase WNK1-like</td><td>---</td><td>0.56</td><td>Regulation of cell signaling, survival, and proliferation</td></tr><tr><td>Reg LOC105321825</td><td>Probable G-protein coupled receptor 139</td><td>0.58</td><td>0.59</td><td>Signal transduction</td></tr><tr><td>Reg LOC105338118</td><td>Ski oncogene-like</td><td>---</td><td>0.62</td><td>Repressor of TGF-beta signaling</td></tr><tr><td>Reg_contig_05776</td><td>Receptor-type tyrosine-protein phosphatase T</td><td>---</td><td>1.63</td><td>Signal transduction</td></tr><tr><td colspan="5"><strong>Ion transport</strong></td></tr><tr><td>Reg LOC101856164.2.2</td><td>Solute carrier family 23 member 1-like</td><td>0.56</td><td>---</td><td>Sodium/ascorbate cotransporter</td></tr><tr><td>Reg LOC105334363</td><td>Monocarboxylate transporter 9-like</td><td>0.57</td><td>---</td><td>Transport</td></tr><tr><td>Reg_contig_33484</td><td>Zinc transporter ZIP1</td><td>0.60</td><td>0.72</td><td>Transport, endogenous zinc uptake</td></tr><tr><td>Reg LOC106156136</td><td>Stromal interaction molecule homolog</td><td>0.72</td><td>0.62</td><td>Calcium transport</td></tr><tr><td colspan="5"><strong>Miscellaneous</strong></td></tr><tr><td>Reg LOC105335663.1.2</td><td>Lysosomal-trafficking regulator-like</td><td>---</td><td>0.62</td><td>Intracellular protein trafficking</td></tr><tr><td>Reg LOC105338935</td><td>Homeodomain-interacting protein kinase 2-like</td><td>---</td><td>0.62</td><td>Cellular response to hypoxia</td></tr><tr><td>Reg LOC105345002.3.3</td><td>Cat eye syndrome critical region protein 5-like</td><td>0.51</td><td>0.60</td><td>Chromatin remodelling</td></tr><tr><td>Reg LOC101857055</td><td>Rab5 GDP/GTP exchange factor-like</td><td>0.53</td><td>---</td><td>Membrane trafficking of recycling endosomes</td></tr><tr><td>Reg LOC105340685</td><td>Neuronal PAS domain-containing protein 4-like</td><td>0.54</td><td>---</td><td>Neuron plasticity</td></tr><tr><td>Reg LOC105330779.3.4</td><td>EMILIN-2-like</td><td>0.59</td><td>---</td><td>Cell adhesion</td></tr><tr><td>Reg_VKT2</td><td>Kunitz serine protease inhibitor Pr-mulgin 2</td><td>0.64</td><td>0.62</td><td>Endopeptidase inhibitor</td></tr><tr><td>Reg LOC105319937</td><td>Beta-TrCP-like</td><td>0.69</td><td>0.62</td><td>Regulation of protein translation, cell growth and survival</td></tr><tr><td>Reg LOC100633144</td><td>Zinc finger BED domain-containing protein 1-like</td><td>1.77</td><td>2.20</td><td>Transcription factor</td></tr><tr><td>Reg LOC105317854.5.5</td><td>Baculoviral IAP repeat-containing protein 7-A-like</td><td>1.90</td><td>---</td><td>Apoptosis regulation</td></tr><tr><td>Reg LOC105348304.2.2</td><td>Heat shock protein 70 B2-like</td><td>1.98</td><td>---</td><td>Protein folding</td></tr></tbody></table>
302 Table 3: List of the most up- and down-regulated genes in WWTP and runoff compared to reference.
303 The 10 most up-regulated (indicated in red) and the 10 most down-regulated (indicated in green) were
304 identified in each site compared to the reference site. (FC, fold-change)

305

306

### 3.5. Common transcriptomic response of *C. fluminea* from both sites

307 Genes involved in xenobiotic detoxification were amongst the most up-regulated genes (Table 3),
308 particularly many isoforms of class III alcohol dehydrogenases (ADH) that are responsible for
309 formaldehyde detoxification by metabolizing glutathione adducts (Haseba et al., 2006). This enzyme
310 involved in Phase 1 xenobiotics metabolism was previously shown to be up-regulated in the digestive
311 gland of bivalves exposed to various organic pollutants, particularly pharmaceuticals (Bebianno et al.,
312 2016; Schmidt et al., 2014). These genes also participate to the KEGG pathway ‘Degradation of
313 aromatic compounds’ in Table 4, suggesting the presence of organic pollutants in WWTP and runoff
314 waters. This hypothesis is partly confirmed by higher concentrations of some pesticides presented in
315 Table 2. Repressed contigs were linked to cell signalling and ion transport. These results were
316 confirmed by GO (Table 4) and KEGG (Table 5) enrichment analysis. Details are given in Table 3 and
317 Appendix B. The mTOR (mammalian target of rapamycin), the FoxO (forkhead box O) and the
318 AMPK (5'AMP-activated protein kinase) pathways (Table 5) are highly conserved across metazoans
319 and intimately cross-talk in response to metabolic stress in order to maintain cell homeostasis. Under
320 stress, it is expected that the mTOR pathway is inhibited leading to a reduction in protein and lipid
321 biosynthesis and an increased rate of autophagy whereas FoxO and AMPK pathways are activated to
322 promote glucose / fatty acid metabolism and autophagy too (Hay, 2011; Mendoza et al., 2011; Zhao et
323 al., 2017). This metabolic dialogue aims to reduce cell growth and proliferation and to sustain
324 anabolic processes such as energy production by providing fatty acid and amino acids to the
325 mitochondria (Laplante and Sabatini, 2009). In a stressful environment, these processes might be used
326 in the digestive gland as an alternative source of energy production by mitochondrial oxidation in
327 order to maintain ATP/ADP homeostasis and the energy demand of other organs. Here, we observed a
328 higher expression of genes involved in energy metabolism (glycolysis, gluconeogenesis and fatty acid
329 degradation) in clams exposed to either WWTP- or urban runoff-impacted waters, which suggests the
330 activation of FoxO (and a repression of mTOR) (Zhao et al., 2017). Yet, genes involved in all the
331 aforementioned pathways were downregulated. We can hypothesize that these signaling disorders and
332 metabolic changes might have been caused directly by the presence of organic pollutants (Regoli and
333 Giuliani, 2014). However, other parameters might affect energy metabolism of which adjustment
334 plays a central role in adaptation and tolerance to environmental stress in aquatic invertebrates
335 (Sokolova et al., 2012). Cellular hypoxia could inhibit mTOR (Laplante and Sabatini, 2009), but
336 appears unlikely as it was evidenced that even under reduced environmental O₂ concentration as
337 observed in WWTP and runoff effluents, *C. fluminea* is able to maintain a remarkably stable level of
338 oxygen in its internal milieu without change in its ventilatory activity (Tran et al., 2000). In the same
339 manner, starvation could be ruled out as stressor since we observed no difference amongst CI after 3
340 months compared to our reference site. Moreover, while food deprivation would inhibit mTOR, it
341 would activate foxO and AMPK, which was not observed here at the transcriptional level.
342 Temperature is also a crucial variable that determines metabolic rate of poikilothermic organisms but
343 as stated previously, temperature differences amongst sites were small (Table 2).

344 **3.6. Transcriptomic response to WWTP effluent**

345 Amongst the most responsive genes in WWTP exposure, it is particularly relevant in stressful
346 environment to find homologues of genes involved in the inhibition of apoptosis, such as BIRC7
347 (Baculoviral IAP repeat-containing protein 7 ; Gyrd-Hansen and Meier, 2010), and protein folding,
348 such as a Heat Shock Protein (HSP) 70 isoform (Table 3). Members of the 70 kDa chaperone family
349 are found in nearly all subcellular compartments of nucleated cells. HSP70 rapidly accumulates after
350 exposure to environmental stress and high cytosolic levels are known to protect from apoptotic cell
351 death (Hartl and Hayer-Hartl, 2002). We also found contigs related to immune response, inflammation
352 and response to pathogens that were either up- or down-regulated almost 2 fold. A potential impact of
353 WWTP exposure on protein metabolism was supported by the significant enrichment of the GO term
354 "protein folding" carried out by 21 up-regulated genes (Table 4) and the KEGG category "protein
355 processing in endoplasmic reticulum (ER)" represented by 26 genes (Table 5). This up-regulation of
356 genes involved in protein processing and folding was further coupled with a higher level of RNA
357 processing (GO:006396, Table 4) and aminoacyl-tRNA biosynthesis (i.e. bounding of a transfer RNA
with its corresponding amino acid, Tables 4 and 5). The functional analysis of DEG also suggested a reduction in auto- and mitophagy (Table 5). Protein synthesis and reduced autophagy generally result from the activation of the mTOR signaling pathway (Mendoza et al., 2011) while, as stated in the previous chapter, the maintenance of cellular homeostasis requires FoxO activation and mTOR inhibition (Zhao et al., 2017). That reinforces the hypothesis of severe metabolic disorders and/or an impairment in signal transduction from the environment to the gene machinery. These pathways were identified as a potential marker for chronic environmental challenge in the Pacific oyster, *Crassostrea gigas*, indicating tension between cell renewal and apoptosis (Clark et al., 2013).

The exposure to WWTP effluent could also be responsible for a higher level of misfolded proteins as suggested by the over-representation of genes involved in proteolysis (GO:0051603, Table 4) amongst DEGs. A comprehensive look at the KEGG orthology revealed that some of these genes were involved in proteasome formation (17 up-regulated genes) and that part of the ER-associated genes were responsible for ER-associated degradation (ERAD, Appendix D). Protein degradation might also supply amino acids for gluconeogenesis and be an evidence of severe energy need / deficit in the organisms exposed to WWTP waters. This impairment of protein synthesis, which can cause cell malfunction, deterioration or even cell death, might also result from modification of RNA and disturbance of the translational process as suggested by the stimulation of genes involved in mRNA surveillance and RNA degradation (Table 5).

Similar to what we found above, genes involved in energy metabolism (8 and 6 DEGs respectively linked to fatty acid metabolism and pyruvate metabolism, Table 5) showed an increased level of transcription. However, WWTP is the only site to exhibit signs of molecular and cellular damages ultimately impacting longevity. This confirms a higher level of stress as previously suggested by the greater number of DEGs.
<table><thead><tr><th></th><th>GO ID</th><th>Description</th><th>FDR</th><th>n DEG</th></tr></thead><tbody><tr><td>common to both sites</td><td>GO:0006094</td><td>Gluconeogenesis</td><td>▲ 3.45E-2</td><td>3</td></tr><tr><td rowspan="10">Enriched in WWTP</td><td>GO:0006457</td><td>Protein folding</td><td>▲ 2.19E-11</td><td>21</td></tr><tr><td>GO:0051603</td><td>Proteolysis involved in cellular protein catabolism</td><td>▲ 5.59E-4</td><td>12</td></tr><tr><td>GO:0006487</td><td>Protein N-linked glycosylation</td><td>▲ 4.60E-2</td><td>2</td></tr><tr><td>GO:0043039</td><td>tRNA aminoacylation</td><td>▲ 1.55E-6</td><td>12</td></tr><tr><td>GO:0006396</td><td>RNA processing</td><td>▲ 5.59E-4</td><td>17</td></tr><tr><td>GO:0016071</td><td>mRNA metabolic process</td><td>▲ 3.60E-2</td><td>7</td></tr><tr><td>GO:0006779</td><td>Porphyrin-containing compound biosynthesis</td><td>▲ 4.60E-2</td><td>3</td></tr><tr><td>GO:0006637</td><td>Acyl-CoA metabolic process</td><td>▲ 2.63E-2</td><td>4</td></tr><tr><td>GO:0045454</td><td>Cell redox homeostasis</td><td>▲ 2.63E-2</td><td>7</td></tr><tr><td rowspan="2">Enriched in runoff</td><td>GO:0007165</td><td>Signal transduction</td><td>▼ 1.80E-5</td><td>58</td></tr><tr><td>GO:0010468</td><td>Regulation of gene expression</td><td>▼ 4.15E-2</td><td>24</td></tr><tr><td></td><td>GO:0006468</td><td>Protein phosphorylation</td><td>▼ 4.17E-04</td><td>17</td></tr><tr><td></td><td>GO:0035556</td><td>Intracellular signal transduction</td><td>▼ 6.93E-7</td><td>28</td></tr></tbody></table>

381

382 Table 4: List of significantly enriched Gene Ontology terms. (n DEG, number of associated differentially expressed genes; FDR, false discovery rate Fisher

383

exact test). Red and green arrows indicate up-regulated and down-regulated pathways, respectively.
<table><thead><tr><th rowspan="2"></th><th rowspan="2">KEGG pathway</th><th rowspan="2">Description</th><th colspan="3">WWTP</th><th colspan="3">Runoff</th><th rowspan="2">function</th></tr><tr><th>pvalue</th><th>up</th><th>down</th><th>pvalue</th><th>up</th><th>down</th></tr></thead><tbody><tr><td rowspan="2">Common to both sites</td><td>map00071</td><td>Fatty acid degradation</td><td>▲</td><td>1.53E-3</td><td>11</td><td>-</td><td>3.38E-2</td><td>6</td><td>2</td><td>energy metabolism</td></tr><tr><td>map00010</td><td>Glycolysis / Gluconeogenesis</td><td>▲</td><td>1.33E-2</td><td>9</td><td>1</td><td>2.80E-2</td><td>9</td><td>-</td><td>detoxification</td></tr><tr><td></td><td>map01220</td><td>Degradation of aromatic compounds</td><td>▲</td><td>1.87E-2</td><td>4</td><td>-</td><td>3.09E-4</td><td>6</td><td>-</td><td>detoxification</td></tr><tr><td rowspan="10">Enriched in WWTP</td><td>map04150</td><td>mTOR signaling pathway</td><td>▼</td><td>1.90E-4</td><td>5</td><td>13</td><td>3.62E-4</td><td>-</td><td>17</td><td>signaling pathways</td></tr><tr><td>map04068</td><td>FoxO signaling pathway</td><td>▼</td><td>5.84E-4</td><td>4</td><td>12</td><td>3.33E-3</td><td>3</td><td>11</td><td>signaling pathways</td></tr><tr><td>map04072</td><td>Phospholipase D signaling pathway</td><td>▼</td><td>1.35E-2</td><td>2</td><td>11</td><td>8.88E-7</td><td>-</td><td>21</td><td>signaling pathways</td></tr><tr><td>map04152</td><td>AMPK signaling pathway</td><td>▼</td><td>2.41E-2</td><td>5</td><td>7</td><td>3.50E-4</td><td>3</td><td>13</td><td>signaling pathways</td></tr><tr><td>map04071</td><td>Sphingolipid signaling pathway</td><td>▼</td><td>2.86E-2</td><td>-</td><td>11</td><td>1.20E-3</td><td>-</td><td>14</td><td>signaling pathways</td></tr><tr><td>map04012</td><td>ErbB signaling pathway</td><td>▼</td><td>9.35E-3</td><td>1</td><td>8</td><td>5.93E-6</td><td>-</td><td>14</td><td>signaling pathways</td></tr><tr><td>map01521</td><td>EGFR tyrosine kinase inhibitor resistance</td><td>▼</td><td>2.70E-2</td><td>1</td><td>7</td><td>5.91E-4</td><td>-</td><td>11</td><td>signaling pathways</td></tr><tr><td>map04630</td><td>Jak-STAT signaling pathway</td><td>▼</td><td>2.75E-2</td><td>2</td><td>4</td><td>2.36E-2</td><td>-</td><td>6</td><td>signaling pathways</td></tr><tr><td>map03040</td><td>Spliceosome</td><td>▲</td><td>3.11E-9</td><td>33</td><td>-</td><td>NS</td><td>-</td><td>-</td><td>gene expression</td></tr><tr><td>map00970</td><td>Aminoacyl-tRNA biosynthesis</td><td>▲</td><td>1.13E-5</td><td>14</td><td>-</td><td>NS</td><td>-</td><td>-</td><td>gene expression</td></tr><tr><td>map03015</td><td>mRNA surveillance pathway</td><td>▲</td><td>3.75E-3</td><td>11</td><td>2</td><td>NS</td><td>-</td><td>-</td><td>gene expression</td></tr><tr><td>map03018</td><td>RNA degradation</td><td>▲</td><td>3.10E-2</td><td>8</td><td>3</td><td>NS</td><td>-</td><td>-</td><td>gene expression</td></tr><tr><td>map04141</td><td>Protein processing in endoplasmic reticulum</td><td>▲</td><td>2.60E-6</td><td>26</td><td>2</td><td>NS</td><td>-</td><td>-</td><td>protein processing</td></tr><tr><td>map03050</td><td>Proteasome</td><td>▲</td><td>5.55E-12</td><td>17</td><td>-</td><td>NS</td><td>-</td><td>-</td><td>protein processing</td></tr><tr><td>map00510</td><td>N-Glycan biosynthesis</td><td>▲</td><td>9.35E-3</td><td>8</td><td>1</td><td>NS</td><td>-</td><td>-</td><td>protein processing</td></tr><tr><td>map01212</td><td>Fatty acid metabolism</td><td>▲</td><td>3.32E-2</td><td>8</td><td>-</td><td>NS</td><td>-</td><td>-</td><td>energy metabolism</td></tr><tr><td>map00620</td><td>Pyruvate metabolism</td><td>▲</td><td>4.32E-2</td><td>6</td><td>1</td><td>NS</td><td>-</td><td>-</td><td>energy metabolism</td></tr><tr><td>map04140</td><td>Autophagy - animal</td><td>▼</td><td>1.81E-2</td><td>3</td><td>13</td><td>NS</td><td>-</td><td>-</td><td>Cell death</td></tr><tr><td>map04137</td><td>Mitophagy - animal</td><td>▼</td><td>9.35E-3</td><td>2</td><td>7</td><td>NS</td><td>-</td><td>-</td><td>Cell death</td></tr><tr><td>map04211</td><td>Longevity regulating pathway</td><td>▼</td><td>2.57E-2</td><td>2</td><td>7</td><td>NS</td><td>-</td><td>-</td><td>Cell death</td></tr><tr><td rowspan="10">Enriched in runoff</td><td>map04010</td><td>MAPK signaling pathway</td><td>▼</td><td>NS</td><td>-</td><td>-</td><td>8.27E-4</td><td>-</td><td>23</td><td>signaling pathways</td></tr><tr><td>map04070</td><td>Phosphatidylinositol signaling system</td><td>▼</td><td>NS</td><td>-</td><td>-</td><td>3.68E-2</td><td>-</td><td>14</td><td>signaling pathways</td></tr><tr><td>map04144</td><td>Endocytosis</td><td>▼</td><td>NS</td><td>-</td><td>-</td><td>1.83E-2</td><td>-</td><td>22</td><td>signaling pathways</td></tr><tr><td>map04062</td><td>Chemokine signaling pathway</td><td>▼</td><td>NS</td><td>-</td><td>-</td><td>5.41E-6</td><td>-</td><td>18</td><td>signaling pathways</td></tr><tr><td>map04666</td><td>Fc gamma R-mediated phagocytosis</td><td>▼</td><td>NS</td><td>-</td><td>-</td><td>5.00E-4</td><td>-</td><td>15</td><td>signaling pathways</td></tr><tr><td>map04664</td><td>Fc epsilon RI signaling pathway</td><td>▼</td><td>NS</td><td>-</td><td>-</td><td>2.93E-5</td><td>-</td><td>12</td><td>signaling pathways</td></tr><tr><td>map04920</td><td>Adipocytokine signaling pathway</td><td>▼</td><td>NS</td><td>-</td><td>-</td><td>4.31E-3</td><td>3</td><td>6</td><td>signaling pathways</td></tr><tr><td>map04670</td><td>Leukocyte transendothelial migration</td><td>▼</td><td>NS</td><td>-</td><td>-</td><td>4.67E-2</td><td>-</td><td>9</td><td>signaling pathways</td></tr><tr><td>map04912</td><td>GnRH signaling pathway</td><td>▼</td><td>NS</td><td>-</td><td>-</td><td>1.41E-2</td><td>-</td><td>14</td><td>signaling pathways</td></tr><tr><td>map04915</td><td>Estrogen signaling pathway</td><td>▼</td><td>NS</td><td>-</td><td>-</td><td>3.70E-2</td><td>-</td><td>13</td><td>endocrine system</td></tr><tr><td>map03320</td><td>PPAR signaling pathway</td><td>▼</td><td>NS</td><td>-</td><td>-</td><td>3.38E-2</td><td>3</td><td>5</td><td>endocrine system</td></tr><tr><td>map04540</td><td>Gap Junction</td><td>▼</td><td>NS</td><td>-</td><td>-</td><td>1.16E-3</td><td>-</td><td>10</td><td>Misc.</td></tr><tr><td>map00562</td><td>Inositol phosphate metabolism</td><td>▼</td><td>NS</td><td>-</td><td>-</td><td>1.23E-2</td><td>-</td><td>9</td><td>Misc.</td></tr></tbody></table>
Table 5: List of significantly enriched metabolic pathways based on KEGG orthology. Red and green arrows indicate globally up-regulated and down-regulated pathways, respectively. Number of up-regulated and down-regulated contigs involved in each pathway are indicated in red and green, respectively. Enrichment was tested with a Fisher exact test (pvalue ≤ 0.05). Only pathways represented by >5 contigs were considered.

### 3.7. Transcriptomic response to urban runoff

The down-regulation of a homeodomain-interacting protein kinase 2 (HIPK2) homologue (Table 3) is particularly relevant in the most hypoxic site (Table 2). HIPK2 was found to repress the transcription of Hypoxia Inducible Factor (HIF) 1α (Nardinocchi et al., 2009). As a transcription factor, HIF1 targets genes involved in glucose metabolism and cell survival in response to low O₂ (Ke and Costa, 2006). This cellular response to hypoxia might lead *C. fluminea* to maintain its homeostasis and ventilatory activity (Tran et al., 2000).

In urban runoff exposure, the repression of DEGs involved in signalling pathways (GO:0035556, Table 4) was also observed. Thanks to KEGG orthology, these pathways could be grouped into immunity and endocrine system related processes (Table 5 and Appendices B and C). The MAPK (Mitogen-activated protein kinases) pathway is amongst the most ancient signal transduction pathways and is widely used throughout evolution in many physiological processes (Cargnello and Roux, 2011; Morrison, 2012). The classical MAPK pathway (also referred to as the ERK1/2 module) is part of several aforementioned pathways and contigs encoding proteins involved in this module were down-regulated at both sites. However, in clams exposed to urban runoff, genes involved in the JNK and p38 modules were also downregulated (Appendix E). These modules are activated in response to various cellular stresses, including heat shock, oxidative stress, cytokines, hypoxia, interleukin-1 (IL-1), lipopolysaccharide (LPS) and growth factor deprivation (Bogoyevitch et al., 2010; Cuadrado and Nebreda, 2010). While the JNK plays an important role in the apoptotic response to cellular stresses (Dhanasekaran and Reddy, 2008), the p38 plays a critical role in normal immune and inflammatory responses. A major function of the p38 pathway is to produce proinflammatory cytokines such as TNF-α, IL-1, IL-6 and IL-12, which are essential for the eradication of infectious microorganisms
413 (Arthur and Ley, 2013; Cuadrado and Nebreda, 2010; Zhang and Dong, 2005). In addition, the
414 ERK1/2 pathway has an important role in macrophages, regulating cytokine production (Arthur and
415 Ley, 2013). Hence, the inhibition of MAPK pathways could result in severe reduction of the immune
416 response in organisms from urban runoff exposure unlike WWTP exposure where the ability to
417 respond to virus was suggested by the presence of immune-related contigs amongst the most up-
418 regulated (Table 3). This inhibition might explain the observed reduction of several immune processes
419 in the present study like endocytosis, phagocytosis and response to and production of chemokines. In
420 bivalves, cellular response is carried out by circulating hemocytes (Allam and Pales Espinosa, 2016;
421 Song et al., 2010). Such effects are similar to what was observed in fish for instance (Segner et al.,
422 2011). In bivalves, different classes of environmental pollutants like metals, pharmaceuticals, biocides,
423 polychlorinated biphenyls (PCBs) and PAHs were shown to affect immunity (Galloway and Depledge,
424 2001). The hypothesis that urban runoff impacts bivalve immunity should be confirmed in future
425 studies by the measurement of hemocytes-related variables such as type counts, size, and production
426 of reactive oxygen species (e.g. Lassudrie et al., 2016). This phenomenon is of particular interest
427 since, in addition to metal and organic pollutants, pathogens can find their way into surface water
428 through leaking sewer systems, sewer pumping station overflows, seepage from septic systems,
429 agricultural runoff and discharge of treated wastewater into aquatic environments. High numbers of
430 enteric viruses, bacteria and protozoa have been reported in stormwater runoff (Sidhu et al., 2012).
431 Such contamination is hence more likely to occur on urban runoff site following heavy rainfall than on
432 WWTP where no leakage was noticed over the studied period. It is interesting to note that MAPKs are
433 targeted by pharmaceuticals such as anti-cancer and anti-rheumatic drugs, and by some pathogens as
434 an immune evasion strategy (Arthur and Ley, 2013). This is the case, for instance, for Vibrio
435 parahaemolyticus and Salmonella enterica, both responsible for gastroenteritis, that is a particularly
436 prevalent disease during the study period, i.e. winter months. Immune system could also be altered by
437 both generalized stress and changes in PPAR pathway (DeWitt et al., 2009). Peroxisome proliferator-
438 activated receptors (PPARs) regulate inflammation, lipid metabolism and energy utilization. Changes
439 in PPAR pathway might result from endocrine disrupting chemicals (Casals-Casas et al., 2008). In
440 fish, such changes were documented in organisms exposed to anthropogenic pollutants like organotins,
oil and nonylphenol (Bilbao et al., 2010; Cocci et al., 2013; Pavlikova et al., 2010). Other disturbances affected GnRH and estrogen signaling pathways. Numerous chemicals, including pharmaceuticals, pesticides, plastics feedstocks and additives, and phytochemicals, can act as estrogen agonists or antagonists, mimicking or blocking the natural effects of estrogens in invertebrates (Keay and Thornton, 2009). Though WWTP are important sources, estrogens and other aforementioned endocrine disruptors can also be found in runoff and surface waters (Adeel et al., 2017). Here, we were only able to measure a limited number of pesticides including glyphosate, AMPA and terbutryn that showed higher concentrations in runoff. Glyphosate is currently the most widely used herbicide in the world. In fish, that received most of the interest amongst aquatic organisms, glyphosate showed potential effects on reproduction through the disruption of steroid hormone synthesis and increased oxidative stress (e.g. Uren Webster et al., 2014). Bivalves received little interest, however recent studies in *Mytilus galloprovincialis* evidenced effects of glyphosate on hemocytes, membrane transport, energy metabolism and neurotransmission (Matozzo et al., 2018; Milan et al., 2018). A more comprehensive chemical analysis of the effluent water looking at various class of compounds would be necessary to confirm the biological outcomes on bivalve's physiology.

### 3.8. Insights into the nature of environmental stress

Taken together, our results suggest that exposure to organic contaminants and hypoxia could lead to energy impairments, cellular damages, immune disorders and endocrine disruption through disturbances in signaling pathways. We observed no difference in CI of organisms from the three sites at the end of the exposure period. In bivalves, it has been shown that environmental stress may affect physiology without impacting CI (e.g. Riascos et al., 2012), then the changes in transcription we observed might represent early warning markers occurring before observable deleterious consequences at the whole-organism level (Forbes et al., 2006).

These outcomes at the gene expression level are compatible with pollution generally associated with urban waters. The presence of emerging contaminants in the environment is mainly attributed to the discharge of treated wastewater from WWTP. Highly prescribed therapeutic substances like anti-inflammatory drugs, β-blockers, antidepressants and antiepileptics, together with personal care products, are the most ubiquitous chemicals to influent wastewaters due to their incomplete removal
during wastewater treatment (Muñoz et al., 2008; Petrie et al., 2015). Urban runoff is a relatively recent concern. It carries multiple pollutants such as metals (mainly Cd, Pb, Zn and Cu), pesticides, PCBs, PAHs, organic material (pet waste, leaves, litter, etc) of which the decay consume oxygen, and bacteria from sewage, pets and urban wildlife (Dong and Lee, 2009; Huber et al., 2016; McLellan et al., 2015; Revitt et al., 2014). Despite being considered as a major source (Heberer, 2002), effects of pharmaceuticals were not observed downstream of the WWTP. We can hypothesize that the effects of the pollutants appeared or were reinforced by additional input in urban runoff. Indeed, studies have shown that mixtures of chemicals exhibit greater effects than the pollutants individually (Cedergreen, 2014). However, cautions must be taken in the prediction and interpretation of biological consequences of the stress associated with both sources. Variability of stormwater quality is high and depends on the accumulation time on diverse impervious surface between two precipitation events (Zgheib et al., 2011) whereas treated wastewater flux and composition are more constant and predictable over a three-month period.

Solely based on a transcriptomic approach and a limited number of chemicals, our study still provided some valuable first clues regarding the nature of the stress experienced by *C. fluminea* and the identification of initiating molecular events of toxicological outcomes. Nevertheless, the complexity of metabolic pathways and their interplay (illustrated by the many signaling pathways cited previously), as well as the diversity and variability in the chemical composition of the two types of effluents make difficult the identification of a specific initiating event and the prediction of the adverse biological consequences. To counter this limitation, the chemical characterization of the pollution, as well as the presence of pathogens appear necessary for future studies in order to confirm or clarify transcriptomic results. In the same manner, the integration of multi-omics technologies, such as the coupling of transcriptomic and metabolomic information (Santos et al., 2010), has the potential to better predict toxicological outcomes and gain further insight into the mechanistic aspects of chemical effects. For instance, such approach could confirm whether transcriptional changes translate into metabolic or physiological adjustments and/or impacts. Freshwater bivalves are amongst the most threatened molluscs worldwide (Lydeard et al., 2004). Such an integrative approach could then be applied to assess the impacts of polluted water in endangered native species.
497

498

## 4. Conclusions

499 Our results provide insights in the molecular impact of low quality waters in Asian clams, *C. fluminea*.

500 The chemical characterization of WWTP and urban runoff exposures over a 3-month period is hardly achievable in field experiments. However, RNA sequencing could be useful for diagnosing sites where
501 the level and nature of contamination are unknown as it allows the investigation of a large collection
502 of biological processes without *a priori*. Our results revealed the up-regulation of genes involved in
detoxification, suggesting the presence of organic pollutants on both sites. The down-regulation of
503 genes involved in intracellular signaling pathways suggested metabolic disorders likely to impact
504 energy metabolism, protein metabolism and cell survival in response to WWTP effluent, and
505 endocrine and immune systems in response to urban runoff exposure. These findings could then allow
506 the identification of relevant biological parameters to survey in similar future studies, such as energy
507 metabolism, endocrine disorders or immunity that should be performed at a higher level of integration,
508 such as histological or enzymatic assays. Moreover, this work highlights the necessity of a
509 comprehensive chemical characterization of the pollution in order to better predict the toxicological
510 outcomes, particularly if applied to endangered freshwater bivalve species.
511
512

512
513 Acknowledgements

514 Authors thank all the partners and people involved in the project. We thank Bordeaux Metropole for
515 allowing the access to sampling sites. This work was supported by the Bordeaux Metropole’s research
516 program REGARD funded by the “Agence Française pour la Biodiversité” and the “Agence de l’Eau
517 Adour Garonne”. Thanks to the anonymous reviewers for their useful comments and the improvement
518 of our manuscript.

519

520 Competing Interests Statement

521 The authors declare no competing interests.

522

523 Appendices

524 Appendix A: Summary of quantitative real-time PCR validation.

525 Appendix B: Complete annotation of the *de-novo* transcriptome of *C. fluminea*'s digestive gland.

526 Appendix C: Results of DESeq2 differential expression analysis.

527 Appendix D: Details of the KEGG pathway 04141.

528 Appendix E: Details of the KEGG pathway 04010.
529 References:

530 Adeel, M., Song, X., Wang, Y., Francis, D., Yang, Y., 2017. Environmental impact of estrogens on human, animal and plant life: A critical review. Environ. Int. 99, 107–119.

531 https://doi.org/10.1016/j.envint.2016.12.010

533 Alexa, A., Rahnenfuhrer, J., 2016. topGO: Enrichment Analysis for Gene Ontology. R Package Version 2280.

535 Allam, B., Pales Espinosa, E., 2016. Bivalve immunity and response to infections: Are we looking at the right place? Fish Shellfish Immunol., Special Issue: ISFSI 2016 53, 4–12.

537 https://doi.org/10.1016/j.fsi.2016.03.037

538 Arini, A., Daffe, G., Gonzalez, P., Feurtet-Mazel, A., Baudrimont, M., 2014. Detoxification and recovery capacities of Corbicula fluminea after an industrial metal contamination (Cd and Zn): a one-year depuration experiment. Environ. Pollut. Barking Essex 1987 192, 74–82.

541 https://doi.org/10.1016/j.envpol.2014.04.012

542 Arlos, M.J., Bragg, L.M., Parker, W.J., Servos, M.R., 2015. Distribution of selected antiandrogens and pharmaceuticals in a highly impacted watershed. Water Res., Occurrence, fate, removal and assessment of emerging contaminants in water in the water cycle (from wastewater to drinking water) 72, 40–50. https://doi.org/10.1016/j.watres.2014.11.008

546 Arthur, J.S.C., Ley, S.C., 2013. Mitogen-activated protein kinases in innate immunity. Nat. Rev. Immunol. 13, 679–692. https://doi.org/10.1038/nri3495

548 Baillon, L., Pierron, F., Coudret, R., Normendeau, E., Caron, A., Peluhet, L., Labadie, P., Budzinski, H., Durrieu, G., Sarraco, J., Elie, P., Couture, P., Baudrimont, M., Bernatchez, L., 2015.

550 Transcriptome profile analysis reveals specific signatures of pollutants in Atlantic eels.

551 Ecotoxicol. Lond. Engl. 24, 71–84. https://doi.org/10.1007/s10646-014-1356-x

552 Bebianno, M.J., Sroda, S., Gomes, T., Chan, P., Bonnafe, E., Budzinski, H., Geret, F., 2016. Proteomic changes in Corbicula fluminea exposed to wastewater from a psychiatric hospital. Environ.

554 Sci. Pollut. Res. 23, 5046–5055. https://doi.org/10.1007/s11356-015-5395-5

555 Belles, A., Pardon, P., Budzinski, H., 2014. Development of an adapted version of polar organic chemical integrative samplers (POCIS-Nylon). Anal. Bioanal. Chem. 406, 1099–1110.

557 https://doi.org/10.1007/s00216-013-7286-2

558 Bertucci, A., Pierron, F., Thébault, J., Klopp, C., Bellec, J., Gonzalez, P., Baudrimont, M., 2017. Transcriptomic responses of the endangered freshwater mussel Margaritifera margaritifera

560 to trace metal contamination in the Dronne River, France. Environ. Sci. Pollut. Res. 1–15.

561 https://doi.org/10.1007/s11356-017-0294-6

562 Bilbao, E., Raingeard, D., de Cerio, O.D., Ortiz-Zarragoitia, M., Ruiz, P., Izagirre, U., Orbea, A., Marigómez, I., Cajaraville, M.P., Cancio, I., 2010. Effects of exposure to Prestige-like heavy fuel oil and to perfluorooctane sulfonate on conventional biomarkers and target gene

565 transcription in the thicklip grey mullet Chelon labrosus. Aquat. Toxicol. Amst. Neth. 98, 282–296. https://doi.org/10.1016/j.aquatox.2010.02.018

567 Bogoyevitch, M.A., Ngoei, K.R.W., Zhao, T.T., Yeap, Y.Y.C., Ng, D.C.H., 2010. c-Jun N-terminal kinase (JNK) signaling: Recent advances and challenges. Biochim. Biophys. Acta BBA - Proteins

569 Proteomics, Inhibitors of Protein Kinases 1804, 463–475.

570 https://doi.org/10.1016/j.bbpap.2009.11.002

571 Cabau, C., Escudié, F., Djari, A., Guiguen, Y., Bobe, J., Klopp, C., 2017. Compacting and correcting Trinity and Oases RNA-Seq de novo assemblies. PeerJ 5, e2988.

573 https://doi.org/10.7717/peerj.2988

574 Cabral, J.P.S., 2010. Water Microbiology. Bacterial Pathogens and Water. Int. J. Environ. Res. Public. Health 7, 3657–3703. https://doi.org/10.3390/ijerph7103657

576 Camacho, C., Coulouris, G., Avagyan, V., Ma, N., Papadopoulos, J., Bealer, K., Madden, T.L., 2009. BLAST+: architecture and applications. BMC Bioinformatics 10, 421.

578 https://doi.org/10.1186/1471-2105-10-421
579 Cargnello, M., Roux, P.P., 2011. Activation and Function of the MAPKs and Their Substrates, the
580 MAPK-Activated Protein Kinases. Microbiol. Mol. Biol. Rev. MMBR 75, 50-83.

581 https://doi.org/10.1128/MMBR.00031-10

582 Casals-Casas, C., Feige, J.N., Desvergne, B., 2008. Interference of pollutants with PPARs: endocrine
583 disruption meets metabolism. Int. J. Obes. 32, S53-S61. https://doi.org/10.1038/ijo.2008.207

584 Cedergreen, N., 2014. Quantifying Synergy: A Systematic Review of Mixture Toxicity Studies within
585 Environmental Toxicology. PLoS ONE 9. https://doi.org/10.1371/journal.pone.0096580

586 Clark, M.S., Thorne, M.A.S., Amaral, A., Vieira, F., Batista, F.M., Reis, J., Power, D.M., 2013.
587 Identification of molecular and physiological responses to chronic environmental challenge in
588 an invasive species: the Pacific oyster, Crassostrea gigas. Ecol. Evol. 3, 3283-3297.
589 https://doi.org/10.1002/ece3.719

590 Cocci, P., Mosconi, G., Palermo, F.A., 2013. Effects of 4-nonylphenol on hepatic gene expression of
591 peroxisome proliferator-activated receptors and cytochrome P450 isoforms (CYP1A1 and
592 CYP3A4) in juvenile sole (Solea solea). Chemosphere 93, 1176-1181.

593 https://doi.org/10.1016/j.chemosphere.2013.06.058

594 Cuadrado, A., Nebreda, A.R., 2010. Mechanisms and functions of p38 MAPK signalling. Biochem. J.
595 429, 403-417. https://doi.org/10.1042/BJ20100323

596 DeWitt, J.C., Shnyra, A., Badr, M.Z., Loveless, S.E., Hoban, D., Frame, S.R., Cunard, R., Anderson, S.E.,
597 Meade, B.J., Peden-Adams, M.M., Luebke, R.W., Luster, M.I., 2009. Immunotoxicity of
598 perfluorooctanoic acid and perfluorooctane sulfonate and the role of peroxisome
599 proliferator-activated receptor alpha. Crit. Rev. Toxicol. 39, 76-94.

600 https://doi.org/10.1080/10408440802209804

601 Dhanasekaran, D.N., Reddy, E.P., 2008. JNK Signaling in Apoptosis. Oncogene 27, 6245-6251.
602 https://doi.org/10.1038/onc.2008.301

603 Dong, T.T.T., Lee, B.-K., 2009. Characteristics, toxicity, and source apportionment of polycylic
604 aromatic hydrocarbons (PAHs) in road dust of Ulsan, Korea. Chemosphere 74, 1245-1253.
605 https://doi.org/10.1016/j.chemosphere.2008.11.035

606 EC, 2000. Directive of the European Parliament and of the Council 2000/60/EC establishing a
607 framework for community action in the field of Water Policy. PE-CONS 3639/1/00.

608 Edgar, R., Domrachev, M., Lash, A.E., 2002. Gene Expression Omnibus: NCBI gene expression and
609 hybridization array data repository. Nucleic Acids Res. 30, 207-210.
610 https://doi.org/10.1093/nar/30.1.207

611 Eriksson, E., Baun, A., Scholes, L., Ledin, A., Ahlman, S., Revitt, M., Noutsopoulos, C., Mikkelsen, P.S.,
612 2007. Selected stormwater priority pollutants — a European perspective. Sci. Total Environ.
613 383, 41-51. https://doi.org/10.1016/j.scitotenv.2007.05.028

614 Fauvelle, V., Nhu-Trang, T.-T., Feret, T., Madarassou, K., Randon, J., Mazzella, N., 2015. Evaluation of
615 Titanium Dioxide as a Binding Phase for the Passive Sampling of Glyphosate and
616 Aminomethyl Phosphonic Acid in an Aquatic Environment. Anal. Chem. 87, 6004-6009.
617 https://doi.org/10.1021/acs.analchem.5b00194

618 Forbes, V.E., Palmqvist, A., Bach, L., 2006. The use and misuse of biomarkers in ecotoxicology.
619 Environ. Toxicol. Chem. 25, 272-280. https://doi.org/10.1897/05-257R.1

620 Fuzzen, M.L.M., Bennett, C.J., Tetreault, G.R., McMaster, M.E., Servos, M.R., 2015. Severe intersex is
621 predictive of poor fertilization success in populations of rainbow darter (Etheostoma
622 caeruleum). Aquat. Toxicol. Amst. Neth. 160, 106-116.
623 https://doi.org/10.1016/j.aquatox.2015.01.009

624 Galloway, T.S., Depledge, M.H., 2001. Immunotoxicity in Invertebrates: Measurement and
625 Ecotoxicological Relevance. Ecotoxicology 10, 5-23.
626 https://doi.org/10.1023/A:1008939520263

627 Gillis, P.L., Higgins, S.K., Jorge, M.B., 2014. Evidence of oxidative stress in wild freshwater mussels
628 (Lasmigona costata) exposed to urban-derived contaminants. Ecotoxicol. Environ. Saf. 102,
629 62-69. https://doi.org/10.1016/j.ecoenv.2013.12.026
630 Göbel, P., Dierkes, C., Coldewey, W.G., 2007. Storm water runoff concentration matrix for urban areas. J. Contam. Hydrol., Issues in urban hydrology: The emerging field of urban contaminant hydrology 91, 26-42. https://doi.org/10.1016/j.jconhyd.2006.08.008

631
632
633 Gonzalez, P., Pierron, F., 2015. Omics in Aquatic Ecotoxicology: The Ultimate Response to Biological Questions?, in: Amiard-Triquet, C., Amiard, J.-C., Mouneyrac, C. (Eds.), Aquatic Ecotoxicology. Academic Press, pp. 183-203. https://doi.org/10.1016/B978-0-12-800949-9.00008-5

636 Grabherr, M.G., Haas, B.J., Yassour, M., Levin, J.Z., Thompson, D.A., Amit, I., Adiconis, X., Fan, L., Raychowdhury, R., Zeng, Q., Chen, Z., Mauceli, E., Hacohen, N., Gnirke, A., Rhind, N., di Palma, F., Birren, B.W., Nusbaum, C., Lindblad-Toh, K., Friedman, N., Regev, A., 2011. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29, 644-652. https://doi.org/10.1038/nbt.1883

641 Gyrd-Hansen, M., Meier, P., 2010. IAPs: from caspase inhibitors to modulators of NF-kappaB, inflammation and cancer. Nat. Rev. Cancer 10, 561-574. https://doi.org/10.1038/nrc2889

643 Hartl, F.U., Hayer-Hartl, M., 2002. Molecular chaperones in the cytosol: from nascent chain to folded protein. Science 295, 1852-1858. https://doi.org/10.1126/science.1068408

645 Haseba, T., Duester, G., Shimizu, A., Yamamoto, I., Kameyama, K., Ohno, Y., 2006. In vivo contribution of Class III alcohol dehydrogenase (ADH3) to alcohol metabolism through activation by cytoplasmic solution hydrophobicity. Biochim. Biophys. Acta BBA - Mol. Basis Dis. 1762, 276-283. https://doi.org/10.1016/j.bbadis.2005.11.008

649 Hay, N., 2011. Interplay between FOXO, TOR, and Akt. Biochim. Biophys. Acta BBA - Mol. Cell Res., PI3K-AKT-FoxO axis in Cancer and Aging 1813, 1965-1970. https://doi.org/10.1016/j.bbamcr.2011.03.013

652 Heberer, T., 2002. Occurrence, fate, and removal of pharmaceutical residues in the aquatic environment: a review of recent research data. Toxicol. Lett. 131, 5-17. https://doi.org/10.1016/S0378-4274(02)00041-3

655 Huber, M., Welker, A., Helmreich, B., 2016. Critical review of heavy metal pollution of traffic area runoff: Occurrence, influencing factors, and partitioning. Sci. Total Environ. 541, 895-919. https://doi.org/10.1016/j.scitotenv.2015.09.033

658 Jasinska, E.J., Goss, G.G., Gillis, P.L., Van Der Kraak, G.J., Matsumoto, J., de Souza Machado, A.A., Giacomin, M., Moon, T.W., Massarsky, A., Gagné, F., Servos, M.R., Wilson, J., Sultana, T., Metcalfe, C.D., 2015. Assessment of biomarkers for contaminants of emerging concern on aquatic organisms downstream of a municipal wastewater discharge. Sci. Total Environ. 530, 140-153. https://doi.org/10.1016/j.scitotenv.2015.05.080

663 Kanehisa, M., Sato, Y., Morishima, K., 2016. BlastKOALA and GhostKOALA: KEGG Tools for Functional Characterization of Genome and Metagenome Sequences. J. Mol. Biol. 428, 726-731. https://doi.org/10.1016/j.jmb.2015.11.006

666 Kayhanian, M., Stransky, C., Bay, S., Lau, S.-L., Stenstrom, M.K., 2008. Toxicity of urban highway runoff with respect to storm duration. Sci. Total Environ. 389, 386-406. https://doi.org/10.1016/j.scitotenv.2007.08.052

669 Kayhanian Masoud, Singh Amardeep, Suverkropp Claus, Borroum Steve, 2003. Impact of Annual Average Daily Traffic on Highway Runoff Pollutant Concentrations. J. Environ. Eng. 129, 975-990. https://doi.org/10.1061/(ASCE)0733-9372(2003)129:11(975)

672 Ke, Q., Costa, M., 2006. Hypoxia-Inducible Factor-1 (HIF-1). Mol. Pharmacol. 70, 1469-1480. https://doi.org/10.1124/mol.106.027029

674 Keay, J., Thornton, J.W., 2009. Hormone-Activated Estrogen Receptors in Annelid Invertebrates: Implications for Evolution and Endocrine Disruption. Endocrinology 150, 1731-1738. https://doi.org/10.1210/en.2008-1338

677 Krein, A., Schorer, M., 2000. Road runoff pollution by polycyclic aromatic hydrocarbons and its contribution to river sediments. Water Res. 34, 4110-4115. https://doi.org/10.1016/S0043-1354(00)00156-1
680 Lagesen, K., Hallin, P., Rødland, E.A., Stærfeldt, H.-H., Rognes, T., Ussery, D.W., 2007. RNAmmer: consistent and rapid annotation of ribosomal RNA genes. Nucleic Acids Res. 35, 3100–3108.
681 https://doi.org/10.1093/nar/gkm160

682 Laplante, M., Sabatini, D.M., 2009. mTOR signaling at a glance. J. Cell Sci. 122, 3589–3594.
684 https://doi.org/10.1242/jcs.051011

685 Lassudrie, M., Soudant, P., Nicolas, J.-L., Miner, P., Le Grand, J., Lambert, C., Le Goïc, N., Hégaret, H., Fabioux, C., 2016. Exposure to the toxic dinoflagellate Alexandrium catenella modulates juvenile oyster Crassostrea gigas hemocyte variables subjected to different biotic conditions. Fish Shellfish Immunol. 51, 104–115. https://doi.org/10.1016/j.fsi.2016.02.017

689 Le Coadou, L., Le Ménach, K., Labadie, P., Dévier, M.-H., Pardon, P., Augagneur, S., Budzinski, H., 2017. Quality survey of natural mineral water and spring water sold in France: Monitoring of hormones, pharmaceuticals, pesticides, perfluoroalkyl substances, phthalates, and alkylphenols at the ultra-trace level. Sci. Total Environ. 603–604, 651–662.
693 https://doi.org/10.1016/j.scitotenv.2016.11.174

694 Li, H., Durbin, R., 2009. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25, 1754–1760. https://doi.org/10.1093/bioinformatics/btp324

696 Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., Marth, G., Abecasis, G., Durbin, R., 2009. The Sequence Alignment/Map format and SAMtools. Bioinformatics 25, 2078–2079.
698 https://doi.org/10.1093/bioinformatics/btp352

699 Love, M.I., Huber, W., Anders, S., 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550. https://doi.org/10.1186/s13059-014-0550-8

701 Lydeard, C., Cowie, R.H., Ponder, W.F., Bogan, A.E., Bouchet, P., Clark, S.A., Cummings, K.S., Frest, T.J., Gargominy, O., Herbert, D.G., Hershler, R., Perez, K.E., Roth, B., Seddon, M., Strong, E.E., Thompson, F.G., 2004. The Global Decline of Nonmarine Mollusks. BioScience 54, 321–330.
704 https://doi.org/10.1641/0006-3568(2004)054[0321:TGDONM]2.0.CO;2

705 Mariette, J., Noirot, C., Nabihoudine, I., Bardou, P., Hoede, C., Djari, A., Cabau, C., Klopp, C., 2014. RNAbrowse: RNA-Seq De Novo Assembly Results Browser. PLOS ONE 9, e96821.
707 https://doi.org/10.1371/journal.pone.0096821

708 Matozzo, V., Fabrello, J., Masiero, L., Ferraccioli, F., Finos, L., Pastore, P., Di Gangi, I.M., Bogialli, S., 2018. Ecotoxicological risk assessment for the herbicide glyphosate to non-target aquatic species: A case study with the mussel Mytilus galloprovincialis. Environ. Pollut. 233, 623–632.
711 https://doi.org/10.1016/j.envpol.2017.10.100

712 McLellan, S.L., Fisher, J.C., Newton, R.J., 2015. The microbiome of urban waters. Int. Microbiol. Off. J. Span. Soc. Microbiol. 18, 141–149. https://doi.org/10.2436/20.1501.01.244

714 Mendoza, M.C., Er, E.E., Blenis, J., 2011. The Ras-ERK and PI3K-mTOR Pathways: Cross-talk and Compensation. Trends Biochem. Sci. 36, 320–328. https://doi.org/10.1016/j.tibs.2011.03.006

716 Mezghani-Chaari, S., Machreki-Ajmi, M., Tremolet, G., Kellner, K., Geffard, A., Minier, C., Hamza-Chaffai, A., 2015. The endocrine-disrupting effect and other physiological responses of municipal effluent on the clam Ruditapes decussatus. Environ. Sci. Pollut. Res. Int. 22, 19716–19728. https://doi.org/10.1007/s11356-015-5199-7

720 Miège, C., Choubert, J.M., Ribeiro, L., Eusèbe, M., Coquery, M., 2009. Fate of pharmaceuticals and personal care products in wastewater treatment plants – Conception of a database and first results. Environ. Pollut., Special Issue Section: Ozone and Mediterranean Ecology: Plants, People, Problems 157, 1721–1726. https://doi.org/10.1016/j.envpol.2008.11.045

724 Milan, M., Dalla Rovere, G., Smits, M., Ferraresso, S., Pastore, P., Marin, M.G., Bogialli, S., Patarnello, T., Bargelloni, L., Matozzo, V., 2018. Ecotoxicological effects of the herbicide glyphosate in non-target aquatic species: Transcriptional responses in the mussel Mytilus galloprovincialis. Environ. Pollut. 237, 442–451. https://doi.org/10.1016/j.envpol.2018.02.049

728 Morris, L., Colombo, V., Hassell, K., Kellar, C., Leahy, P., Long, S.M., Myers, J.H., Pettigrove, V., 2017. Municipal wastewater effluent licensing: A global perspective and recommendations for best practice. Sci. Total Environ. 580, 1327–1339. https://doi.org/10.1016/j.scitotenv.2016.12.096
731 Morrison, D.K., 2012. MAP Kinase Pathways. Cold Spring Harb. Perspect. Biol. 4, a011254.

732 https://doi.org/10.1101/cshperspect.a011254

733 Mouthon, J., 2001. Life cycle and population dynamics of the Asian clam Corbicula fluminea (Bivalvia:

734 Corbiculidae) in the Saone River at Lyon (France). Hydrobiologia 452, 109–119.

735 https://doi.org/10.1023/A:1011980011889

736 Muñoz, I., José Gómez, M., Molina-Díaz, A., Huijbregts, M.A.J., Fernández-Alba, A.R., García-Calvo, E.,

737 2008. Ranking potential impacts of priority and emerging pollutants in urban wastewater

738 through life cycle impact assessment. Chemosphere 74, 37–44.

739 https://doi.org/10.1016/j.chemosphere.2008.09.029

740 Nardinocchi, L., Puca, R., Guidolin, D., Belloni, A.S., Bossi, G., Michiels, C., Sacchi, A., Onisto, M.,

741 D'Orazi, G., 2009. Transcriptional regulation of hypoxia-inducible factor 1α by HIPK2 suggests

742 a novel mechanism to restrain tumor growth. Biochim. Biophys. Acta BBA - Mol. Cell Res.

743 1793, 368–377. https://doi.org/10.1016/j.bbamcr.2008.10.013

744 Pavlikova, N., Kortner, T.M., Arukwe, A., 2010. Peroxisome proliferator-activated receptors,

745 estrogenic responses and biotransformation system in the liver of salmon exposed to

746 tributyltin and second messenger activator. Aquat. Toxicol. Amst. Neth. 99, 176–185.

747 https://doi.org/10.1016/j.aquatox.2010.04.014

748 Petrie, B., Barden, R., Kasprzyk-Hordern, B., 2015. A review on emerging contaminants in

749 wastewaters and the environment: Current knowledge, understudied areas and

750 recommendations for future monitoring. Water Res., Occurrence, fate, removal and

751 assessment of emerging contaminants in water in the water cycle (from wastewater to

752 drinking water) 72, 3–27. https://doi.org/10.1016/j.watres.2014.08.053

753 Pierron, F., Normandeau, E., Defo, M.A., Campbell, P.G.C., Bernatchez, L., Couture, P., 2011. Effects

754 of chronic metal exposure on wild fish populations revealed by high-throughput cDNA

755 sequencing. Ecotoxicology 20, 1388–1399. https://doi.org/10.1007/s10646-011-0696-z

756 Poynton, H.C., Zuzow, R., Loguinov, A.V., Perkins, E.J., Vulpe, C.D., 2008. Gene Expression Profiling in

757 Daphnia magna, Part II: Validation of a Copper Specific Gene Expression Signature with

758 Effluent from Two Copper Mines in California. Environ. Sci. Technol. 42, 6257–6263.

759 https://doi.org/10.1021/es800262k

760 Quayle, D.B., Newkirk, G.F., 1989. Farming bivalve molluscs: methods for study and development,

761 Advances in world aquaculture. World Aquaculture Society in association with the

762 International Development Research Centre, Baton Rouge, LA.

763 Quevillon, E., Silventoinen, V., Pillai, S., Harte, N., Mulder, N., Apweiler, R., Lopez, R., 2005.

764 InterProScan: protein domains identifier. Nucleic Acids Res. 33, W116–W120.

765 https://doi.org/10.1093/nar/gki442

766 Regier, N., Baerlocher, L., Münsterkötter, M., Farinelli, L., Cosio, C., 2013. Analysis of the Elodea

767 nuttallii Transcriptome in Response to Mercury and Cadmium Pollution: Development of

768 Sensitive Tools for Rapid Ecotoxicological Testing. Environ. Sci. Technol. 47, 8825–8834.

769 https://doi.org/10.1021/es401082h

770 Regoli, F., Giuliani, M.E., 2014. Oxidative pathways of chemical toxicity and oxidative stress

771 biomarkers in marine organisms. Mar. Environ. Res., Sensing the marine environment using

772 different animal models and levels of complexitySensing the marine environment using

773 different animal models and levels of complexity 93, 106–117.

774 https://doi.org/10.1016/j.marenvres.2013.07.006

775 Revitt, D.M., Lundy, L., Coulon, F., Fairley, M., 2014. The sources, impact and management of car

776 park runoff pollution: A review. J. Environ. Manage. 146, 552–567.

777 https://doi.org/10.1016/j.jenvman.2014.05.041

778 Riascos, J.M., Avalos, C.M., Pacheco, A.S., Heilmayer, O., 2012. Testing stress responses of the bivalve

779 Protothaca thaca to El Niño–La Niña thermal conditions. Mar. Biol. Res. 8, 654–661.

780 https://doi.org/10.1080/17451000.2011.653367

781 Rozen, S., Skaletsky, H., 2000. Primer3 on the WWW for general users and for biologist programmers.

782 Methods Mol. Biol. 132, 365–386.
783 Salazar, M.H., Salazar, S.M., 2005. Field experiments with caged bivalves to assess chronic exposure and toxicity, in: Techniques in Aquatic Toxicology. Ostrander, G.K., pp. 117–136.

785 Santos, E.M., Ball, J.S., Williams, T.D., Wu, H., Ortega, F., van Aerle, R., Katsiadaki, I., Falciani, F., Viant, M.R., Chipman, J.K., Tyler, C.R., 2010. Identifying health impacts of exposure to copper using transcriptomics and metabolomics in a fish model. Environ. Sci. Technol. 44, 820–826.
https://doi.org/10.1021/es902558k

789 Sato, T., Qadir, M., Yamamoto, S., Endo, T., Zahoor, A., 2013. Global, regional, and country level need for data on wastewater generation, treatment, and use. Agric. Water Manag. 130, 1–13.
https://doi.org/10.1016/j.agwat.2013.08.007

792 Schmidt, W., Rainville, L.-C., McEneff, G., Sheehan, D., Quinn, B., 2014. A proteomic evaluation of the effects of the pharmaceuticals diclofenac and gemfibrozil on marine mussels (Mytilus spp.): evidence for chronic sublethal effects on stress-response proteins. Drug Test. Anal. 6, 210–219.
https://doi.org/10.1002/dta.1463

796 Segner, H., Wenger, M., Möller, A.M., Köllner, B., Casanova-Nakayama, A., 2011. Immunotoxic effects of environmental toxicants in fish - how to assess them? Environ. Sci. Pollut. Res. Int. 19, 2465–2476.
https://doi.org/10.1007/s11356-012-0978-x

799 Sidhu, J.P.S., Hodgers, L., Ahmed, W., Chong, M.N., Toze, S., 2012. Prevalence of human pathogens and indicators in stormwater runoff in Brisbane, Australia. Water Res., Special Issue on Stormwater in urban areas 46, 6652–6660.
https://doi.org/10.1016/j.watres.2012.03.012

802 Simão, F.A., Waterhouse, R.M., Ioannidis, P., Kriventseva, E.V., Zdobnov, E.M., 2015. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs.
Bioinformatics 31, 3210–3212.
https://doi.org/10.1093/bioinformatics/btv351

805 Sokolova, I.M., Frederich, M., Bagwe, R., Lannig, G., Sukhotin, A.A., 2012. Energy homeostasis as an integrative tool for assessing limits of environmental stress tolerance in aquatic invertebrates. Mar. Environ. Res. 79, 1–15.
https://doi.org/10.1016/j.marenvres.2012.04.003

808 Song, L., Wang, L., Qiu, L., Zhang, H., 2010. Bivalve Immunity, in: Invertebrate Immunity, Advances in Experimental Medicine and Biology. Springer, Boston, MA, pp. 44–65.
https://doi.org/10.1007/978-1-4419-8059-5_3

811 Sousa, R., Nogueira, A.J.A., Gaspar, M.B., Antunes, C., Guilhermino, L., 2008. Growth and extremely high production of the non-indigenous invasive species Corbicula fluminea (Müller, 1774): Possible implications for ecosystem functioning. Estuar. Coast. Shelf Sci. 80, 289–295.
https://doi.org/10.1016/j.ecss.2008.08.006

815 Supek, F., Bošnjak, M., Škunca, N., Šmuc, T., 2011. REVIGO Summarizes and Visualizes Long Lists of Gene Ontology Terms. PLOS ONE 6, e21800.
https://doi.org/10.1371/journal.pone.0021800

817 Ternes, T.A., Stumpf, M., Mueller, J., Haberer, K., Wilken, R.-D., Servos, M., 1999. Behavior and occurrence of estrogens in municipal sewage treatment plants — I. Investigations in Germany, Canada and Brazil. Sci. Total Environ. 225, 81–90.
https://doi.org/10.1016/S0048-9697(98)00334-9

821 Tran, D., Boudou, A., Massabuau, J.-C., 2000. Mechanism for maintaining oxygen consumption under varying oxygenation levels in the freshwater clam *Corbicula fluminea*. Can. J. Zool. 78, 2027–2036.
https://doi.org/10.1139/z00-139

824 Uren Webster, T.M., Laing, L.V., Florance, H., Santos, E.M., 2014. Effects of Glyphosate and its Formulation, Roundup, on Reproduction in Zebrafish (Danio rerio). Environ. Sci. Technol. 48, 1271–1279.
https://doi.org/10.1021/es404258h

827 Vidal, M.-L., Bassères, A., Narbonne, J.-F., 2002. Seasonal variations of pollution biomarkers in two populations of Corbicula fluminea (Müller). Comp. Biochem. Physiol. Part C Toxicol. Pharmacol. 131, 133–151.
https://doi.org/10.1016/S1532-0456(01)00291-5

830 Zgheib, S., Moilleron, R., Chebbo, G., 2011. Influence of the land use pattern on the concentrations and fluxes of priority pollutants in urban stormwater. Water Sci. Technol. J. Int. Assoc. Water Pollut. Res. 64, 1450–1458.
https://doi.org/10.2166/wst.2011.770

833 Zhang, Y.L., Dong, C., 2005. MAP kinases in immune responses. Cell. Mol. Immunol. 2, 20–27.
834 Zhao, Y., Hu, X., Liu, Y., Dong, S., Wen, Z., He, W., Zhang, S., Huang, Q., Shi, M., 2017. ROS signaling

835 under metabolic stress: cross-talk between AMPK and AKT pathway. Mol. Cancer 16.

836 https://doi.org/10.1186/s12943-017-0648-1

837

AUTHORS VERSION