Proteomic analysis and biochemical alterations in marine mussel gills after exposure to the organophosphate flame retardant TDCPP

1 Proteomic analysis and biochemical alterations in marine mussel gills after exposure to the organophosphate flame retardant TDCPP Paula Sánchez-Marín1,2,*, Leticia Vidal-Liñán1,2,a, Laura Emilia Fernández- González1,3,4,a, Rosa Montes5, Rosario Rodil5, José Benito Quintana5, Mónica Carrera6, Jesús Mateos6, Angel P. Diz3,4,b and Ricardo Beiras1,4,b 1Department of Ecology and Animal Biology, University of Vigo, 36310 Vigo, Galicia, Spain. 2Centro Oceanográfico de Vigo, Instituto Español de Oceanografía, 36390 Vigo, Spain. 3Department of Biochemistry, Genetics and Immunology, University of Vigo, 36310 Vigo, Spain. 4Marine Research Centre, University of Vigo (CIM-UVIGO), Isla de Toralla, Vigo, Spain. 5Department of Analytical Chemistry, Nutrition and Food Sciences, Institute of Research in Chemical and Biological Analysis (IAQBUS), Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Spain. 6Institute of Marine Research (IIM), Spanish National Research Council (CSIC) , 36208 Vigo, Spain. a Authors with equal contribution. b Authors with equal contribution. * Corresponding author e-mail: paula.[email protected] This is a post-peer-review, pre-copyedit version of an article published in [Aquatic Toxicology]. The final authenticated version is available online at: http://dx.doi.org/[10.1016/j.aquatox.2020.105688] 2 Abstract Organophosphate flame retardants (OPFRs) are (re-)emergent environmental pollutants increasingly being used because of the restriction of other flame retardants. The chlorinated OPFR, tris(1,3-dichloro-2-propyl) phosphate (TDCPP) is among those of highest environmental concern, but its potential effects in the marine environment have rarely been investigated. We exposed a widely used sentinel marine mussel species, Mytilus galloprovincialis, to 10 µg L-1 of TDCPP during 28 days and studied: (i) the kinetics of bioaccumulation and elimination of the compound, (ii) the effect on two molecular biomarkers, glutathione S-transferase (GST) and acetylcholinesterase (AChE) activities, and (iii) proteomic alterations in the gills, following an isobaric labeling quantitative shotgun proteomic approach, at two exposure times (7 and 28 days). Uptake and elimination of TDCPP by mussels were very fast, and the bioconcentration factor of this compound in mussels was 147 L kgww-1, confirming that this compound is not very bioaccumulative, as predicted by its chemical properties. GST activity was not affected by TDCPP exposure, but AChE activity was inhibited by TDCPP at both 7 and 28 days of exposure. Proteomic analysis revealed subtle effects of TDCPP in mussel gills, since few proteins (less than 2% of the analysed proteome) were significantly affected by TDCPP, and effect sizes were low. The most relevant effects detected were the up-regulation of epimerase family protein SDR39U1, an enzyme that could be involved in detoxification processes, at both exposure times, and the down-regulation of receptor-type tyrosine-protein phosphatase N2-like (PTPRN2) after 7 days of exposure, which is involved in neurotransmitter secretion and might be related to the neurotoxicity described for this compound. Exposure time rather than TDCPP exposure was the most important driver of protein abundance changes, with 33% of the proteome being affected by this factor, suggesting that stress caused by laboratory conditions could be 3 an important confounding factor that needs to be controlled in similar ecotoxicology studies. Proteomic data are available via ProteomeXchange with identifier PXD019720. Keywords: Proteomics, isobaric labelling, tandem mass tag (TMT), biomarkers, flame retardants, marine mussels 4 1. Introduction Flame retardants are widely used in several consumer and industrial products like textile, polyvinyl chloride plastics, cellulose, coatings, polyester resins, polyurethane foams, etc. These compounds are designed to prevent combustion and to delay the spread of fire after ignition, and can reach the environment by volatilization, leaching and abrasion during the entire lifespan of the products to which they are added (Betts, 2008; Van der Veen and de Boer, 2012). Several compounds that have been used as flame retardants since the 60s, such as PCBs or PBDEs, have been banned or restricted afterwards because of their persistence in the environment, and their bioaccumulative and toxic effects (Beiras, 2018). Alternative compounds are being used in substitution, such as organophosphate flame retardants (OPFRs), which are re-emerging as environmental pollutants in recent years (Wei et al., 2015; Pantelaki and Voutsa, 2019). Within existing OPFRs, chlorinated ones, such as tris(2-chloroethyl)phosphate (TCEP), tris(chloropropyl)phosphate (TCPP) and tris(1,3-dichloro-2-propyl)phosphate (TDCPP), are within the most frequently used, but are also the ones of highest environmental concern, due to their persistence and reported toxicological effects, such as neurotoxicity and carcinogenicity (Reemtsma et al., 2008; Van der Veen and de Boer, 2012). Within them, TDCPP was considered of highest priority for aquatic life risk assessment by the European Union (EU, 2008), it is persistent and has been detected in rivers, coastal waters and marine biota (Aznar-Alemany et al., 2018; Pantelaki and Voutsa, 2019). Toxicological studies with mammals and cell lines have pointed towards its carcinogenicity and neurotoxicity (WHO, 1998; ATSDR, 2009; Dishaw et al., 2011). Despite its potential effects in the aquatic environment are less known, studies with 5 zebrafish have shown developmental and endocrine alterations upon exposure to this compound (McGee et al., 2012; Liu et al., 2013; Wang et al., 2015). Blue mussels (Mytilus sp.) are frequently used in marine pollution monitoring (Sericano, 2000; Beyer et al., 2017). Because they are sessile, widespread, easy to sample, and filter-feeders that bioaccumulate pollutants in their tissues, they are a preferred biomonitoring species not only for the measurement of pollutants in their tissues, but also for the evaluation of biomarkers of exposure or effect, as early warning tools for the detection of deleterious effects of pollutants in the marine environment (Beliaeff and Burgeot, 2002; Vidal-Liñán et al., 2010; Davies and Vethaak, 2012). In the last years, the development of proteomic tools has allowed the study of regulation level of hundreds or even thousands of proteins in mussel tissues upon exposure to contaminants, which has opened the door to the discovery of new biomarkers and to increase knowledge about the molecular mechanisms of toxicity of contaminants (Campos et al., 2012; Gouveia et al., 2019). Proteomic studies of the effects of pollutants are usually performed at a fixed exposure time (e.g. (Apraiz et al., 2006; Duroudier et al., 2019). However, the exposure time needed for a contaminant to elicit a response is not always known a priori, and even the effects of the contaminant can be different depending on the exposure time (Oliveira et al., 2016). Laboratory experiments with mussels have revealed that protein abundance was more affected by exposure time (including control organisms) than by exposure to a pollutant (Oliveira et al., 2016). However, the specific proteomic changes caused by exposure time, and their possible confounding effects on the outcomes of this type of experiments have not yet been evaluated. As mussels pump water for feeding and breathing purposes, their gills represent the primary organ in contact with pollutants, and are a sensitive tissue for detecting 6 toxicological effects, as previously demonstrated with the use of molecular biomarkers (Vidal-Liñán et al., 2010; Vidal-Liñán and Bellas, 2013; Vidal-Liñán et al., 2015). The objective of this work was to study the bioaccumulation of TDCPP and the effects in the gill proteome produced by this re-emerging pollutant in the marine mussel Mytilus galloprovincialis, at two different exposure times, 7 and 28 days, to test for short vs. long term exposure, and also to test for the effects of exposure time itself on mussel proteome. Likewise, this study focused on traditional biomarkers assessed in gills of exposed mussels to detect the early onset of adverse effects, including a neurotoxic response as acetylcholinesterase (AChE) and a Phase II detoxification enzyme: glutathione S-transferase (GST). 2. Material and methods 2.1. Mussel exposure Mussels (Mytilus galloprovincialis) were obtained from aquaculture from a clean area in the Ría de Vigo (Galicia, Spain) in November 2015. Obtained specimens corresponded to a common stock of reared mussels, of a similar age, and were additionally selected using the narrowest size range that was possible from the available mussels, resulting in mussels from 38 to 47 mm length used for the experiments. This size range is not expected to introduce much heterogeneity in mussels physiology since the age of the mussels was similar. Mussels were acclimated to laboratory conditions in an open system with running seawater at 13 ºC during one week. During the acclimation period mussels were fed twice per week with a phytoplankton mixture of three microalgae species: Chaetoceros neogracilis, Rhodomonas lens and Phaeodactylum 7 tricornutum, in a proportion 50% flagellates 50% diatoms, representing 1% of mussel tissue dry weight. Two experiments were performed, a preliminary experiment, performed to check for the stability of different OPFRs in seawater under laboratory conditions, and a long term exposure experiment, testing for the effects of TDCPP on mussels. In the preliminary experiment, six 30-L glass aquaria were filled with 1 µm-filtered and UV treated seawater and spiked with 160 µL of a stock solution of TDCPP, TCEP and TCPP (Sigma-Aldrich) dissolved in acetone (HPLC grade, Sigma-Aldrich), to give a concentration of 50 nM (corresponding to 21.5, 14 and 16 µg L-1 of each compound, respectively). After 24 h, the water was renewed and water samples were taken in 100 mL glass bottles (t0). Then, 15 mussels were introduced in three of the aquaria and water samples were taken after 24, 48 and 72 h to test the stability of the compounds in water under the experimental conditions, both in the presence and absence of mussels. For the exposure experiment, mussels were exposed to TDCPP at a concentration of 10 µg L-1, to a solvent control (SC) (acetone) or to control (unspiked) seawater (1 µmfiltered and UV treated seawater) during 28 days. The concentration of TDCPP was chosen on the basis of described chronic effects for other aquatic species (Liu et al., 2013; Wang et al., 2015; Wang et al., 2015) (see discussion for details). Mussels exposed to TDCPP were additionally subjected to a depuration period of 14 days in clean seawater. TDCPP was added from a stock solution prepared in acetone, the volume added being 100 µL per 30 L aquaria. The same volume of acetone was added to SC aquaria. Mussel density was 20 mussels per each 30 L aquaria at the beginning of the experiment and varied between 15 ‒ 20 mussels/30 L aquaria during the 28 days 8 exposure period. Mussels were exposed at 15 ºC and using a 10:14 h light:dark photoperiod, and filtered air was gently bubbled to assure good oxygen levels. Water was renewed twice per week after feeding the mussels, using the following procedure: First, T, S, pH and O2 were measured before every water renewal (and also at the beginning of the experiment), resulting in stable values: T = 15.6 ± 0.7 ºC, S = 35.5 ± 0.5 ppt, pH = 7.9 ± 0.1, O2 = 7.9 ± 0.6 mg L-1. Then mussels were fed during 1 h before every water renewal with a mixture of phytoplankton representing 1% of mussel tissue dry weight. After feeding, old water and faeces were removed by siphoning, and mussels were then carefully taken off the aquaria before refilling with fresh seawater that was spiked with SC and TDCPP additions. After letting 5 min for homogenization of the spiked solution, mussels were introduced back in the aquaria. Water samples were taken before and after every water renewal to check for actual TDCPP concentrations during exposure. Eight mussels exposed to TDCPP were collected at time 0, 3, 7, 14, 21 and 28 days during exposure and at time 3, 7 and 14 during depuration for chemical analysis of TDCPP in their tissues. For proteomic and biochemical analysis, 30 mussels were collected at the beginning of the experiment (T0), and after 7 (T7) and 28 (T28) days of exposure from TDCPP and SC treatments. In addition, 19 individuals were sampled after 28 days of exposure to control seawater. Mussel sampling was homogeneously distributed among all aquaria available at each sampling time for each treatment, to avoid for "aquarium effects" that could bias the results. Due to space limitations, a single aquarium was used for the control treatment, so results from this treatment were included for comparison but were not given the same grade of importance in data reporting as TDCPP and SC ones, because of potential flaws due to underreplication (Krull et al., 2013; Tincani et al., 2017). The gonads and gills were dissected, and gills were divided in two sub-samples, quick-frozen in liquid 9 nitrogen and kept at -80 ºC until further processing for proteomic and biochemical analysis. Gonads were fixed in Davidson formaldehyde solution for 24 h and processed for histology analysis. Water contaminated with OPFRs produced during the experiment was treated by recirculation through an activated carbon filter during at least 48 h. Water samples of treated wastewater were taken to check for the effectiveness of the decontamination process. 2.2. Chemical analyses The analysis of OPFRs in the different matrices was performed by liquid chromatography-tandem mass spectrometry (LC-MS/MS) in a Varian (ProStar 210 and 320-MS) triple quadrupole instrument. Deuterium labelled analogs, i.e. TCEP-d12, TCCP-d18 and TDCPP-d15 (Wellington Laboratories) were used as surrogate internal standards (IS). The separation was carried out in a Synergy fusion RP 4 μm C18 column (100×2 mm) connected to a C18 (2×4 mm) guard cartridge (Phenomenex). The mobile phases consisted of Milli-Q water (0.1% formic acid) (A) and methanol (0.1% formic acid) (B) and the flow rate and column temperature were 0.2 mL min-1 and 35 °C, respectively. The gradient elution started with 5% B (0-1min), increasing to 100% B (1- 6 min), held for 2 min. Subsequently, it returned to initial conditions (5% B), held for 5 min for column back-conditioning. Nitrogen was used as nebulizing (55 psi) and drying gas (200 °C, 18 psi) in the electrospray (ESI) source and Argon (99.999%) was employed as collision gas (2 mTorr). The voltage of the ESI needle was 5,000 V and the temperature of the ESI housing was set at 55 °C. Analysis was performed in the ESI positive and multiple reaction monitoring (MRM) modes, using two (precursor > fragment) transitions per analyte and one transition per IS. The most intense transition 16 experiment auto-normalization). Proteins that were only detected in one of the three proteomic experiments were excluded (586 proteins), resulting in 4214 proteins in the dataset. A second normalization was performed in order to compare results from different proteomic experiments (inter-experiment normalization), consisting in calculating the proportion of intensities of each reporter ion in relation to the sum of all six reporter ions for each protein. This normalization strategy was inspired on the CONSTANd method described in Maes et al. (2016). Proteins with missing values for one or more reporter ions within each proteomic experiment were rejected, resulting in a final list of 4191 proteins. In the absence of inter-experiment normalization, results were biased by the proteomic experiment, as observed in the clustering of samples (figure S1a in supporting information). After performing this inter-experiment normalization, sample clustering was no longer affected by the proteomic experiment but by experimental conditions (especially exposure time), and the two technical replicates were clustered together, indicating that the normalization method was adequate (figure S1b) and that the biological variability was higher than the technical variability. Indeed, mean CV% was 6% for pooled biological replicates and 4% for technical replicates. In order to meet normality and homocedasticity of the data, normalized protein abundance values were log transformed for all subsequent statistical treatment, including tests for differential expression, principal component analysis (PCA), hierarchical clustering and heatmap analysis. Differential expression of each protein was tested by ANOVA followed by post-hoc test Tukey HSD using the aov (analysis of variance) and TukeyHSD functions in R. Different corrections to account for the multiple hypothesis testing problem were calculated by using the SGoF+ software v.3.8 (Carvajal-Rodriguez and de Uña-Alvarez, 2011) following the procedure and rationale discussed in Diz et al. (2011). 17 Hierarchical Clustering analysis was performed using the dist and hclust functions in R, using the euclidean distance and average linkage method, and heatmaps and PCA analysis were obtained with ClustVis (Metsalu and Vilo, 2015). Detected proteins were re-annotated using the program OmicsBox (version 1.3) against NCBI non-redundant protein sequences (Update date: 7 May 2020; Number of sequences:281,252,422), using an E-value threshold of 10-5. Annotation of the 17 differentially expressed (DE) proteins selected for presentation in the results section was checked manually and updated if necessary by using the NCBI Blastp search web tool on May 2020. Enrichment analysis was performed in OmicsBox using the one-tailed Fisher exact test to look for GO terms significantly enriched in the test set (DE proteins at 28 days of exposure compared to 0 and 7 days) compared to the reference set (all other proteins). Default settings were used and the enrichment results were reduced to most specific GO terms. 3. Results 3.1. Stability of selected OPFRs in seawater Results from the preliminary experiments showed that the three compounds (TDCPP, TCEP and TCPP) were kept at constant concentrations in the water in the absence of mussels in the aquaria during the three days tested (concentrations were maintained over 93% of initial values). In aquaria that contained mussels (15 mussels / 30 L), TCEP and TCPP values after three days were above 98% of initial nominal concentrations, and only TDCPP showed a slight (16%) decrease in water concentrations, corresponding to 18 a rate of decrease of TDCPP in water of 1.15 µg L-1 per day. On this basis, TDCPP was chosen for the bioaccumulation and toxicity experiment, as it was the chemical with the highest bioaccumulation potential in mussels, while a frequency for water changes of twice per week was considered enough to maintain reasonable (less than 20% change) stability in water concentrations. This was later confirmed by checking TDCPP concentrations in the exposure experiment taking samplings of old water (before renewal) and new water (with fresh spike) during every water change event. Average concentration of TDCPP measured in the 10 µg L-1 treatment was 8.7 ± 0.5 µg L-1 (n = 17), and was kept at concentrations between 7.75 and 9.50 µg L-1 throughout the experiment. TDCPP in control, solvent control and wastewater treated with activated carbon was below the limit of detection of the technique (0.03 µg L-1). 3.2. TDCPP bioaccumulation Only three mussels died during exposure, one after 1 day in the SC treatment, another after 7 days of exposure to TDCPP and another after 22 days of exposure to control water. Mussels exposed to TDCPP accumulated the compound at a rate of 789 ± 83 L kgdw-1 d-1 and showed very fast elimination kinetics (ke = 1.07 ± 0.11 d-1), reaching the steady state already after 5 days of exposure and showing rapid elimination, with 98% of the compound being eliminated during the first 3 days of depuration (Figure 1). The bioconcentration factor (BCF) was 147 L kgww-1, calculated on the basis of nominal water content of mussels (corresponding to an 80% of humidity in mussels’ tissues) to convert dry to wet weight values. 3.3. Histological analysis of the gonads 19 Histological analysis of the gonads showed the expected sex proportion of 1:1 in all treatments (Binomial exact test, p > 0.05, see supporting information, Table S1 for detailed p-values of the test for each treatment), despite there were some individuals that could not be sexed because they presented a very initial maturation stage (stage 0 or I), without follicles or with small follicles that did not present gametogenic cells. Individuals with unknown sex represented < 10% in mussels sampled at 0 and 7 days, and from 20 to 43% in mussels sampled after 28 days of exposure (Table S1, supporting information). In general, individuals were in gametogenesis stage, with the majority of individuals (43%) being in stage I of gonadal development, or early gametogenesis, followed by st. II of gonadal development (33%), while only a 5% were mature. No significant differences were found in the gametogenesis stage distribution among treatments (Fisher exact test; p = 0.283 for females and p = 0.291 for males). Individuals used for proteomic analysis were selected on the basis of histological analysis of the gonads. Despite the proteomic analysis was made in gills, there is evidence of sex-specific responses to contaminants (Atasaral-Şahin et al., 2015; Ji et al., 2016). Therefore, in order to avoid possible bias in the results caused by dissimilar proportion of males and females in the pooled samples, proteomic analysis were performed in pools of gills of four individuals, including two females and two males in similar stage of gametogenesis. 3.4. Proteomic analysis of Mytilus galloprovincialis gills 3.4.1. Exploratory assessment of regulation of gill proteome The dataset used for differential regulation analysis comprised 4191 proteins that were identified and quantified in at least two of the three 6-plex proteomic experiments. Two way ANOVA was performed independently for each protein in order to check the effect 20 of “time” (0, 7 and 28 days) and “chemical” (Control, SC, TDCPP) factors in protein abundance levels, resulting in 1389, 355 and 306 proteins with significant (p < 0.05) differential expression due to “time”, “chemical” and “time × chemical” factors respectively (see supplementary spreadsheet "Annex_data_pvalues" with p-values for each one of the analyzed factors for each protein). After multitest correction, few proteins remained differentially expressed for the “chemical” and “interaction” factors, with only one and none of the proteins presenting q-values < 0.05 respectively. For factor “time”, on the contrary, the majority (1300) proteins presented a q-value < 0.05 and the highest q-value was only 0.055, meaning that only 5.5% of those 1389 significant DE proteins are expected to be false positives (Diz et al., 2011). Hierarchical clustering of samples (figure S1(b) in the supporting information) confirms the above ANOVA results showing a higher and broader effect at proteome level of exposure time on the basis of the observation of sample clustering according to the exposure time rather than to the chemical (Control, SC or TDCPP). That is, mussels exposed during a short term (0 and 7 days) clustered together in a different branch from those exposed during 28 days. A plot of the first two principal components of PCA analysis (jointly explaining 33.6% of data variability) also shows data grouping according to exposure time, with T28 mussels grouped apart from mussels exposed during a short term (Figure S2). Hierarchical clustering and heatmap analysis was repeated with a shorter dataset, corresponding to 1148 proteins that were differentially expressed at any treatment (after one-way ANOVA, see section 3.4.2). A clear distinction was again observed between mussels exposed during a short term (0 and 7 days) and a longer term (28 days). Additionally, samples were grouped according to the chemical (Control, SC, TDCPP) within each exposure time (Figure 2a). The two first PCs of PCA analysis explained 21 57.6% of data variability, and also showed that data were mainly grouped according to exposure time, with three distinct groups formed with mussels from T0, T7 and T28, respectively. Within T28, TDCPP-exposed samples were also separated from SC and Control samples (Figure 2b). The above results revealed that exposure time was a more important driver of changes in protein abundance than chemical exposure to TDCPP or SC, but revealed also some differences between individuals exposed to TDCPP, SC or control water, that depended on the exposure time as well, hence these proteins must be evaluated separately for each exposure time. 3.4.2. Effect of exposure time on gill proteins In order to specifically study the effects of exposure time on the gill proteome, samples were grouped according to exposure time in two groups: short exposure (including samples exposed 0 and 7 days), and long exposure (including samples exposed during 28 days), on the basis of the clear distinction of these two groups of samples observed in cluster analysis and PCA (see previous section). A total of 1450 proteins were differentially expressed (DE) at the long exposure time in comparison with the short exposure time, representing a 35% of the whole detected proteome (Figure 3). After performing multitest correction we decided to keep all 1450 DE proteins for further exploration, since the highest q-value of these proteins was rather low, q = 0.063. Enrichment analysis was performed with this dataset of 1450 DE proteins (test set) in comparison with the rest (2741 non-DE proteins, reference set). The specific biological processes, cellular components and molecular functions of proteins affected by exposure time are shown in Table 2. Within these three categories, it can be observed that all GO terms can be assigned to three main processes to which they are related: 22 protein synthesis and translation processes, energy metabolism, and motility (Table 2). In order to see more precisely if proteins related with these processes were up or downregulated in response to exposure time, the fold change (FC) was calculated. 93% of DE proteins related to motility and 70% of DE proteins related to energy metabolism were down-regulated in the long exposure time. On the contrary, the majority of DE proteins related to protein synthesis and translation processes (75%) were up-regulated in the long exposure time. 3.4.3. Effect of TDCPP on gill proteins Given the strong effect of exposure time on the gill proteome, and considering that for 306 proteins their abundance levels were significantly affected by the interaction term (time × chemical), the effect of TDCPP was evaluated separately for each exposure time. This was done by performing one-way ANOVA followed by a post-hoc test (Tukey HSD) on the six treatments (T0, T7-SC, T7-TDCPP, T28-Control, T28-SC, T28-TDCPP), and focusing the description of results on the two comparisons of interest (T7-TDCPP vs. T7-SC and T28-TDCPP vs. T28-SC), so that sample variability across all experimental treatments could be taken into account in statistical analysis (Underwood, 1997). The results of the one-way ANOVA resulted in 1148 DE proteins, and we decided to keep this number after multitest correction since the highest q-value among all DE proteins was 0.069. Differential expression of proteins in the different treatments was explored using volcano plots (Figure 4), with p-values obtained from multiple comparisons (Tukey HSD). After 7 days of exposure, mussels exposed to TDCPP presented 37 DE proteins (p < 0.05) in comparison with mussels exposed to the solvent control. After 28 days of exposure, 39 proteins were differentially expressed in the TDCPP treatment in 23 comparison with the SC (see supporting file “Annex_data_pvalues”). The number of proteins showing a significant effect after exposure to TDCPP was low, representing less than a 2% of the analyzed proteome, and effect sizes were also low (FC of DE proteins varied between 0.13 and 2.3, 90% of FC were between 0.6 and 1.9). These results suggest that any strong effect of TDCPP on gill proteome can be discarded at least under the conditions of the present experiment. A complementary analysis was performed to double check that above detected significant differences are confident findings. One-way ANOVAs were repeated after permutation of sample labels from the original dataset in order to check whether the number of DE proteins detected was higher or similar to the expected by chance (supporting information, table S2). When two treatments were compared, the number of DE proteins in the original dataset was in the same range that those observed for data with shuffled sample-id labels, with some comparisons showing significantly higher number of DE proteins in the original dataset but others not (for details, see table S2). According to this, it is possible that only few of the DE proteins observed correspond to real effects of TDCPP in mussel gills, while most of them might be due only to chance. For this reason, we decided to report and further explore only the DE proteins that showed differences at the p < 0.01 level, that have a lower probability of being false positives. Exceptionally, for proteins that were differentially regulated at both exposure times in TDCPP vs. SC comparisons, the significance level for reporting was raised to p < 0.05 (Table 3). Only 9 and 7 proteins were differentially expressed at the p < 0.01 level in TDCPP vs. SC after 7 days and 28 days of exposure, respectively (Figure 4), and only two additional proteins were differentially expressed at the p < 0.05 level in TDCPP vs. SC comparisons at both exposure times. FC and p-values of different pairwise comparisons 24 of these 17 selected proteins are shown in Table 3, and protein abundance levels are represented in Figure 5. After 7 days of exposure, three proteins were significantly up-regulated (at the p < 0.01 level) in the gills exposed to TDCPP compared to those exposed to SC: apolipophorins (FC = 1.45), mitochondrial transcription factor A (FC = 1.10) and toll-interacting protein (FC = 1.09). Six proteins were down-regulated in T7-TDCPP in comparison with T7-SC: A disintegrin and metalloproteinase with thrombospondin motifs (ADAMTS) (FC = 0.93), two dynein heavy chains (FC = 0.91 and 0.67), acetylcholinesterase (FC = 0.84), adenosine deaminase (FC = 0.72) and receptor-type tyrosine-protein phosphatase N2 (FC = 0.62). Concerning the longer exposure time (T28), three proteins were significantly (p < 0.01) and slightly up-regulated (FC < 1.3) in T28-TDCPP compared to T28-SC, and one of them (cilia- and flagellaassociated protein 47) was also up-regulated (p = 0.013) in T28-TDCPP compared to T28-Control (see Table 3). Four proteins were downregulated in the T28-TDCPP vs. T28-SC comparison: AChE, ribosomal protein S6 kinase beta-1 and two other proteins that could not be identified (Table 3). However, ribosomal protein S6 kinase beta-1 and the unknown protein of sequence ID Unigene6637_All_5 were up-regulated in T28-SC treatment in comparison with all other treatments, with effect sizes of a factor of 2 and 6, respectively (Figure 5M,N), so effects might be attributable to SC, and not to TDCPP. Only one protein, AChE, was differentially expressed at the p < 0.01 level at both sampling times, resulting in lower abundance in TDCPP exposed mussels in comparison with SC ones. However, inspection of the T28-SC vs. T28-Control comparison, reveals that at T28 this was due to up-regulation in T28-SC in comparison 25 with T28-TDCPP and T28-Control (Figure 5O). The other two proteins differentially expressed at the p < 0.05 level in TDCPP vs. SC at both sampling times were DNA topoisomerase II and an epimerase family protein. For the first one, the effect was the opposite at T7 (FC = 0.9) compared to T28 (FC = 1.11), while the epimerase was upregulated in the TDCPP treatment at both sampling times. 3.5. Molecular biomarkers Figure 6 shows AChE and GST activities in gills of mussels exposed to 10 µg L-1 of TDCPP for 7 and 28 days. For the two biomarkers analysed, baseline levels are comparable with those obtained in wild mussels from pristine areas of the NW Iberian coast (Vidal-Liñán et al., 2010; Vidal-Liñán and Bellas, 2013). According to ANOVA, both TDCPP exposure (p < 0.001) and time (p < 0.01) significantly affected AChE activity. The effect of TDCPP, when compared with SC for each exposure time, was more remarkable after 7 days (p < 0.01) than after 28 days (p < 0.05), and consisted in a decrease in AChE activity (Figure 6). Exposed animals showed no significant variation in GST activity after 7 and 28 days of TDCPP exposure compared to solvent controls (Figure 6). For this biomarker, a general decrease was observed after 28 days, disregarding treatment, but this decrease did not reach statistical significance, according to the two-way ANOVA. 4. Discussion 4.1. TDCPP bioaccumulation and elimination The BCF of TDCPP (147 L kgww-1) was within the expected range according to its water solubility and octanol-water partition coefficient (Kow) (with reported log Kow values 32 masked by these confounding factors, the effects of TDCPP would need to be still subtle (low-medium effect size) to have not been detected in this study. 4.3. Inhibition of AChE activity by TDCPP Given the similarities in chemical structure of OPFRs and organophosphate pesticides ‒known as inhibitors of AChE activity‒, the neurotoxicity of TDCPP through this mechanism was investigated and it was observed that this compound did also inhibit AChE activity in mussel gills at both exposure times (Figure 6). Inhibition of AChE activity is due to the selective binding of the xenobiotic to AChE, so that the enzyme becomes phosphorylated and inactivated for a long time, even irreversibly in the case of most organophosphorous insecticides (Moriarty, 1999; Beiras, 2018). Unselective inhibition by other type of pollutants, such as metals, surfactants, PAHs, polybrominated diphenyl ethers (BDE-47), and nonylphenol has also been observed (Guilhermino et al., 1998; Vidal-Liñán et al., 2015; Vidal-Liñán et al., 2015). Effects of contaminants on the abundance level of the enzyme have been less frequently reported than its activity, although several studies have shown covariation of both responses (Xing et al., 2010; Rhee et al., 2013; Sun et al., 2016), indicating that in many cases enzyme inhibition is caused by down-regulation of the enzyme abundance, and not by its inactivation. In the present study, we have also observed the downregulation of an AChE-like enzyme in TDCPP exposed mussels as compared to the SC at both sampling times (Table 3 and Figure 5). However, inspection of protein abundance levels in all treatments reveals that down-regulation of this enzyme is more clear at T7, while at T28 results might be confounded by up-regulation of this protein in the SC treatment (Figure 5O). There was another sequence identified as an AChE-like isoform detected in the proteomic analysis, but its abundance level was not significantly 33 (p < 0.05) affected by TDCPP (sequence CL8744.Contig3_All_6, see supporting spreadsheet "Annex_data_pvalues"). Similar studies in different organisms have shown contrasting results. No evidence of AChE inhibitory activity by TDCPP was observed in zebrafish (Wang et al., 2015). Other OPFRs have been shown to have different effects on AChE activity and AChE mRNA expression in the Japanese medaka (Oryzias latipes), with TPHP inhibiting both responses, TNBP inducing them, and TBOEP and TCEP showing no significant effects (Sun et al., 2016). Similar results to the present ones were observed upon exposure of the clam Corbicula fluminea to a mixture of three OPFRs (TCEP, TPhP and TBP) in sediments, where inhibition of AChE activity, but not of AChE mRNA expression, was observed (Li et al., 2018). 4.4. Reported TDCPP concentrations causing biological effects Acute toxicity of TDCPP in aquatic organisms has been reported to occur at concentrations of mg L-1 (WHO, 1998), and developmental malformations have been reported at concentrations as low as 0.4 mg L-1 in zebrafish (McGee et al., 2012; Fu et al., 2013; Alzualde et al., 2018). In chronic studies with fish species, lower concentrations of TDCPP, between 4 and 40 µg L-1, have been shown to cause several deleterious effects, such as inhibition of female growth and fecundity (Wang et al., 2015; Zhu et al., 2015; Yu et al., 2017), endocrine disruption through altered steroidogenesis (Liu et al., 2012; Liu et al., 2013), neurotoxic effects (Wang et al., 2015) and hepatotoxicity and carbohydrates metabolism disruption (Hong et al., 2019). TDCPP at 6 µg L-1 also led to slight reductions in fecundity for the water flea Daphnia magna, and exposure to environmentally relevant concentrations of TDCPP (65 ng L-1 of TDCPP) caused slight (4%) reductions in growth of offspring (Li et al., 2015). 34 The concentration used in the present study (10 µg L-1) is within the lower range of concentrations were subtle effects in reproduction or neurotoxicity have been observed for other species. However, this concentration is one order of magnitude higher than measured concentrations of TDCPP in the aquatic environment, including seawater, with maximum reported concentrations being 378 ng L-1 in seawater of impacted coasts in East China (Hu et al., 2014), 855 ng L-1 in urban rivers and up to 1.4 µg L-1 in wastewater (Pantelaki and Voutsa, 2019). 4.5. Effect of exposure time on mussel gill proteome In the present experiment, the effects of exposure time on mussel gill proteome were by far much more important than the effects of the “chemical” exposure to TDCPP. Also Oliveira et al (2016), who exposed mussels to the antifouling biocide tralopyril during 30 days, observed that the variability in proteins abundance was mostly explained by the time of sampling rather than by different chemical treatments. This result is interesting because it provides evidence of a potential methodological limitation, yet largely ignored, in long term exposure experiments focussed on gene expression analyses, as far as it would be desirable to have minimal changes in proteins abundance in control organisms over time. In our case, 33% of the detected proteins were affected by exposure time (according to two-way ANOVA, see section 3.4.1), and enrichment analysis revealed that these changes were mainly related to a general decrease in energy metabolism and motility, and an increase in protein biosynthetic process. Decrease in energy metabolism and motility also occur during metabolic rate depression (MRD) observed in anoxia-tolerant or hibernating organisms (Storey and Storey, 2004). MRD or hypometabolism is characterized by decrease in ATP consumption that is attained by a lack of voluntary muscle movements, reduction in heart beat and 35 respiration rate, intermittent breathing, reduced kidney filtration rate, and slow down of synthesis and degradation of proteins and mRNA so that their half live is increased (Storey and Storey, 2004). This response has also been observed as a consequence of environmental stress caused by high and low temperature, oxygen deprivation, food restriction and water limitation (Storey and Storey, 1990). In our case, this might have been caused by the feeding regime used for mussels maintenance, consisting in feeding twice per week during a period of approximately 2 hours. A previous study using a similar feeding regime showed that a mild energetic disturbance is provoked in the mussels, and after 24 days with this regime, mussels showed a slightly but significantly lower condition index, with around 10% loss in tissue dry weight (Fernández-González et al., 2020). The increase in protein synthesis observed at T28 is however contrary to the expectations, given that protein synthesis is well known to be sensitive to the availability of energy and the rate of protein synthesis is suppressed during starvation and in response to stress (Storey and Storey, 2004; Larade and Storey, 2009). Indeed, González-Fernández et al (2017) showed that mussels fed a low (near-starving) feeding regime during a longer time (79 days) experienced a global suppression of transcriptional activity in gills (an indication of decreased protein synthesis as well). We can only hypothesize that despite being acclimated during one week to laboratory conditions, mussels were still in the process to adapt to these laboratory conditions, that include food restriction, during the course of our experiment, and for this they needed a new machinery operating in gills that required of an increased rate in protein production. For instance, it has been seen that mussels living in food restricted areas or acclimated to a low feeding regime in laboratory for a long time, show higher clearance rates (CR) than mussels living in food abundant areas/acclimated to high food 36 availability (Albentosa et al., 2012; Bellas et al., 2014; González-Fernández et al., 2015; González-Fernández et al., 2015). The increase in CR might involve not only protein production but also ATP consumption, that was not down-regulated with exposure time in our experiment, as occurred for ATP production (Table 2; see number of up- and down-regulated proteins for "ATPase activity" vs. "ATP biosynthesis process" and "mitochondrial/electron transport chain"). More research is needed in order to confirm whether the effects of exposure time on protein expression patterns, also observed by other researchers, are only due to the low and intermittent feeding regime usually employed in this type of experiments or if they are also caused by the general stress produced by captivity itself. In the first case, this could be solved or ameliorated by feeding the organisms in continuum or at least every day, requiring continuous-flow experiments or more frequent water changes in semistatic experiments. In any case, the detected effects of the contaminant/toxicant are still valid if compared with a control sampled at the same exposure time, although some effects might be masked, since most of the times the same processes are affected by both chemical and physical stressors (Groh and Suter, 2015). In fact, using metabolomics, Campillo et al (2019) showed that the effect of the toxicant fluoranthene was similar to that produced under conditions of food restriction in mussels, which masked the effect of the toxicant under these conditions, while when the feeding conditions became more favorable, the toxic effect of fluoranthene was more apparent. Also, physiological responses to pollutants can be affected by the energetic status of the mussels, and even defense mechanisms can be activated in response to low food stress, increasing mussels ability to cope with the pollutant, due to activation of the autophagic system and increased antioxidant protection (González-Fernández et al., 2015). 37 5. Conclusions The present study showed that TDCPP is quickly bioaccumulated and eliminated in mussels, reaching steady state already after 5 days of exposure. Tissue concentrations of TDCPP at equilibrium are around 150 times higher than those in the surrounding seawater. Detected effects in gills of mussels exposed to 10 µg L-1 of TDCPP (one order of magnitude higher than measured concentrations of TDCPP in the aquatic environment) were the inhibition of AChE activity and subtle proteomic changes, that give clues on the toxicity mechanisms (such as neurotoxicity through inhibition of receptor-type tyrosine-protein phosphatase N2-like) and detoxification of the compound (possibly involving Epimerase family protein SDR39U1). However, any strong effect of TDCPP on gill proteome can be discarded, at least under conditions of the present experiment. Methodological issues in toxicological proteomics were also highlighted. First of all, we observed that exposure time had widespread effects on the proteome analyzed, possibly related with stress produced by laboratory conditions, including intermittent feeding. Secondly, it was observed that changes in proteins abundance observed when comparing only the contaminant-exposed and solvent control treatments are not necessarily caused by the contaminant and can be caused by changes occurring in the solvent control treatment. This can only be detected if the data is observed in a broader context (including control individuals, or other exposure times). These methodological issues should be taken into account in the design of future similar studies. Stress factors that could mask proteomic changes produced by the assessed compound should be minimized as much as possible, and analysis of control organisms in addition to solvent control treatments, in combination with an adequate experimental replication, is recommended in order to be able to minimize and/or detect false positives. 38 Appendix A. Supplementary data The following is supplementary data to this article: File S1. Microsoft word file "Supporting_information.docx" including additional Tables S1 and S2 and additional Figures S1 to S3. File S2. Microsoft excel file "Annex_data_pvalues.xlsx" including two spreadsheets with normalized data (normalized reported ion intensities of TMT-labeled peptides) and p-values of the different two-way and one-way ANOVAs and Tukey HSD tests performed for each one of the 4191 proteins quantified. File S3. Text file "Annot_nrNCBI_May20.txt" including the complete annotation, with E-values and GO terms of the 4191 proteins detected. File S4. Compressed files "PEAKS_report.rar" with an html report exported from PEAKS Studio v8.0 software including extended information on protein identifications, their coverage, supporting peptides, etc. The mass spectrometry proteomics data (raw data and protein identifications) have been deposited to the ProteomeXchange Consortium via the PRIDE (Perez-Riverol et al., 2019) partner repository with the dataset identifier PXD019720 and 10.6019/PXD019720. Acknowledgements Exposure experiments were made in the ECIMAT (Universidad de Vigo). LC-MS/MS proteomic analyses were performed in the CACTI (Universidad de Vigo) by Manuel 39 Marcos and histological sample processing were performed in the ECIMAT by Rosana Rodríguez. We thank Pilar Feijóo for technical assistance provided during exposure experiments and Valentín Trujillo for valuable help in creating R scripts. L.E. F.-G. is supported with a predoctoral fellowship from Xunta de Galicia (Consellería de Cultura, Educación e Ordenación Universitaria grant ED481A-2017/298) partially co-founded by operative program FSE Galicia 2014-2020. M.C. is supported by the Ramón y Cajal contract (Ministry of Science and Innovation of Spain). This study was funded by Xunta de Galicia (refs. ED431C 2017/36, ED431C 2017/46, ED431C 2020/05 and IN607B 2018/10) and the Spanish Agencia Estatal de Investigación (refs. CTM2017- 84763-C3-R-2 and CTM2016-77945-C3), partially co-funded by the European Regional Development Fund (FEDER/ERDF). 40 References Albentosa, M., Viñas, L., Besada, V., Franco, A. and González-Quijano, A., 2012. First measurements of the scope for growth (SFG) in mussels from a large scale survey in the North-Atlantic Spanish coast. Sci. Tot. Environ. 435-436, 430-445. Alzualde, A., Behl, M., Sipes, N.S., Hsieh, J.-H., Alday, A., Tice, R.R., Paules, R.S., Muriana, A. and Quevedo, C., 2018. Toxicity profiling of flame retardants in zebrafish embryos using a battery of assays for developmental toxicity, neurotoxicity, cardiotoxicity and hepatotoxicity toward human relevance. Neurotoxicol. Teratol. 70, 40-50. Apraiz, I., Mi, J. and Cristobal, S., 2006. Identification of proteomic signatures of exposure to marine pollutants in mussels (Mytilus edulis). Mol. Cell. Proteomics 5, 1274-1285. Atasaral-Şahin, Ş., Romero, M.R., Cueto, R., González-Lavín, N., Marcos, M. and Diz, A.P., 2015. Subtle tissue and sex-dependent proteome variation in mussel (Mytilus galloprovincialis) populations of the Galician coast (NW Spain) raised in a common environment. Proteomics 15, 3993-4006. ATSDR (2009). Draft Toxicological Profile for Phosphate Ester Flame Retardants (September). United States Department of Health and Human Services. Aznar-Alemany, Ò., Aminot, Y., Vilà-Cano, J., Köck-Schulmeyer, M., Readman, J.W., Marques, A., Godinho, L., Botteon, E., Ferrari, F., Boti, V., Albanis, T., Eljarrat, E. and Barceló, D., 2018. Halogenated and organophosphorus flame retardants in European aquaculture samples. Sci. Tot. Environ. 612, 492-500. Babin, P.J., Bogerd, J., Kooiman, F.P., Van Marrewijk, W.J.A. and Van der Horst, D.J., 1999. Apolipophorin II/I, apolipoprotein B, vitellogenin, and microsomal triglyceride transfer protein genes are derived from a common ancestor. J. Mol. Evol. 49, 150-160. Beiras, R., 2018. Marine Pollution, Sources, Fate and Effects of Pollutants in Coastal Ecosystems. Elsevier, Amsterdam. Beliaeff, B. and Burgeot, T., 2002. Integrated biomarker response: A useful tool for ecological risk assessment. Environ. Toxicol. Chem. 21, 1316-1322. Bellas, J., Albentosa, M., Vidal-Liñán, L., Besada, V., Franco, M.T., Fumega, J., González-Quijano, A., Viñas, L. and Beiras, R., 2014. Combined use of chemical, biochemical and physiological variables in mussels for the assessment of marine pollution along the N-NW Spanish coast. Mar. Env. Res. 96, 105-117. Betts, K.S., 2008. New thinking on flame retardants. Environ. Health Persp. 116, A210- 213. Beyer, J., Green, N.W., Brooks, S., Allan, I.J., Ruus, A., Gomes, T., Bråte, I.L.N. and Schøyen, M., 2017. Blue mussels (Mytilus edulis spp.) as sentinel organisms in coastal pollution monitoring: A review. Mar. Env. Res. 130, 338-365. 41 Bocquené, G. and Galgani, F., 1998. Biological effects of contaminants: Cholinesterase inhibition by organophosphate and carbamate compounds. ICES Techniques in Marine Environmental Sciences 22, 1-12. Bradford, M.M., 1976. A rapid and sensitive method for the quantitation of microgram quantities of protein utilizing the principle of protein-dye binding. Anal. Biochem. 72, 248-254. Campillo, J.A., Sevilla, A., González-Fernández, C., Bellas, J., Bernal, C., Cánovas, M. and Albentosa, M., 2019. Metabolomic responses of mussel Mytilus galloprovincialis to fluoranthene exposure under different nutritive conditions. Mar. Env. Res. 144, 194-202. Campone, L., Piccinelli, A.L., Östman, C. and Rastrelli, L., 2010. Determination of organophosphorous flame retardants in fish tissues by matrix solid-phase dispersion and gas chromatography. Anal. Bioanal. Chem. 397, 799-806. Campos, A., Tedesco, S., Vasconcelos, V. and Cristobal, S., 2012. Proteomic research in bivalves: Towards the identification of molecular markers of aquatic pollution. J. Proteomics 75, 4346-4359. Campos, A., Danielsson, G., Farinha, A.P., Kuruvilla, J., Warholm, P. and Cristobal, S., 2016. Shotgun proteomics to unravel marine mussel (Mytilus edulis) response to long-term exposure to low salinity and propranolol in a Baltic Sea microcosm. J. Proteomics 137, 97-106. Carvajal-Rodriguez, A. and de Uña-Alvarez, J., 2011. Assessing significance in highthroughput experiments by sequential goodness of fit and q-value estimation. PLoS ONE 6, e24700. Davies, I.M. and Vethaak, A.D., 2012. Integrated marine environmental monitoring of chemicals and their effects. ICES Cooperative Research Report No. 315. Dishaw, L.V., Powers, C.M., Ryde, I.T., Roberts, S.C., Seidler, F.J., Slotkin, T.A. and Stapleton, H.M., 2011. Is the PentaBDE replacement, tris (1,3-dichloro-2- propyl) phosphate (TDCPP), a developmental neurotoxicant? Studies in PC12 cells. Toxicol. Appl. Pharm. 256, 281-289. Diz, A.P., Carvajal-Rodríguez, A. and Skibinski, D.O.F., 2011. Multiple hypothesis testing in proteomics: A strategy for experimental work. Mol. Cell. Proteomics 10. Duroudier, N., Cardoso, C., Mehennaoui, K., Mikolaczyk, M., Schäfer, J., Gutleb, A.C., Giamberini, L., Bebianno, M.J., Bilbao, E. and Cajaraville, M.P., 2019. Changes in protein expression in mussels Mytilus galloprovincialis dietarily exposed to PVP/PEI coated silver nanoparticles at different seasons. Aquat. Toxicol. 210, 56-68. EU (2008). European Union risk assessment report: Tris[2-chloro-1- (chloromethyl)ethyl] phosphate (TDCP). European Communities. CAS No. 13674-87-8, EINECS No. 237-159-2. Luxembourg. 48 FIGURE CAPTIONS Figure 1. Bioaccumulation of TDCPP in mussels exposed during 28 days to 10 µg L-1 TDCPP (exposure phase, black dots) and to control seawater during additional 14 days (depuration phase, white dots). Curves represent the fitting of the data to equations 1 and 2, with fitted parameters ku = 789 ± 83 L kgdw-1 d-1 and ke = 1.07 ± 0.11 d-1. Figure 2. Hierarchical cluster with heatmap (a) and PCA (b) analysis of 1148 proteins that were differentially expressed (p < 0.05) in mussel gills at any treatment. Both rows and columns were clustered using Euclidean distance and average linkage methods. PCA method was SVD (singular value decomposition) with iterative imputation for missing value estimation (missing values represented 1.8%). T0 (black circles) = mussels collected at the beginning of the experiment; T7-TDCPP and T28-TDCPP (black and white triangles, respectively) = mussels exposed during 7 and 28 days respectively to 10 µg L-1 of TDCPP; T7-SC and T28-SC (black and white squares, respectively) = mussels exposed during 7 and 28 days to the solvent control; T28- Control (white circles) = mussels exposed during 28 days to control seawater. Replicates 1, 2, 3 correspond to TMTsixplex proteomic experiments 1, 2 and 3, respectively. Figure 3. Volcano plot of differences in protein abundance in mussel gills according to factor "time" (28 days compared to 0 and 7 days of exposure time, all chemical treatments combined). FC = fold change (T28 / T0-7). Horizontal lines represent pvalues of 0.01 (short dash) and 0.05 (long dash). Figure 4. Volcano plots of differences in protein abundance in pairwise comparisons of TDCPP vs. Solvent Control (SC) treatments at 7 days (T7) or 28 days (T28) of exposure time. FC = fold change. p-value was obtained from post-hoc Tukey HSD test after 49 ANOVA. Horizontal lines correspond to p-values of 0.01 (short dash) and 0.05 (long dash). Figure 5. Relative protein abundance (based on normalized reported ion intensities of TMT-labeled peptides) of proteins differentially expressed in mussel gills in TDCPP treatments vs. SC after 7 days of exposure (p < 0.01, 8 proteins), after 28 days of exposure (p < 0.01, 6 proteins), or at both exposure times (p < 0.05, 3 proteins). Mean ± SD (n = 3 pools of 4 individuals) is represented, with the exception of T28-Control (where n = 2 pools of 4 individuals). Asterisks above a bar denote significant differences at the p < 0.05 level (*) or at the p < 0.01 level (**) compared to the solvent control at the same sampling time. Figure 6. AChE and GST activity in gills of mussels at the beginning of the experiment (T0) and after being exposed during 7 (T7) or 28 (T28) days to control seawater, an acetone solvent control (SC), or 10 µg L-1 TDCPP (TDCPP). Mean ± SD (n = 12) is represented. Asterisks denote significant differences at the p < 0.05 level (*) or at the p < 0.01 level (**) of TDCPP treatments in comparison to SC at the same sampling time. 50 TABLES Table 1. Set up of the three proteomic experiments (Exp.) labeled with TMTsixplexTM (Tags: TMT6-126, TMT6-127, TMT6-128, TMT6- 129, TMT6-130 and TMT6-131). Each pool consisted in combining the same protein amount of four gills corresponding to two males and two females exposed to the same treatment (Control = filtered seawater, SC = solvent control, TDCPP = 10 µg L-1 TDCPP) during a given time (T0, T7 and T28 correspond to 0, 7, and 28 days of exposure, respectively). T28-Control_pool 1bis corresponds to a technical replicate of T28-Control_pool 1. TMT6-126 TMT6-127 TMT6-128 TMT6-129 TMT6-130 TMT6-131 Exp. 1 T0_pool 1 T7-SC_pool 1 T7-TDCPP_pool 1 T28-Control_pool 1 T28-SC_pool 1 T28-TDCPP_pool 1 Exp. 2 T0_pool 2 T7-SC_pool 2 T7-TDCPP_pool 2 T28-Control_pool 2 T28-SC_pool 2 T28-TDCPP_pool 2 Exp. 3 T0_pool 3 T7-SC_pool 3 T7-TDCPP_pool 3 T28-Control_pool 1bis T28-SC_pool 3 T28-TDCPP_pool 3 51 Table 2. GO (Gene Ontology) enrichment analysis of the proteins affected by exposure time. Test set corresponds to the list of proteins that were differentially expressed (p < 0.05) in the gill proteome of mussels exposed during 28 days in comparison with those that were exposed during 0 or 7 days (regardless treatment applied). Reference set corresponds to all other non-differentially expressed proteins in the dataset. The number (and percentage) of proteins related with each GO Term that were significantly up- or down-regulated in the test set is also shown. FDR = False discovery rate. GO Category GO ID GO Term Related process Test set Reference set FDR up-regulated (T28/T0-7) down-regulated (T28/T0-7) BIOLOGICAL PROCESS GO:0006412 translation Protein synthesis 5.8% 2.6% 1.8E-04 81 (96%) 3 (4%) GO:0034660 ncRNA metabolic process 1.8% 0.5% 6.1E-03 26 (100%) 0 GO:0022618 ribonucleoprotein complex assembly 2.1% 0.7% 1.2E-02 29 (97%) 1 (3%) GO:0022900 electron transport chain Energy metabolism 1.86% 0.55% 9.3E-03 4 (15%) 23 (85%) GO:0006754 ATP biosynthetic process 0.62% 0.04% 3.1E-02 0 9 (100%) GO:0070286 axonemal dynein complex assembly Motility 0.62% 0.04% 3.1E-02 0 9 (100%) CELLULAR COMPONENT GO:0005930 axoneme Motility 1.8% 0.5% 1.3E-02 1 (4%) 25 (96%) GO:0030286 dynein complex 4.1% 2.2% 2.8E-02 4 (7%) 56 (93%) GO:0005840 ribosome Protein synthesis 3.0% 1.3% 1.6E-02 42 (95%) 2 (5%) GO:0005747 mitochondrial respiratory chain complex I Energy metabolism 0.83% 0.11% 2.8E-02 0 (0%) 12 (100%) MOLECULAR FUNCTION GO:0042623 ATPase activity, coupled Energy metabolism 2.6% 0.9% 6.1E-03 20 (54%) 17 (46%) GO:0003735 structural constituent of ribosome Protein synthesis 2.8% 1.1% 6.8E-03 40 (98%) 1 (2%) GO:0019843 rRNA binding 0.76% 0.04% 8.9E-03 11 (100%) 0 GO:0035639 purine ribonucleoside triphosphate binding 14% 11% 1.1E-02 134 (66%) 70 (34%) GO:0032555 purine ribonucleotide binding 14% 10% 3.0E-02 134 (66%) 70 (34%) GO:0003777 microtubule motor activity Motility 4.1% 2.2% 3.6E-02 5 (8%) 54 (92%) 52 Table 3. Differentially expressed (DE) proteins in mussels exposed to TDCPP in comparison with those exposed to solvent control after (A) 7 or (B) 28 days of exposure (p < 0.01), and (C) after 7 and 28 days of exposure (p < 0.05). FC = fold change. The FC and p-value for these (TDCPP/SC) and other (TDCPP/Control; SC/Control) pairwise comparisons and for the factor time are also shown. T28 = 28 days of exposure; T0-7 = 0 and 7 days of exposure. Control = filtered seawater. SC = solvent control. TDCPP = 10 µg L-1 TDCPP. Protein ID corresponds to the given in the transcriptomic database in Moreira et al. (2015), with an added number indicating the translation frame, protein annotation according to NCBI Blastp search (May 2020), overall function according to UniProt. Cell color legend: light grey (pvalue < 0.05), dark grey (p-vaue < 0.01), red (up-regulated proteins, FC > 1), blue (down-regulated proteins, FC < 1). 53 Protein ID T7-TDCPP / T7- SC T28-TDCPP / T28-SC T28-TDCPP / T28-Control T28-SC / T28- Control T28 / T0-7 Annotation Overall function p-value FC pvalue FC p-value FC p-value FC p-value FC A) DE proteins in T7-TDCPP vs. T7-SC CL716.Contig2_All_1 0.002 1.45 0.489 1.13 0.937 1.03 0.975 0.92 0.771 1.02 Apolipophorins-like Lipid transport Unigene19196_All_4 0.007 1.10 0.559 0.97 0.005 0.89 0.051 0.92 0.044 1.07 Mitochondrial transcription factor A Mitochondrial transcription regulation CL12088.Contig1_All_5 0.006 1.09 0.292 1.04 0.998 1.01 0.629 0.97 0.566 0.99 Toll-interacting protein Negative regulation of immune response CL9887.Contig2_All_3 0.004 0.93 1.000 1.00 0.449 1.03 0.372 1.03 0.044 1.03 Uncharacterized transmembrane protein -- Unigene42653_All_4 0.008 0.91 0.120 1.06 0.998 0.98 0.104 0.93 0.403 0.98 Axonemal dynein heavy chain 5 Force generating protein of cilia CL9603.Contig1_All_5 0.002 0.72 0.116 0.88 0.559 1.07 0.019 1.22 0.957 1.00 Adenosine deaminase Positive regulation of immune response, Cellular signaling events CL4659.Contig1_All_1 0.002 0.67 0.213 0.86 0.049 0.81 0.763 0.93 0.183 0.87 Axonemal dynein heavy chain Force generating protein of cilia CL9301.Contig2_All_2 0.009 0.62 0.980 0.95 1.000 0.98 0.998 1.03 0.666 0.95 Receptor-type tyrosine-protein phosphatase N2-like Required in vesicle-mediated secretory processes B) DE proteins in T28-TDCPP vs. T28-SC Unigene96059_All_5 0.877 0.94 0.009 1.28 0.013 1.32 1.000 1.04 0.003 0.81 Cilia- and flagella-associated protein 47 Component of motile cilium CL4722.Contig2_All_3 0.992 0.99 <0.001 1.21 0.149 1.11 0.057 0.92 0.002 1.13 Protein FAM114A2-like Purine nucleotide binding CL10240.Contig2_All_6 0.182 0.91 0.008 1.19 0.116 1.16 0.843 0.98 0.457 0.96 Plasminogen activator inhibitor 1 Serine protease inhibitor, regulator of cell migration Unigene73717_All_2 0.857 1.10 0.006 0.69 0.115 0.87 0.761 1.26 0.007 1.25 -- -- Unigene5922_All_2 0.988 1.08 <0.001 0.51 1.000 1.02 <0.001 2.07 0.017 1.48 Ribosomal protein S6 kinase beta-1 Promote cell proliferation, cell growth and cell cycle progression Unigene6637_All_5 0.921 1.33 0.005 0.13 0.996 0.82 0.016 6.11 0.045 4.09 -- -- C) DE proteins in T7-TDCPP vs. T7-SC and T28-TDCPP vs. T28-SC Unigene40987_All_5 <0.001 0.84 0.006 0.88 0.297 1.05 0.001 1.19 0.002 1.22 Acetylcholinesterase-like Terminate neurotransmission by rapidly hydrolyzing acetylcholine CL15287.Contig2_All_4 0.028 0.90 0.021 1.11 0.041 1.092 1.000 0.982 0.244 1.04 DNA topoisomerase II Control of topological states of DNA, Chromosome segregation Unigene31732_All_1 0.014 1.23 0.034 1.20 0.063 1.177 1.000 0.989 0.033 0.89 Epimerase family protein SDR39U1 Putative NADP-dependent oxidoreductase 54 Figures Figure 1. Bioaccumulation of TDCPP in mussels exposed during 28 days to 10 µg L-1 TDCPP (exposure phase, black dots) and to control seawater during additional 14 days (depuration phase, white dots). Curves represent the fitting of the data to equations 1 and 2, with fitted parameters ku = 789 ± 83 L kgdw-1 d-1 and ke = 1.07 ± 0.11 d-1. Exposure time (days) 010 20 30 40 TDCPP bioaccumulation (ng / g dw ) 0 2000 4000 6000 8000 55 Figure 2. Hierarchical cluster with heatmap (a) and PCA (b) analysis of 1148 proteins that were differentially expressed (p < 0.05) in mussel gills at any treatment. Both rows and columns were clustered using Euclidean distance and average linkage methods. PCA method was SVD (singular value decomposition) with iterative imputation for missing value estimation (missing values represented 1.8%). T0 (black circles) = mussels collected at the beginning of the experiment; T7-TDCPP and T28-TDCPP (black and white triangles, respectively) = mussels exposed during 7 and 28 days respectively to 10 µg L-1 of TDCPP; T7-SC and T28-SC (black and white squares, respectively) = mussels exposed during 7 and 28 days to the solvent control; T28- Control (white circles) = mussels exposed during 28 days to control seawater. Replicates 1, 2, 3 correspond to TMTsixplex proteomic experiments 1, 2 and 3, respectively. 56 Figure 3. Volcano plot of differences in protein abundance in mussel gills according to factor "time" (28 days compared to 0 and 7 days of exposure time, all chemical treatments combined). FC = fold change (T28 / T0-7). Horizontal lines represent pvalues of 0.01 (short dash) and 0.05 (long dash). 57 Figure 4. Volcano plots of differences in protein abundance in pairwise comparisons of TDCPP vs. Solvent Control (SC) treatments at 7 days (T7) or 28 days (T28) of exposure time. FC = fold change. p-value was obtained from post-hoc Tukey HSD test after ANOVA. Horizontal lines correspond to p-values of 0.01 (short dash) and 0.05 (long dash).