Assessing the impact of Benzo[a]pyrene on Marine Mussels: Application of a novel targeted low density microarray complementing classical biomarker responses

Despite the increasing use of mussels in environmental monitoring and ecotoxicological studies, their genomes and gene functions have not been thoroughly explored. Several cDNA microarrays were recently proposed for Mytilus spp., but putatively identified partial transcripts have rendered the generation of robust transcriptional responses difficult in terms of pathway identification. We developed a new low density oligonucleotide microarray with 465 probes covering the same number of genes. Target genes were selected to cover most of the well-known biological processes in the stress response documented over the last decade in bivalve species at the cellular and tissue levels. Our new ‘STressREsponse Microarray’ (STREM) platform consists of eight sub-arrays with three replicates for each target in each sub-array. To assess the potential use of the new array, we tested the effect of the ubiquitous environmental pollutant benzo[a]pyrene (B[a]P) at 5, 50, and 100 μg/L on two target tissues, the gills and digestive gland, of Mytilus galloprovincialis exposed invivo for three days. Bioaccumulation of B[a]P was also determined demonstrating exposure in both tissues. In addition to the well-known effects of B[a]P on DNA metabolism and oxidative stress, the new array data provided clues about the implication of other biological processes, such as cytoskeleton, immune response, adhesion to substrate, and mitochondrial activities. Transcriptional data were confirmed using qRT-PCR. We further investigated cellular functions and possible alterations related to biological processes highlighted by the microarray data using oxidative stress biomarkers (Lipofuscin content) and the assessment of genotoxicity. DNA damage, as measured by the alkaline comet assay, increased as a function of dose.DNA adducts measurements using 32P-postlabeling method also showed the presence of bulky DNA adducts (i.e. dG-N2-BPDE). Lipofiscin content increased significantly in B[a]P exposed mussels. Immunohistochemical analysis of tubulin and actin showed changes in cytoskeleton organisation. Our results adopting an integrated approach confirmed that the combination of newly developed transcriptomic approcah, classical biomarkers along with chemical analysis of water and tissue samples should be considered for environmental bioimonitoring and ecotoxicological studies to obtain holistic information to assess the impact of contaminants on the biota.


Introduction
In the last two decades, genomics-based approaches have provided insights into the mechanisms underlying the biological responses in natural biota following exposures to environmental pollutants. Recent developments in 'omics' technologies have raised our expectations to elucidate crucial biological processes following exposure to pollutants to which biota are particularly exposed chronically at low dose/ concentrations. Notably, transcriptomics, the simultaneous measurement of thousands of mRNAs in a biological sample, has been a dependable tool, increasing our understanding of many important physiological processes in marine biosystems [1,2,3,4,5]. The accessibility of new genomic resources, high-throughput molecular technologies such as microarray, pyrosequencing and next generation sequencing have provided information and helped to identify molecular targets of pollutants (e.g. genes) that contribute to stress in natural populations [6,7].
Marine mussels belonging to the genus Mytilus are commonly found in the coastal areas of the Mediterranean and Atlantic Ocean and play a significant role in coastal ecology apart from being economically important for shellfish industries [8]. Mytilus spp. have been used extensively in biomonitoring projects through the application of a battery of physiological and cellular markers, providing evidence of a stress syndrome and demonstrating the biological risk associated with polluted environments [9,10]. Mussels are particularly useful in this context because they inhabit regions of different pollution status, accumulate xenobiotics, and are sessile. However, despite the increasing use of these organisms in an ecotoxicological context, their genomes and gene functions have not been thoroughly investigated. Although several cDNA microarrays were developed recently in Mytilus galloprovincialis, Mytilus edulis, and Mytilus trossulus, putatively identified partial transcripts limit their use in generating reliable gene expression profiles in terms of pathway identification [11].
Benzo[a]pyrene (B[a]P) is a known genotoxic and carcinogenic agent and a ubiquitous environmental contaminant. It is one of the US EPA's priority pollutants. As one of the most commonly tested compounds in ecotoxicology [12,13], B[a]P has been shown to produce toxic effects in aquatic organisms across a range of biological levels. For example, B[a]P induces xenobiotic metabolism in fish and marine polychaetes [14,15], causes developmental abnormalities and genotoxic effects in the embryos of oysters, annelids and zebrafish [16,17,18]. B[a] P also has the ability to form genotoxic intermediates, such as B[a]P-7,8-dihydrodiol-9,10-epoxide (BPDE), which form bulky DNA adducts that may result in mutation [19,20], with potential impacts on the individual organism or population [21,22]. B[a]P has also been reported to induce genotoxicity, oxidative stress [12,23,24] and endocrine disruption [25] in marine bivalves. In addition, B[a]P has been shown to interact with other pollutants (e.g. engineered nanoparticles) that could alter the potential toxic effects [13,26,27]. As such, mixture effects become more prominent in ecotoxicology and environmental protection [28].A better understanding of the mechanisms underlying the toxicity of this ubiquitous contaminant would allow us to assess its role in more realistic, multi-component toxicity.
In the backdrop of the above information, our study adopts an integrated approach, presenting for the first time, a new low density oligonucleotide microarray in mussel species complementing other novel biomarker responses in an ecotoxicological context. Target genes used in the new platform were selected to cover most of the well-known biological processes including covering stress related responses in bivalve species(e.g. lysosomal activity, DNA metabolism, translation, mitochondrial activity, carbohydrate metabolism, heat shock response, oxidative stress, proteolysis, phase I and II xenobiotic metabolism, metal resistance, endocrine disruption, immune response, pumps, cytoskeleton, and adhesion to substrate).The new platform was used to evaluate the biological effects of exposure to three concentrations of B[a]P on the gills and digestive gland (DG) of M. galloprovincialis. Q-RT-PCR analysis was performed to confirm the over-represented biological processes in response to B[a]P exposure (i.e., DNA metabolism, oxidative stress, and cytoskeleton). An important aspect of this study was to investigate the consequences of transcriptomic changes depicted by microarray data on relevant biological processes in gills and DG at the cellular/tissue level. To complement the microarray and transcriptomics responses, genotoxic responses were assessed by determining DNA strand breaks using the comet assay and determination of bulky DNA adduct formation by 32 P-postlabelling method. The presence of oxidative stress was evaluated by lysosomal lipofuscin accumulation and the alteration of cytoskeletal structures by immuno-histochemistry. Concentration of B[a]P in water and tissue samples were also determined using analytical chemistry techniques.

Chemicals and organisms
Unless otherwise stated, all chemicals were purchased from Sigma-Aldrich (UK/Italy). Adult M. galloprovincialis (50.7 ± 2.8 mm) were collected from the intertidal zone at Trebarwith Strand, Cornwall, UK (50˚38' 40" N, 4˚45' 44" W) in October 2014. No specific permission is required to collect marine mussels which are abundant at this collection site in Cornwall, UK. Collection of the limited number of mussels for laboratory based experiments does not have any detrimental effects on the environment at the site and does not involve any ethical issues. It does not require any licence from any authorities in the UK as they are not protected or endangered species. This is in line with our previous studies using this species collected from this site [12,13].This site has been used previously as a reference location for ecotoxicological studies [31]. The site is relatively free of disease and is remotely located [32]. Mussels were transported back to the laboratory in cool boxes and allowed to depurate for a minimum of 1 week in natural seawater from Plymouth Sound. During this period, mussels were fed a suspension of Isochrysis galbana every 3 days (1.05 × 10 6 cells/mL), with a 100% water change 2 hours after each feeding [33].

Experimental design and sampling
After depuration, the mussels were transferred to 2-L glass beakers containing 1.8 L of the same seawater as above and allowed to acclimatize for 48 h. The experiment began after this period and consisted of a 3-day static exposure with no water changes, during which the mussels were not fed. Two mussels were used per beaker. A photoperiod of 12h light: 12h dark was maintained throughout the experiment. Good seawater oxygenation was provided by a bubbling system. Seawater quality was monitored in each of the beakers by measuring salinity (35.4 ± 0.09‰), pH (7.9 ± 0.06), %dissolved oxygen (97.9 ± 3.22%) and temperature (15.3 ± 0.68˚C).Groups of mussels were exposed to five treatments: a solvent control (0.02% dimethyl sulfoxide [DMSO]; 36 mussels); 5 μg/L B[a]P (36 mussels); 50 μg/L B[a]P (36 mussels); 100 μg/L B[a]P (36 mussels); and 1000 μg/L B[a]P (9 mussels). After 3 days, mussel sex was determined by mantle smear and light microscopy. Females were processed for transcriptomics (microarray and qRT-PCR) and immunohistochemistry (DG only). Individuals of both sexes were used for the assessment of genotoxicity (i.e. DNA damage measured by the comet assay or bulky DNA adducts by 32 P-postlabelling; gills and DG). Tissue samples were also taken from the DG and gills for chemical analysis of the B[a]P concentrations using gas chromatography-mass spectrometry (GC-MS). Samples were transferred into pre-weighed glass vials and frozen for storage awaiting analysis.

Determination of B[a]P concentrations in water and tissues by GC-MS
Water samples (9 mL) were collected in glass vials and dichloromethane (1 mL, HPLC grade, Rathburn Chemicals Ltd., UK) was added before storage in the dark at 4˚C [13]. Prior to analysis, B[a]P d12 (15 μL of a 3 mg/L solution in acetone) was spiked as an internal standard. Following thorough shaking, the dichloromethane layers were removed into glass injection vials.
Before analysing the tissue samples(Gills and DG), water was removed by freeze drying (Christ Alpha LDplus freeze-dryer), typically for 24 h, and the water content determined by differential weighing. Dried tissue was ground in the vials using a razor blade and was spiked with 50 μL of the internal standard. Dichloromethane (2 mL) was added and the tissues extracted in a sonication bath for 20 min. The extracts were filtered through pre-cleaned glass wool. The lipid content was determined by sub-sampling a weighed aliquot of 500 μL and drying it in a nitrogen evaporator. The dry weight was considered to be the lipid content [34].
Both water and tissue extracts were analysed using an Agilent Technologies 7890A Gas Chromatography (GC) system interfaced with an Agilent 5975 series Mass Selective(MS) detector. A Restek Rxi-1MS (crosslinked dimethyl polysiloxane) capillary column (30 m) with a film thickness of 0.25 μm and internal diameter of 0.25 mm was used for separation, with helium as a carrier gas (maintained at a constant flow rate of 1 mL/min). Extracts were injected splitless, with the injector maintained at 300˚C. The oven temperature was 40˚C for 1 minute and then increased at 15˚C/minute to a final temperature of 300˚C, which was held for 4 minutes. The mass spectrometer was operated in electron impact mode at 70 eV with the ion source and quadrupole analyzer temperatures fixed at 230˚C and 150˚C, respectively. Samples were screened for B[a]P and B[a]P d12 using selected ion monitoring in which the target ions were 252 (253 and 126 for confirmation purposes) and 264, respectively. Within each batch of samples, multiple solvent blanks and standard calibrating mixtures were run regularly to ensure the stability of the system performance. Procedural blanks were also included in the sequence for quality assurance purposes. The B[a]P concentration was calculated based on the internal standard. Often, due to their high lipid content, mussel extracts can require a preliminary sample clean-up prior to GC-MS. In our case, the B[a]P signal intensities for the exposed mussels were sufficiently high to avoid this additional step.

Microarray hybridization and analysis
Competitive dual-color microarray hybridization was performed with the new STREM platform; fluorescence-labeled cDNA probes were obtained by direct labeling in the presence of modified Cy3-and Cy5-dCTP (Perkin Elmer). The procedure was carried out as described previously [3] using 0.5 μg of an anchored oligodT (19)VN. Total RNA was extracted from pieces of individual tissues(gills and DG) from females using acid phenol-chloroform precipitation according to [35] with TRI-Reagent (Sigma-Aldrich).Four biological replications were considered for the microarray and for the qRT-PCR analysis. The RNA was further purified by precipitation in the presence of 1.5 M LiCl 2 , and the quality of each RNA preparation was confirmed by UV spectroscopy and TBE agarose gel electrophoresis in the presence of formamide [3]. Laser scanning of microarrays was performed with an Agilent G2565CA scanner (Agilent Technologies, Inc., USA) at 5-μm resolution. Sixteen-bit TIFF images were analyzed using Genepix 6.0 (Axon) to extract raw fluorescence data from each spot.
We performed global mean normalization across microarray surfaces and local mean normalization across element signal intensity. Log2 transformation was performed for each expression level. Each experimental condition was evaluated as 12 separate values (3 technical and 4 biological replicates). Differentially expressed genes were identified by Significance Analysis of Microarray (SAM, http://statweb.stanford.edu/~tibs/SAM/). Gene expression was considered to be significantly different in the test condition versus the reference condition when the log value was > 0.7.
The experimental design accounted for a complete "triangular loop" in which each RNA sample from the tissue taken from mussels exposed to the three B[a]P conditions was hybridized with RNA from control mussels (DMSO-exposed mussels). Each experimental condition had at least four biological replicates from individual female animals using the day-swap procedure, for a total of 24 experiments. MIAMI-compliant microarray data, including a detailed description of the experimental design and each hybridization experiment, were deposited in the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/query/) under identifier GSE84605.

Functional genomics analysis
Functional characterization of mussel genes represented on a microarray was based on Gene Ontology (GO) annotation and carried out in Blast2GO [36] using default parameters. However, in the case of the STREM platform, target genes were putatively annotated and ranked under established biological processes, making the generation of robust processes easier and faster.

qRT-PCR
The qRT-PCR was carried out using the same RNA extract as microarray hybridization. Relative mRNA abundance of the mussel genes encoding 10 Probes and primer pairs (S1 Table) were designed using Beacon Designer v3.0 (Premier Biosoft International, Inc.). All primers and dual-labeled TaqMan probes were synthesized by MWG-Biotech GmbH (Germany).
Relative expression data were geometrically normalized to 18S rRNA (L33452), an invariant actin isotype (AJ625116), and ribosomal protein riboL27 (AJ625928), which were selected from a list of genes for which expression does not vary over more than 50 conditions, including toxic treatments, stages of the life cycle, and various tissues [3].
A specific duplex TaqMan assay was developed to amplify 0.25 ng of RNA reverse-transcribed to cDNA in the presence of 0.1 μmol/L of each dual-labeled probe (hexachlorofluorescein/BH1 for actin and Texas Red/BH2 for 18S rRNA) and 0.1 μmol/L and 0.4 μmol/L of forward and reverse primer, respectively, for 18S rRNA and actin (S1 Table). For all TaqMan assays, the reaction occurred as follows: 30 seconds at 95˚C, followed by 40 cycles of 10 seconds at 95˚C and 20 seconds at 60˚C. The qRT-PCR was performed with four biological replicates and three technical replicates. For the assays P53,Topoisomerase and DNA ligase, the thermal protocol was as presented by [3]. Statistical analyses were carried out on the group mean values using a random reallocation test [37].

Determination of DNA strand breaks using the comet assay
Gill and DG tissues were stored on ice prior to processing for the comet assay. After rough chopping, each piece of tissue was added to 1 mL of 1.6 mg/mL dispase II solution (in Hank's buffered saline) pre-warmed to 37˚C. After digesting in the dark for 30 min, the resulting suspension was coarsely filtered through gauze and spun at 200 g to remove any debris. The supernatant was checked for cell viability using trypan blue (0.04%), and only samples with >90% unstained cells were used. A sub-sample of 100 μL of the cell suspension was pelleted at 350 g and 180 μL low melting point agarose added. Two replicate microgels (75 μL) were formed by coverslipping and allowed to set at 4˚C for a minimum of 15 min. After the cover slips were removed, cells were lysed in 2.5 M NaCl, 100 mM EDTA, 10 mMTris, 1% N-laurylsarcosine, 1% Triton X-100, and 10% DMSO (pH adjusted to 10 with NaOH) for 1 h at 4˚C in the dark. Slides were then transferred to an electrophoresis chamber containing 1 mM EDTA and 0.3 M NaOH (pH 13). DNA was allowed to unwind in the dark for 20 min, followed by electrophoresis at 1 V/cm (*300 mA) for 20 min in the dark. Gels were preserved with ice cold 100% ethanol and then scored after the addition of 20 μg/mL ethidium bromide using a Leica epifluorescence microscope (Leica Microsystems, Milton Keynes, UK) and Comet 4 image analysis system (Perceptive Instruments, Bury St Edmunds, UK). One hundred nucleoids were assessed per slide; all samples were measured blind. Tail intensity (% tail DNA), defined as the percentage of DNA migrated from the head of the comet into the tail, was used as a measure of DNA damage induced.

Lysosomal lipofuscin content
Lysosomal lipofuscin content in the cells of the digestive glands was estimated on cryostat sections (10 μm) obtained as previously reported [42]. The lipofuscin content was assessed using the Schmorl reaction [43]. Lipofuscin accumulation was quantified by image analysis as described above and expressed as a percentage variation with respect to controls. Data were analyzed by the non-parametric Mann-Whitney U-test.
Immunofluorescence analysis. For histochemistry, tissues were mounted on aluminum chucks, frozen in super-cooled n-hexane and stored at-80˚C. Frozen digestive gland sections (10 μm) of ten mussels from each exposure condition were cut by cryostat (Leica CM3050) and flash-dried by transferring them onto poly-L-lysine-coated microscope slides at room temperature. Digestive gland tissue sections were fixed in PFA solution (4% in PBS, pH 7.2) for 20 min at 20 ± 1˚C.
Immunofluorescent anti-tubulin staining: after fixation, sections were washed three times in PBS (5 min) and incubated in a permeabilisation and blocking solution (0.5% Triton X-100, 2% bovine serum albumin-BSA, 0.5% goat serum in PBS) for 1 h at 20 ± 1˚C. After rinsing, sections were incubated with the primary antibody (rabbit polyclonal to tubulin, Abcam, 1/100 in PBS containing 1% BSA and 0.05% Triton X-100) overnight at 4˚C, and then in the secondary antibody i.e. goat anti-rabbit IgGH&L (Chromeo™ 488) (Abcam) (1/100 in 1% BSA and 0.05% Triton X-100 in PBS) for 1 h at 20 ± 1˚C in the dark. Sections were then rinsed in PBS and then mounted in Mowiol mounting medium (Cold Spring HarbProtoc, 2006). Controls for non-specific staining included sections that were processed in the absence of the secondary antibodies. Slides were viewed under 400 x magnification by an inverted photo-microscope (Zeiss Axiovert 100M connected to a digital camera Zeiss AxioCamMRm) equipped for fluorescence microscopy using FITC emission filter [44]. Images were analysed using an image analysis system (Scion Image) that allowed for the quantification of the mean fluorescence intensity. Data were analysed by the non-parametric Mann-Whitney U-test.
F-actin staining: after fixation, sections were washed three times in PBS (5 min) and incubated in a permeabilisation solution (0.5% Triton X-100) for 1 h at 20 ± 1˚C. After rinsing, sections were incubated with Green Fluorescent Phalloidin Conjugate (CytoPainter F-actin Labeling Kit-Green Fluorescence-Abcam) for 1 h at 20 ± 1˚C. Sections were then rinsed in PBS and then mounted in Mowiol mounting medium. Slides were viewed under 400 x magnification by an inverted photo-microscope (Zeiss Axiovert 100M connected to a digital camera Zeiss AxioCamMRm) equipped for fluorescence microscopy using a FITC emission filter.

B[a]P in water and tissue samples
The performance of the chemical analyses was determined by calculating the recovery after extracting spiked samples. Water samples were spiked with B[a]P at 4.3 and 33 μg/L in triplicate and extracted following the protocol described in the Methods section. The method was shown to be accurate and repeatable with recoveries of 99 ± 2% and 102 ± 4%, respectively. Regarding tissue analyses, gills and DG from unexposed mussels were spiked with B[a]P at 10 and 50 μg/g dry weight, respectively. Recoveries of 99 ± 3% and 104 ± 3% for the gills and digestive gland (DG), respectively, indicate accuracy and repeatability, and that the mussel's high lipid content had no matrix effects on the analyses at these selected concentrations.
Measured water concentrations of B[a]P at time 0were in good agreement with nominal values at 63%, 104% and 121% for 5, 50, and 100 μg/L, respectively, The B[a]P concentration in water did not change significantly between 0 and 1 h, but decreased to approximately 10% of their initial values within 24 h (S1 Fig). The B[a]P concentration in both DG and gill tissues was below the limit of detection for control mussels (<0.5 μg g -1 dry weight)and increased with the concentration of exposure (Table 1). In addition, the lipid concentration also increased with increasing exposure concentrations for both tissues. In gills, the water content was higher and the lipid content was lower than in DG.

Microarray architecture
A final set of 15 biological processes at the cellular and tissue levels were retained: lysosomal activity, DNA metabolism, translation, mitochondrial activity, carbohydrate metabolism, heat shock response, oxidative stress, proteolysis, phase I and II xenobiotic metabolism, metal resistance, endocrine disruption, immune response, pumps, cytoskeleton, and adhesion to substrate. A minimum of 10 and a maximum of 80 target genes were considered for each biological process depending on the available annotated sequences (S2 Table). All probes were designed starting from complete or partial cDNA sequences. Oligo probes (60-mer) were designed for each target gene, covering the 3' region of the mRNA sequence and with particular attention on guaranteeing optimal gene specificity in the case of homologous genes, such as metallothionein isotypes. A total of 465 probes (S2 Table) that were triplicated randomly to fill a 1395 feature sub-microarray for a total of eight sub-arrays per slide (

Microarray test and applications
Dual colour hybridization revealed evident alterations in gene expression patterns in the investigated tissues due to the B[a]P treatments ( Table 2, Fig 1 and S3 Table). Our data revealed distinct patterns for 205 and 109 differentially expressed genes (DEGs) in the gills and DG, respectively. Of the 205 DEGs in gills, only 14 genes were shared among the three experimental condition datasets, whereas no target was shared between the three conditions in DG tissues (Fig 2 and S3 Table).As our new platform is based on annotated genes that are already ranked into well-known biological processes (S2 Table), thus it was easier to generate the list of New insight into Benzo[a]pyrene effects on marine mussels   biological processes involved in the mussel's response to B[a]P-induced stress from the DEGs obtained for each condition (Figs 3 and 4). Our data highlighted 13 contributing biological processes in gills in response to B[a]P concentrations. The response was largely dominated by the down-regulation of translation, the immune response, mitochondrial activities, and cytoskeleton-related genes in both 50 and 100 μg/L exposures (Fig 3, S3 Table). In mussels exposed to 5 μg/L B[a]P, only seven biological processes, represented by a maximum of five targets (translation), were depicted. Among the 15 biological processes used in the STREM, only "metal resistance" was not involved in the response to B[a]P challenge at the three investigated concentrations.
The main contributing biological processes in the mussel's response to B[a]P in DG are reported in Fig 4 (S4 Table). Our data indicate that cytoskeleton, translation, mitochondrial activities, carbohydrate metabolism, immune response, and DNA metabolism (over 14 biological processes) are the most dominating processes in response to increasing concentrations of  Table and S4 Table. Venn diagram representation of gene expression patterns (Panel B) clearly depicted that no DEGs are shared between the three experimental conditions. All DEGs are obtained respect to the control condition. Data used to generate the Venn-diagram were obtained from microarray analysis (S4 Table). Four biological replications of were considered. https://doi.org/10.1371/journal.pone.0178460.g002 New insight into Benzo[a]pyrene effects on marine mussels B[a]P. Furthermore, the lowest B[a]P concentration rendered only six biological processes, largely characterized by the up-regulation of the immune response, oxidative stress, and cytoskeleton-related genes and to lesser extend DNA metabolism. Over-representation analysis of DEGs in the gills tissue of mussels exposed to B[a]P (5μg/L, 50μg/L and 100 μg/L). Showed are: experimental condition; biological processes; Number of up-regulated genes; Number of downregulated genes. The over-represented biological processes in B[a]P-exposed animals versus control. (S3 Table). https://doi.org/10.1371/journal.pone.0178460.g003 New insight into Benzo[a]pyrene effects on marine mussels Over-representation analysis of DEGs in the digestive gland tissue of mussels exposed to B[a]P (5μg/L, 50μg/L and 100 μg/L). Showed are: experimental condition; biological processes; Number of up-regulated genes; Number of down-regulated genes. The over-represented biological processes in B[a]P-exposed animals versus control. (S4 Table). https://doi.org/10.1371/journal.pone.0178460.g004 New insight into Benzo[a]pyrene effects on marine mussels qRT-PCR We carried out qRT-PCR to confirm and refine the relative expression levels of 10 homologous genes belonging to the most important biological processes, including the genes encoding caspase, p53, DNA ligase and topoisomerase (DNA metabolism), tubulin (cytoskeleton), hsp70, hsp27, (heat shock response), catalase, Zn-Cu superoxide-dismutase (oxidative stress), and the gst(phase II detoxification). Gill tissue was used for qRT-PCR analysis because it exhibited the most important responses to B[a]P and stronger bioaccumulated B[a]P than DG tissue (Table 1) Biomarkers of genotoxicity DNA damage using the comet assay. The alkaline version of the comet assay was used to detect DNA damage including single-and double-strand breaks and alkali-labile (e.g. apurinic) sites. Induction of DNA damage was greater in DG than gills. In the DG, significant effects occurred at both 50 and 100 μg/L, compared to only 100 μg/L for the gills (one-way ANOVA, p<0.05) (Fig 5).
Bulky DNA adduct formation by 32 P-postlabelling. There was a distinct tissue-specific pattern in the occurrence of DNA adducts; all but one gill tissue sample had none detectable. This finding is in sharp contrast to DG tissues, 11 samples of which showed the presence of B [a]P-derived DNA adducts. On the TLC sheets the adduct pattern consisted of a single adduct spot (spot 1) (Fig 6), previously identified by mass spectrometry as dG-N 2 -BPDE [40]. Analysis of DNA adduct data using glm revealed significant effects of B[a]P concentration (p<0.0001) and tissue type (p<0.0001) on the level of adducts, as well as significant interactions between B [a]P concentration and tissue (p<0.05). The two higher B[a]P concentrations comprised all of the adduct-containing samples (Fig 6A).Sex-specific differences between these two concentrations (100 and 1000 μg L -1 B[a]P) were evident, in that female mussels showed significantly higher levels of adducts at both concentrations, although n numbers were < 5 for two groups (Fig 6B).

Effects at tissue level
Lipofuscin accumulation. As shown in Fig 7, the amount of the lipofuscin accumulation in the digestive gland lysosomes of B[a]P-exposed mussels increased significantly with respect to controls, reaching a plateau already in the animals exposed to the lowest chemical concentration (5 μg/L).
Immunofluorescence analysis of tubulin and F-actin in digestive gland tissue sections. The consequences on cell physiology of the changes observed at transcriptomic level have been also investigated by looking at the cytoskeletal organization.
In particular, the possible effects on tubulin, which mRNA content increases in B[a]Pexposed mussels, and actin showing no change of the transcription of its coding gene, were evaluated using an immunohistochemical approach. The results reported in Fig 8 clearly show that the tubulin amount in the digestive gland cells increases in mussels exposed to B[a]P increasing concentrations with a bell-shape trend reaching a maximum in the animals exposed to 5-50 μg/L. In these animals, the cytoplasm of the cells contains numerous highly fluorescent granules that contribute to determine the observed tubulin increase.
In the digestive gland cells of the mussels exposed to B[a]P, the amount of F-actin seems to be only slightly decreased but the cytoskeletal architecture is, at least in part, lost. At the higher B[a]P concentrations, actin cytoskeletal alterations are particularly evident for what concerns the cortical cell compartment and its role in the structure of the digestive gland tubule organization (Fig 9).

Discussion
We have presented the development of a new low density oligonucleotide microarray to study alterations in mussel gene expression profiles after exposure to anubiquitous environmental stressor. Genes represented on the microarray were identified from public databases, mainly NCBI Genbank. The gene list included 465 annotated genes known to be implicated in the stress response of bivalves [3,11,4,29,2]. The annotated genes were collected in different gene data banks and represent a more complete list among genes involved in the selected pathways involved in the mussel's stress response. Some features of the novel 'STREM-chip' that guarantee a good performance and reliable data include the simultaneous presence of eight replicated sub-arrays on the same slide and the use of a single probe covering the 3' region of the mRNA, Data are shown with the sexes combined (a) and separated (for those concentrations that resulted in detectable numbers of adducts)(b). All treatments had an n of 9 in total but the distribution of sexes was not even, therefore group n values are indicated in (b). Mismatching lowercase letters indicate differences between groups according to Tukeys HSD applied as post hoc tests following a generalised linear model (p< 0.05; n = 9). Insert: Representative autoradiographic profile of DNA adducts, measured by thin-layer chromatography 32 P-postlabelling, in digestive gland tissue exposed to B[a]P; no B[a]P-derived DNA adducts were detected in control (untreated) tissue (data not shown). The origin (OR), at the bottom left-hand corners, was cut off before exposure. The arrow indicates 10-(deoxyguanosin-N 2 -yl) -7,8,9-trihydroxy-7,8,9,10- New insight into Benzo[a]pyrene effects on marine mussels improving the accuracy of measurements. Moreover, to optimize data confidence, dual color hybridizations were performed in at least quadruplicate. For each sub-array, each target was triplicated and four biological replicates considered, resulting in 12 values for each target. In addition, the eight sub-arrays may assure at least two experimental conditions (considering four biological replicates per condition) in the same slide.
To identify DEGs, we did not set a default threshold for the fold change a priori, but instead the statistical criterion for significance was adopted, in order to appreciate moderate gene modulations. The new platform presents an advantage of being open to any additional target based on the increase in bivalve species sequencing, especially Mytilus sp.
To assess the potential use of the new array platform, we tested the effect of B[a]P, a relevant environmental contaminant, in two mussel tissues, the gills and DG. Chemical data suggested that the concentration of B[a]P was higher in the gills than the DG. The B[a]P concentration in the DG after exposure to 50 μg/L B[a]P was of a similar order of magnitude as in [24] after exposure to 19 μg/L/animal (10.29 μg/g). Such levels are also consistent with the whole body burden of total polycyclic aromatic hydrocarbons (PAHs) in field sampled mussels from the Irish Sea (0.13-22.48 μg/g) [44], suggesting that they are environmentally relevant.
Interesting information could be gleaned from the transcriptomic data. As expected, we observed a marked regulation of genes involved in DNA metabolism in both tissues at B[a]P concentrations of 50 and 100 μg/L. Our data suggest the occurrence of a DNA repair or apoptosis event with respect to the applied B[a]P concentration. Investigation of DNA damage New insight into Benzo[a]pyrene effects on marine mussels response genes coding for caspase (HQ424451.1) revealed a marked up-regulation of apoptotic genes in mussels exposed to 50 or 100 μg/L B[a]P. However, DNA repair-related enzymes (DNA ligase:AJ624686.1; p53:KC545827.1) were markedly up-regulated in animals exposed to 5 μg/L B[a]P (data confirmed by qRT-PCR), except for topoisomerase, which was down-regulated (AF227976.1).
Mussel's cells may respond to DNA damage via p53-mediated cell cycle arrest or programmed cell death [45]. The p53 family of proteins has been demonstrated to be involved in molluscan leukemic hemocytes, and to be inactivated by sequestration in the cytoplasm [46,47]. Furthermore, bulky B[a]P-DNA adducts (i.e. dG-N 2 -BPDE) are reported to be New insight into Benzo[a]pyrene effects on marine mussels removed by nucleotide excision repair (NER) enzymes [20], including DNA ligase [48]. NER is dependent on the base pair conformation, suggesting a direct implication of DNA ligaserelated proteins in such processes. Finally, topoisomerase, a DNA replication/repair-related enzyme has been reported to be highly active in mammalian cells exposed to PAHs [49]. PAH-DNA adducts alter the execution or fidelity of enzymatic DNA transcription and may cause mutations and cancer [50,20]. The dissimilar data obtained in this work for topoisomerase may be explained by the higher replication rate of mammalian cells maintained in culture conditions [50] compared to in vivo-sampled cells, as well as the difference in concentrations used in our study compared to the exposure conditions for model mammalian cells.
Investigations of DNA damage (comet assay) and DNA adduct formation ( 32 P-postlabelling) showed higher genotoxicity in mussels exposed to 50 and 100 μg/L B[a]P, which is in line with the gene expression data. Moreover, our results suggest a higher genotoxicity in DG cells compared to gill cells. Previous work on the location and concentration of mono-oxygenase enzymes required for the transformation of B[a]P into its reactive intermediates in mussels, indicates much higher concentrations in DG [51], which supports our finding of higher % tail DNA values in this tissue compared to the gills, despite less B[a]P accumulation. Oxidative DNA lesions (i.e., the formation of 8-oxo-dG) have also been found to be more prevalent in DG cells than the gills after B[a]P exposure [21,24]. Livingstone et al. [51] reported cytochrome P450 enzyme activity in the DG only and the highest activities of P450-related enzymes (including B[a]P hydoxylase) in this organ (i.e.DG). This is consistent with higher levels of B [a]P-derived DNA adducts (i.e. dG-N 2 -BPDE) in this tissue / organ reported here. The cytochrome P450 monooxygenase system also provides a potential explanation for the sex-specific difference in BPDE-adducts that we observed, as the activities of these enzymes were also New insight into Benzo[a]pyrene effects on marine mussels higher in females than males, implying that more of the genotoxic intermediates may be present in their digestive gland [51].
Another marked response of mussels to increasing concentrations of B[a]P in gills is the down-regulation of genes related to the translation process. Thirty-five of 38 DEGs and 31 of 36 DEGs were down-regulated in mussels exposed to 50 and 100 μg/L B[a]P, respectively, though the interpretation of this may be complex because these genes also encode proteins that are implicated in multiple cellular pathways. Of these down-regulated DEGs, 30 and 38 encode ribosomal protein subunits (S4 Table) in mussels exposed to 50 and 100 μg/L B[a]P, respectively. This down-regulation may indicate that mRNA-directed protein synthesis is reduced in mussels exposed to higher B[a]P loads. PAHs are known to interfere with protein synthesis, as the abundance of proteins involved in protein synthesis and degradation, ATP supply and structural proteins changes in response to B[a]P in M. Galloprovincialis [26]. Ribosome biogenesis has been reported to be affected by environmental stressors in bacteria [52]. In this context, it has also been suggested that both transcription and translation are regulated in zebrafish in response to environmental contaminants [53].
In this study, mitochondrial activities and carbohydrate metabolism are reported to be among the processes contributing to the mussel's response to B[a]P exposure in two important target tissues. Recently, it has been reported that the regulation of ATP synthase F0 subunit 6 in mussels exposed to 10 μg/L B[a]P using a proteomics approach [26]. The regulation of ATP concentration is linked to inner mitochondrial membrane, through which the respiration process supplies energy to cells. Similarly, proteomic analysis of freshwater mussel D. polymorpha exposed to B[a]P for 7 days identified proteins associated with cellular metabolism (N-acetyltransferase 8-like and probable aspartate aminotransferase) [54].
B[a]P exposure can be considered as a potent inducer of oxidative stress in mussel tissues [25,55,26,56]. The production of oxy-radicals in mussels has been reported to be mediated by a higher range of contaminants, including B[a]P [57,58]. In this work, even though the number of DEGs involved in the oxidative stress response present in the new platform (i.e., 10) was less than other marked processes (e.g. translation or immune response), the results indicate clear effects as 80% of the genes involved in the oxidative stress response were affected in the gills and DG of mussels exposed to higher B[a]P concentrations. In this case, the observed transcriptomic changes were indicative of oxidative damage in the cells, and the lipofuscin lysosome content was significantly increased. In mussels exposed to higher B[a]P concentrations, the heat shock response was also depicted as a contributing process in the gill and DG tissues. Among the stress proteins, the heat shock proteins are molecular chaperones involved in the assembly, folding, and intracellular transport of proteins, protecting against stress-associated cellular damage [59].
Among the genes present in our microarray that code for protein folding are calreticulin and FK506-binding protein. Calreticulin is known to bind mis-folded proteins, preventing their export from the endoplasmic reticulum to the Golgi apparatus [60]. Small heat shock proteins, such as HSP27, HSP26, and FK506-binding protein, are likely to act as chaperones for the maintenance of cytoskeletal structural elements during stress [61]. Calreticulin and FK506-binding protein were recently reported to be highly up-regulated in M. galloprovincialis with respect to M. trossulus exposed to higher heat stress. In addition, oxidative stress may induce heat shock protein phosphorylation, resulting in actin binding and the induction of apoptosis [62]. Thus, chaperone induction may represent an effort to protect the cells against the stress effect caused by impaired cellular responses associated with transcriptional regulation.
Exposure to PAHs has a number of effects on the immune function in marine organisms. Decreased lymphocyte proliferation has been observed in Japanese medaka (Oryziaslatipes) [63] and bluegill sunfish (Lepomismacrochirus) following dietary exposure to three different PAHs [64]. Altered in vitro proliferation of leukocytes from common carp has been reported after exposure to 3-methylcholanthrene. In our case, the immune response process was affected at the transcriptional level in response to B[a]P exposure in the two investigated tissues. This may attempt to speculate a possible effect ofB[a]P also in cells involved in immune response which are always present in the investigated tissues. The down-regulation trend was most dominant in the DEGs involved in immune response, suggesting a possible immune-toxicity of the applied concentrations. The study of the B[a]P effect on the immune cells could contribute to clarify this aspect.
The transcriptional response of mussels exposed to increasing concentrations of B[a]P was also marked by the down-regulation of genes encoding proteins involved in cytoskeleton maintenance and repair (Figs 3 and 4). Cytoskeletal defence was previously proposed as a potential mechanism for increased stressors, such as thermo-tolerance in mussels [6]. Increased production of reactive oxygen species (ROS) has been reported to severely affect the actin cytoskeleton in mussel hemocytes following exposure to B[a]P [56]. Recently, it has also been induction of α-tubulin in B[a]P-exposed mussels [27]. This is a key cytoskeletal component responsible for microtubule polymerization, cell transport, and motility. The study therefore reinforces the influence of all PAHs in cytoskeleton organization [27].
In the present study, tubulin appears to be one of the few cytoskeleton genes up-regulated following exposure to different concentrations of B[a]P. To complement B[a]P induced transcriptomic changes, the effects of B[a]P on the cytoskeleton structure of the digestive gland cells was also investigated. The increase in the tubulin mRNA content showed an enhancement of this protein in the cytoplasm of the digestive gland cells. Using an immuno-histochemical approach, it was possible to verify that in B[a]P exposed mussels, it induces only limited effects on the tubulin filament structure but a large amount of tubulin was found to be sequestered in cytoplasmic granules. The formation of granules containing tubulin has been previously described in mammalian and molluscan cells due to different stress/pathological conditions [65,43]. It could therefore hypothesized that B[a]P affect tubulin filaments via different mechanisms(e.g. by inducing oxidative stress) with the sequestration of the damaged tubulin monomers and part of the filaments in the granules. The decrease of free tubulin monomers should stimulate the tubulin neosynthesis [66] through the enhancement of the protein gene transcription in an attempt to maintain in the cells the cytoskeleton architecture needed for important cell functions such as vesicles transport and autophagic activity.
We did not observe any change on the actin mRNA levels in B[a]P exposed mussels. This indicates the absence of effects of the chemical on the actin cytoskeleton. In contrast, the immuno-fluorescence analysis of F-actin clearly demonstrated a relevant alteration of the actin structure in the digestive glands of treated mussels. The noxious effects of B[a]P on the F-actin architecture have been previously reported by various authors in mammalian cells [67] and in mussel's haemocytes [56]. In this last case, the B[a]P effects were mainly related to the cortical component of the actin cytoskeleton structures. Our results demonstrated that B[a]P affects the actin ring in the cortex area involved in cell-cell contacts. These effects, which are particularly relevant at the higher B[a]P concentrations (50-100 μg/L), also imply an alteration of the organization of the digestive gland tubules.
It should be noted that the exposure to H 2 O 2 of mammalian monolayer intestinal cells showed similar effects on the actin cytoskeleton with the loss of the actin structure present in the cortex area [68]. In our study, mussel exposed to B[a]P induced oxidative stress as demonstrated by the significant increase of lysosomal lipofuscin accumulation and by the activation of the transcription of genes related to the antioxidant responses (e.g.catalase and SOD). It therefore seems possible that oxidative stress may represent an important component in the mechanism of action of B[a]Pin mussels, at least for what concerns the observed negative effects on the cytoskeletal structure of the digestive gland cells.
In the present study, adhesion to substrate was modulated at the transcriptomic level in both tissues upon exposure to B[a]P. Similar results were observed from proteomic characterization of the DG from M. galloprovincialis after exposure to sub-lethalB[a]P concentrations. The down-regulation of this process by organic compounds could lead to impairment of the force-maintaining structures during the catch state, compromising the adhesion (byssus strength) and motility processes in exposed mussels. Investigation of this parameter would provide a link between transcriptomic studies and organism responses which could be extrapolated to population dynamics studies.
Finally, it is important to note that transcriptomic data may provide an overview of the cell's response to stress in term of changes in the genes products. Changes in the mRNA profiles may be related to the maintenance of cell/tissue homeostasis. In this case, the organism appears to be healthy but has a decreased capacity to adapt to further environmental changes. However, alteration of the mRNA profile could also indicate that the organisms are not more capable of coping with the toxic environment [69]. To clarify these different situations, it is important to associate the transcriptomic data with the results of functional studies of cells/tissues related to different stress response pathways. One advantage of the new platform is that each biological process represented in the array could be related to the alteration of a particular biological function using previously developed biomarkers. In this work, gene expression changes were related to cytoskeleton organization (tubulin immuno-histochemistry), oxidative stress (lipofuscin content), and innovative studies on DNA damage and autophagy regulation (Sforzini et al., unpublished data).

Conclusions
In conclusion, we have demonstrated that the novel low density DNA microarray is a robust and valid tool for obtaining insights of the molecular responses of Mytilussp., a model organism widely used in ecotoxicological studies. The new platform allows easy retrieval of the involved biological processes in response to stress because targets are already clustered according to their biological implications. The new platform remains open to eventual implementation of any new targets or processes in view of the full genome sequencing of Mytilus sp. Moreover, the presence of eight sub-arrays per platform is associated with a reduction in the cost/experimental effort and will allow future routine use of this microarray in field or biomonitoring studies. The integration of transcriptomic data with the effects on cell/tissue functions as determined by other novel sub-lethal biomarker responses provide a more biologically meaningful interpretation of the physiological alteration, e.g. rearrangement of nucleus activity on the capability of the mussels to adapt to toxic environments. In particular, despite an increase in activity of the genes coding for proteins involved in DNA repair, DNA damage was observed. The changes in DEGs coding for cytoskeleton proteins reflect an alteration in its structure in DG cells and the increase in mRNA coding for antioxidant enzymes is not sufficient to protect the cells from oxidative stress.

S1 Fig. Concentration of benzo[a]pyrene (B[a]P) in seawater over the first 24 h of 3 d expo-
sure at 0-100 μg L -1 . Data are means ± SD (n = 3). Dashed lines represent the nominal concentrations at 5, 50 and 100 μg L -1 . At 24h, B[a]P concentrations in the 5 μg L -1 group were below the limit of detection (0.25 μg L-1). Although a 1000 μg L -1 exposure was also performed (for DNA adduct analysis only), water concentrations were not measured for this treatment for logistical reasons. (DOCX)

S2 Fig. A typical dual color hybridization analysis of Cy3/Cy5-labelled cDNAs from B[a]Ptreated vs control mussels, obtained by means of a dual laser source microchip scanner.
For data normalization, the two channels were balanced on RNA intensities. (DOCX) S3 Fig. qRT-PCR confirmation of microarray data. Targets expressions have been analyzed by real-time PCR, using a 18S rRNA, Beta actin and Ribol27 as reference genes for data normalization. Microarray data for the investigated genes were confirmed. These results indicate a high accuracy and sensitivity of the STREM-Ship that was able to detect even small change of expression. Data represent the mean of at least four independent experiments. Calculation of relative expression levels and statistics (pairwise randomization test, p < 0.05) were obtained using the REST software [37]. Experimental coefficient of variation (CV) was below 5% for all the investigated targets. (DOCX) S1 Table. qRT-PCR primers and Taqman probes.