Transcriptomic response and hydrocarbon accumulation in the eastern oyster (Crassostrea virginica) exposed to crude oil

Comparative Biochemistry and Physiology, Part C 225 (2019) 108571

Contents lists available at ScienceDirect

Comparative Biochemistry and Physiology, Part C journal homepage: www.elsevier.com/locate/cbpc

Transcriptomic response and hydrocarbon accumulation in the eastern oyster (Crassostrea virginica) exposed to crude oil

T

Edgar A. López-Landaverya, Gerardo Amador-Canob, Naholi Alejandria, Nancy Ramirez-Álvarezc, ⁎ Isidro Montelongob, Fernando Díaza, Clara E. Galindo-Sáncheza, a

Department of Marine Biotechnology, Centro de Investigación Científica y Educación Superior de Ensenada (CICESE), Ensenada, BC, Mexico Universidad Tecnológica del Mar de Tamaulipas (UTMART), La Pesca, Soto La Marina, Tamaulipas, Mexico c Instituto de Investigaciones Oceanológicas (IIO), Universidad Autónoma de Baja California (UABC), Ensenada, BC, Mexico b

A R T I C LE I N FO

A B S T R A C T

Keywords: Gene expression RNA-Seq Transcriptome Xenobiotic

The adverse effect of crude oil on marine invertebrates is well known. To have a better understanding of its effects on marine invertebrates, Crassostrea virginica was exposed to different concentrations (50, 100 and 200 μg/L) of a mixture of super-light and light crude oil for two weeks, evaluating the transcriptomic response of the digestive gland using RNA-Seq and their accumulation in soft tissues. A total of 33,469,374 reads were assembled, which resulted in 61,356 genome assemblies (‘Genes’). Trinotate was used for transcript annotation. At the end of this process, 86,409 transcripts were maintained, comprising a broad set of enzymes from xenobiotics metabolism, oxidative stress, stress and immune responses, and energetic metabolism. The enrichment analysis revealed a change in biological processes and molecular functions, finding from 100 to 200 μg/L. Moreover, the differential gene expression analysis showed a dose-dependent transcriptional response, generally up to 100 μg/L and in some cases up to 200 μg/L, which suggested that oysters' response decreased after 100 μg/ L; the analysis of crude oil presence in soft tissues indicated that C. virginica is a suitable candidate for ecotoxicology. Finally, these results should contribute to expanding current genomic resources for C. virginica. Furthermore, they will help to develop new studies in aquatic toxicology focused on knowledge in depth of metabolic pathways, jointly with other approaches (such as proteomics) to allow obtaining a complete idea about the eastern oyster response to crude oil.

1. Introduction In recent years, the interaction between crude oil and the marine environment has been studied widely due to the negative impact on its biological component (Enwere, 2009). The adverse effect of an oil spill on marine life varies with location, spill magnitude, life stage, habitat, sensitivity and the capacity of organisms to prevent or process contaminants (Eisler, 1987; Blackburn et al., 2014). Crude oil is composed of a mixture of aliphatic hydrocarbons (AHs) and monocyclic and polycyclic aromatic hydrocarbons (MAHs and PAHs) that vary in density, chemical structure, molecular weight and interaction with sediment and water (Fingas, 2001). PAHs are among the most potent carcinogenic, teratogenic and mutagenic components; they are toxic to the immune system in a diversity of aquatic species and humans (Eisler, 1987; Reynaud and Deschaux, 2006; Omar-Ali et al., 2015), and they have also shown a high level of bioaccumulation (Schäfer and Köhler, 2009).

⁎

One way to assess the environmental impact of an oil spill is through the study of exposure effects of crude oil on organisms that live throughout the water column and on the bottom, as indicators of the health status of coastal ecosystems (Neff et al., 2000). However, different responses have been reported when vertebrates and invertebrates are compared. Bioaccumulation in marine invertebrates depends on the type of crude oil released, as well as the organism's feeding strategy and ability to metabolize and excrete oil components (Livingstone, 1998; Neff et al., 1985; Meador, 2003). For example, filter-feeding oysters and mussels may incorporate hydrocarbons at a higher rate than invertebrate grazers or predators and have higher levels of bioaccumulation, especially compared to marine vertebrates (Lech and Bendt, 1980; Widdows et al., 1982; Neff, 2002; Martínez-Gómez et al., 2010). Organisms have elaborate systems of biotransformation with enzymes that can metabolize organic compounds, such as PAHs. The metabolites generated by biotransformation are less complex, easier to excrete and usually associated with endogenous substrates by enzymes,

Corresponding author at: Department of Marine Biotechnology, CICESE, Mexico. E-mail address: [email protected] (C.E. Galindo-Sánchez).

https://doi.org/10.1016/j.cbpc.2019.108571 Received 28 February 2019; Received in revised form 5 July 2019; Accepted 9 July 2019 Available online 12 July 2019 1532-0456/ © 2019 Elsevier Inc. All rights reserved.

Comparative Biochemistry and Physiology, Part C 225 (2019) 108571

E.A. López-Landavery, et al.

oysters, one pool per experimental unit). Treatments corresponded to control (without crude oil) and 50, 100 and 200 μg/L, respectively.

such as glutathione S-transferases (Hu et al., 2016). Studies have shown that both adaptive and innate immune systems of vertebrates respond to crude oil exposure with apparent mechanisms of action that include P450-mediated biotransformation, aryl hydrocarbon receptor binding, calcium mobilization, and heat-shock protein 90 (Reynaud and Deschaux, 2006; Roméo and Wirgin, 2011). Mollusks exposed to pollutants, such as crude oil, have shown low growth and reduced reproductive success (Sindermann, 2006). That condition has generated impairment of the immune system making the organisms more vulnerable to stress, diseases, and predators (Grundy et al., 1996) and generally reducing their fitness, generating adverse effects on the marine environment and leading to economic losses (Vignet et al., 2016). For example, oyster reefs serve as habitat for other aquatic invertebrates and fish (Karnauskas et al., 2013) associated with nitrification and denitrification processes and deposition of phosphorous (Newell et al., 2005; Caffrey et al., 2016). Moreover, oyster fishery production in Mexico ranged from 42,945 to 50,715 t/year with the states of Veracruz, Tabasco, and Tamaulipas providing 87.5% of the total volume in average from 2010 to 2013 (SAGARPA, 2013). The eastern oyster (Crassostrea virginica Gmelin, 1791) is a sessile organism that can be found in coastal and estuarine habitats from the Gulf of Saint Lawrence (Canada) to subtropical Mexico, spanning a wide range of physical and biotic conditions (Varney et al., 2009). It is an economically and ecologically valuable resource (Anderson et al., 2015), which has shown to be, along with other organisms, a suitable species to assess the level of environmental contamination (Sericano et al., 1996; Aguilar et al., 2012). In Mexico, increase in crude oil exploration, production and transport represents a risk of accidental oil spills with potentially adverse effects on the marine ecosystem. For example, an increase in PAH levels from 80,000 and 230,000 t/year has been reported in estuarine and marine environments as a result of anthropogenic activities (Kennish, 1991; Wrigth and Welburn, 2002). Furthermore, the Deepwater Horizon oil spill released approximately three million barrels of crude oil into the Gulf of Mexico (Finch et al., 2016). Based on the above, the objective of this study was to assess the effects of crude oil exposure under experimental conditions on the transcriptomic response of the eastern oyster C. virginica.

2.3. Oyster acclimatization At the UTMART laboratory, oysters were distributed into two rectangular 250-L tanks (260 and 265 oysters/tank) and kept with constant aeration. Water exchanges were carried out daily, and salinity was gradually increased from 19 to 33 g/L (3 g/L/day) diluting seawater to simulate environmental salinity conditions of most of the year (Sistema de Venecia, 1958; Shumway, 1996; Varinka Sáenz, UTMART, pers. comm.). When salinity reached 33 g/L, animals were maintained at this salinity with water temperature from 22 ± 1 °C for two weeks.

2.4. Crude oil exposure, experimental design and feeding Two types of crude oil were provided by Petróleos Mexicanos (PEMEX): super-light oil (40 American Petroleum Institute gravity or 40° API, 0.82 g/cm3) and light oil (35° API, 0.85 g/cm3). The experiment consisted of a completely randomized design, which included four treatments of a crude oil mixture (control, 50, 100 and 200 μg/L) and three replicates (n = 40), generating twelve experimental units. The crude oil mixture was obtained by combining super-light and light oils in a 1: 1 ratio (w:w). Then, it was homogenized in acetone by pipetting and shaking 200 rpm at dark for 24 h in a 2:1 ratio (v:v, crude oil: acetone), as reported in Jeong and Cho (2007). Crude oil mixture was added into the experimental units, and concentrations were adjusted based on the volume of the experimental units. Only acetone was used in the control treatment. Water (from experimental units) was exchanged every two days (50% of the total volume), and crude oil mixture was added after water exchange and proportional to the volume removed to maintain experimental concentrations. Exposure time to the crude oil mixture was 14 days. The oysters were fed with a concentrated solution of microalgal mixture (Shellfish Diet 1800™, 4–20 μm, ~2.0 × 109 cells mL−1, Reed Mariculture, Campbell, CA, USA). Daily ration was calculated based on oysters' dry weight (6% per day), and they were fed at 08:00 h and 16:00 h.

2. Materials and methods 2.5. Detection of crude oil compounds in oysters, water and on sediment 2.1. Oyster samples 2.5.1. Extraction of aliphatic hydrocarbons and polycyclic aromatic hydrocarbons in oysters and on sediment The extracts for determining AHs and PAHs in whole soft tissues and on sediment were obtained based on Murphy et al. (2012) with modifications. Briefly, A Dionex ASE 350 solvent-accelerated extraction system (Thermo Scientific, USA) equipped with 34-mL stainless steel cells was used. Each cell was packaged sequentially with a filter GF/F, 5 g of florisil and 6 g of basic alumina, both fully activated. In the case of sediment, silica was used instead of florisil. Next, a white sand layer was added before adding mixture homogenization resulting from 1.5 g of the dry tissue sample and 3 g of diatomaceous earth. Successively, the remaining space in the cell was filled with white sand (before use, resins and sand were cleaned at 400 °C for 4 h). Then, surrogate standards were added. The extraction method consisted of four static 5min cycles at 100 °C at a constant pressure of 1500 psi with a replacement volume of 60% of the cell capacity and using dichloromethane as extraction solvent. After that, the extracts obtained were evaporated by a Rocket evaporator, transferred to concentration tubes and evaporated under Nitrogen flow until almost dry. The extracts were carried to a volume of 1 mL (oysters) and 0.5 mL (sediment) with hexane. Cleaning of extracts was done by liquid chromatography by column. Finally, the extracts were resuspended up to a final volume of 0.5 mL with isooctane, to which the internal standards were added before their analysis in gas chromatography.

Oysters were obtained from Laguna de Morales located at La Pesca, Tamaulipas, Mexico (Fig. 1). This lagoon was granted food safety certification in 2014 (Ivonne Padrón, Fisheries and Aquaculture Office, Tamaulipas, pers. comm.). Oysters (n = 525) were collected from three sampling points (175 oysters per point) in a shallow subtidal zone. Organisms were kept in a rectangular tank with brackish water to a salinity of 19 g/L, temperature of 22 ± 1 °C and constant aeration. They were transferred to the laboratory of shellfish production at Universidad Tecnológica del Mar de Tamaulipas (UTMART). Oysters' shell length was from 49 to 60 mm. 2.2. Presence of aliphatic hydrocarbon and polycyclic aromatic hydrocarbon in oysters, water and on sediment To determine the presence of AHs and PAHs from sources different from crude oil, whole soft tissue (three pools of 20 oysters), sediment (2 × 100 g) and water samples (2 × 1000 mL) were collected during field sampling at the lagoon. To determine bioaccumulation of crude oil components during the experiment, whole soft tissue samples were collected from two-week acclimatized oysters (three pools of 15 oysters, one pool per experimental unit), representing the control at day 0. Then, at 7 and 14 days after crude oil exposure, whole soft tissue samples from 45 oysters per treatment were collected (three pools of 15 2

Comparative Biochemistry and Physiology, Part C 225 (2019) 108571

E.A. López-Landavery, et al.

Fig. 1. Location of the sampling area and points of the eastern oyster Crassostrea virginica in Laguna de Morales, Tamaulipas, Mexico. The sampling points are located in the area that was certified in terms of sanitary regulations (Ivonne Padrón, certifier, pers. comm.)

2.5.3. Aliphatic hydrocarbons and polycyclic aromatic hydrocarbon detection in oysters, water and on sediment For AHs and PAHs determination, the extracts obtained were analyzed by gas chromatography (Hewlett-Packard, model 6890, Wilmington, DE, USA) coupled to a flame ionization detector (GC-FID) and a mass selective detector (GC-MSD) in electronic impact mode (EI), respectively. For AHs, during the separation of the components of interest, a column DB-5MS (Agilent Technologies Inc., Wilmington, DE, USA) of 30 m × 0.32 μm × 0.25 μm was used with a constant flow of 1.6 mL/min and N2 as a gas hauler. For PAHs, the column dimensions were 30 m × 250 μm × 0.25 μm with a constant flow of 1 mL/min and Helium (He) as a gas hauler. The program in the oven for AHs consisted of an initial temperature of 60 °C for 3 min and a ramp of 8 °C/min until a final temperature of 300 °C was reached. The initial temperature for PAH was 70 °C for 0 min, a ramp of 5 °C/min and a final temperature of 325 °C. The injection mode was splitless, and quantification was performed through the standard internal method.

2.5.2. Extraction of aliphatic hydrocarbons and polycyclic aromatic hydrocarbons in water The extracts for determining AHs and PAHs in water were obtained based on guidelines from Glaser et al. (1981). AHs and PAHs extraction was standardized based on a volume of 1800 mL of the collected water sample. Before extraction, 50 μL of the standard was added to each sample for monitoring aliphatic compounds. Next, the process consisted of three extractions with dichloromethane as follows: for the first extraction, 60 mL of dichloromethane were added and stirred manually for degasification. Then, it was placed in stirring for 5 min, allowing it to stand vertically to separate the phases in a universal plate stand and extract dichloromethane in a collector flask. The same procedure was performed for the second and third extraction. Next, sodium sulfate anhydrous was added to catch water remnants, and the samples were concentrated in glass containers from the Rocket Synergy evaporator (Thermo Scientific, USA) to 2 mL, almost to dryness. The extract obtained was transferred to concentration tubes by several washes of the receiver container with dichloromethane, and the tubes were concentrated by a Nitrogen grade HP (high purity) flow to a volume of 0.5 mL and finally resuspended to 1 mL with hexane. 3

Comparative Biochemistry and Physiology, Part C 225 (2019) 108571

E.A. López-Landavery, et al.

performed by applying the Benjamini-Hochberg method (Benjamini and Hochberg, 1995). Transcripts represented with an adjusted P-value (FDR) < 1E-02 and at least a two-fold change (FC ≥ 2) were considered significantly differentially expressed in the pairwise comparison of the samples.

2.6. Tissue sampling, RNA extraction and DNase digestion Samples from the digestive gland were dissected out from oysters at times T0 (0 days, 0 μg/L) and T2 (14 days, 50–200 μg/L). Four organisms were selected per experimental unit, generating a n = 12 per treatment. All samples were fixed in RNALater (Ambion, Austin, TX, USA) and stored at −80 °C until RNA purification. Total RNA was extracted with the RiboPure™ Kit (Ambion, Austin, TX, USA) as described by the manufacturer. Briefly, 100 mg of sample were homogenized in 1 mL of Tri Reagent® (Sigma, St. Louis, MO, USA), centrifuged and supernatant precipitated at 4 °C with chloroform. Next, the aqueous phase with absolute ethanol was passed through a filter cartridge, which was washed twice with wash solution. Total RNA was eluted with 50 μL of elution buffer. RNA quantity was verified through absorbance (NanoDrop ND-2000, Thermo Scientific, USA), and RNA integrity was confirmed with a denatured 1.2% agarose-formaldehyde gel (Sambrook et al., 2001). Five μg of total RNA from each sample were treated with RQ1 RNase-Free DNase (Promega, Madison, WI, USA) to avoid DNA contamination. RNA quantity was verified with a Qubit® 3.0 fluorometer using a Qubit™ RNA Broad Range assay (Life Technologies, USA) while RNA quality was checked with the Agilent 2100 Bioanalyzer and RNA 6000 Nano Kit (Agilent Technologies, Wilmington, DE, USA). All the samples had an RNA Integrity Number (RIN) > 7 and were considered for library preparation.

2.9. Data accessibility All the raw reads were deposited into the Sequencing Read Archive (SRA) of NCBI with the accession number SRP183024. The BioProject ID of our data is PRJNA517955 and the BioSample accession is SAMN10848635. 2.10. Statistical analysis To perform the statistical analysis, Shapiro-Wilk and Barlett's tests were used to evaluate normality and homoscedasticity of oil bioaccumulation data on soft tissues of C. virginica. To determine statistical differences among groups, a one-way ANOVA with a significant level of 0.05 was used. When significant differences were observed, a Tukey's honestly significant difference (HSD) test was applied to determine which group had the highest oil concentration. Analyses were performed with SPSS software v23 (SPSS Inc., Chicago, IL, USA). 3. Results

2.7. Library preparation and sequencing

3.1. Super-light and light crude oil characterization

Prior to library preparation, 1 μg of total RNA from every sample was pooled for each treatment (Control (n = 12) and 50–200 μg/L (n = 12 per treatment) at day 14). Based on these pools, the biological averaged method (Kendziorski et al., 2003; Kendziorski et al., 2005) was used to compare gene expression levels between the four treatments (Control and 50, 100, and 200 μg/L, respectively). When biological variation is high and biologically averaged individuals are evenly pooled, this method reduces the biological variability and not only increases the power to detect differentially expressed genes, but it also retains biological information at a lower cost of sequencing (Biswas et al., 2013; Rajkumar et al., 2015). Libraries were prepared with the TruSeq® RNA Seq Sample Prep Kit v2 Set A (Illumina, San Diego, CA, USA). The mRNA purification and fragmentation, cDNA synthesis, end repair, adenylation, adapters ligation, PCR amplification, library validation, normalization, and pooling were performed with the TruSeq® RNA Seq Sample Prep v2 Guide (Illumina, Part # 15026495 Rev F), using the Low Sample Protocol. Libraries and PhiX control concentrations were 8 pM and 1%, respectively. Sequencing was performed using the Illumina MiSeq platform to generate about 76 bp pair-end raw reads.

Aliphatic hydrocarbons for super-light crude oil included components from nine to 32 carbons (C9–C32) while light crude oil was composed of smaller components (up to C27) (Fig. 2A). Pristane and Phytane were also detected in both oil compositions; C9, C13, and C28 were the most abundant components for super-light crude oil and C12, C15 and C27 were for light crude oil. For the PAHs, Dibenzothiophene and 1-Methylphenanthrene were the most abundant components for super-light crude oil while 1-Methylphenanthrene and Chrysene were for light crude oil (Fig. 2B). 3.2. Presence of aliphatic hydrocarbons and polycyclic aromatic hydrocarbons on sediment, in water and oysters On sediment, 15 components from C17 to C37 were detected for AHs, with a total concentration < 2.1 μg/g dry weight (Fig. 3A). For PAHs, seven components were detected with a total concentration < 0.028 μg/g dry weight and with Benz(a)anthracene and Chrysene as the most abundant (Fig. 3B). These AH and PAH concentrations would indicate that the zone had not been affected by human activities because PAH concentrations had been reported from 1 to 6.6 μg/g dry weight in contaminated sediments (Tolosa et al., 2005; Men et al., 2009). In water, AHs and PAHs were non-detectable or they were under the detection limit. Oysters sampled at the lagoon (in situ) and those of the control treatment through different times (0, 7 and 14 days) had no significant differences (P > 0.05) to AHs and PAHs concentrations. The mean concentration of AHs fluctuated from 10.45 to 11.49 μg/g d.w while that of PAHs fluctuated from 0.23 to 0.30 μg/g d.w. In the different experimental units (50–200 μg/L), AHs concentration showed an increase at day 7 to finally decrease significantly (P < 0.05) towards the end of the experiment (Fig. 4A). However, PAHs concentration in the experimental units showed a significant increase over time (P < 0.05) only to 100 and 200 μg/L. Unlike AHs, PAHs are difficult to process at cellular and tissue levels. For PAHs, the mean fold change for each comparison between T2 (14 days) and T0 (control) was: (a) 14.96 (50 μg/L); (b) 32.77 (100 μg/L); and (c) 39.92 (200 μg/L). Mean fold change between final and initial control was 1.17 (Fig. 4B). Finally, the most abundant PAHs in experimental treatments during exposure were

2.8. Bioinformatic analysis Adapters and reads with low quality (ambiguous bases (“N”) c < 30) were trimmed using the Trimmomatic software 0.32 (Bolger et al., 2014). Then, all clean reads were assembled with the de novo assembly software Trinity 2.2.0 (Grabherr et al., 2011). Functional annotation of assembled transcripts was performed with Trinotate software 3.0.1 (Haas et al., 2013) and databases Swissprot/Uniprot and PFAM-A, included in the Trinity package. Trinotate assigns the best BLAST result to each protein against the SwissProt/UniProt database and predictions for PFAM domains. The transcripts obtained were filtered to eliminate all those that did not have similarity with genes of eukaryotic origin. The resulting transcripts were used as reference for the analysis of differential expression (Robinson et al., 2010) to obtain putative products of transcription in silico using script run_Trinity_edgeR_pipeline.pl (dispersion value = 0.1 in edgeR). To control the expected false discovery rate (FDR), a multiplicity correction was 4

Comparative Biochemistry and Physiology, Part C 225 (2019) 108571

E.A. López-Landavery, et al.

Fig. 2. Percentage of compounds present in super-light and light crude oil used in Crassostrea virginica exposure experiment. (A) Aliphatic compounds (AHs) and (B) Aromatic compounds (PAHs).

58,217,157 of total assembled bases and N50 of 960 bp. Based on the longest isoform per gene, the median contig length was 353 bp with 36,238,438 of total assembled bases and N50 of 773 bp (Table 2).

phenanthrene, 1-methyl phenanthrene and pyrene (data not shown). 3.3. Illumina sequencing and de novo transcriptome assembly

3.4. Functional annotation of Crassostrea virginica transcriptome

The digestive gland transcriptome of the eastern oyster C. virginica was characterized with Illumina MiSeq platform. A total of 34,054,207 raw reads were obtained. After the quality control process with the software Trimmomatic 0.32 that included trimming of adapters and filtering of sequences with low quality, a total of 33,469,374 reads were retained. The clean reads obtained represented 98.28% of raw reads (Table 1). De novo assembly of clean reads was performed with Trinity software using default parameters to pair-end reads. Through assembly, check and elimination of redundancy, 61,356 genes and 86,947 transcripts were generated with 39.19% of GC content. Statistics based on all the transcripts reported a median contig length of 420 bp with

Trinotate pipeline was used to perform the functional annotation of the complete and assembled transcriptome of C. virginica. With this pipeline, Transdecoder identified coding regions within assembled transcriptome and searched an open reading frame (minimum length = 100 amino acids). The output of Trinotate pipeline was a file (Supplementary material 1) that showed different columns, such as gene_id, transcript_id, sprot_Top_BLASTX_hit, prot_coords, Pfam, eggNOG (Evolutionary Genealogy of Genes: Non-supervised Orthologous Groups), gene_ontology_blast, among others. EggNOG 5

Comparative Biochemistry and Physiology, Part C 225 (2019) 108571

E.A. López-Landavery, et al.

Fig. 3. Levels of aliphatic (AHs) (A) and aromatic (PAHs) (B) compounds detected on sediment from Laguna Morales.

and modification, 1.09%”, “Defense mechanisms, 0.72%” and “Cell motility, 0.24%”. Functional annotation from KEGG Pathway database showed that 5.98% of transcripts had significant matches with database and were assigned to 362 KEGG pathways based on KEGG Orthology identities (KO). KEGG Pathway database was composed of six categories; the most abundant pathways per category were: Metabolic pathways (Metabolism, 23.93%), ubiquitin-mediated proteolysis (Genetic information processing, 11.49%), PI3K-Akt signaling pathway (Environmental information processing, 6.51%), endocytosis (Cellular processes, 7.89%), NOD-like receptor signaling pathway (Organismal systems, 2.90%) and pathways in cancer (Human diseases, 5.43%). The top-five pathways per category are shown in Table 3. The Gene Ontology identities (GO ID) based on nr (Non-redundant) annotation database were used to assign and classify the biological process, molecular functions and cellular components associated with C. virginica transcripts. A total of 46,327 hits were obtained,

column included KEGG id (Kyoto Encyclopedia of Genes and Genomes) and COG id (Cluster of Orthologous Groups) terms while gene_ontology_blast column showed the relationship with the Gene Ontology (GO) terms. Based on sprot_Top_BLASTX_hit, 12.59% of transcripts were annotated with a minimum, maximum and average length of 297, 7896 and 816.51 bp, respectively. Also, 34.90% of BLASTX_hits were associated with length sequences from 297 to 500 bp, 38.15% between 500 and 1000 bp and 26.95% with sequences higher than 1000 bp. A total of 5.95% of transcripts were annotated from eggNOG database. Of these, 2.75% had hits with ENOG terms and 3.20% with COG terms. Based on twenty-five COG categories (Fig. 5), the most abundant groups were “General function prediction only, 22.97%”, “Posttranslational modification, protein turnover, chaperones, 7.90%” and “Cell cycle control, cell division, chromosome partitioning, 7.20%”. However, no transcript was associated with the “Extracellular structures” category. The least represented categories were “RNA processing 6

Comparative Biochemistry and Physiology, Part C 225 (2019) 108571

E.A. López-Landavery, et al.

Table 2 Statistical summary of the transcriptome assembly of the eastern oyster Crassostrea virginica. Total Trinity ‘genes’ Total Trinity transcripts GC (%) Stats based on all transcript contigs Contig N10 Contig N20 Contig N30 Contig N40 Contig N50 Median contig length Average contig Total assembled bases Stats based on only longest isoform per ‘gene’ Contig N10 Contig N20 Contig N30 Contig N40 Contig N50 Median contig length Average contig Total assembled bases

Table 1 RNA-Seq reads obtained before and after quality control.

Before After

50 μg/L

100 μg/L

200 μg/L

Total (%)

9,586,130 9,417,388

8,595,889 8,460,586

8,850,088 8,693,919

7,022,100 6,897,481

34,054,207 (100.00) 33,469,374 (98.28)

2532 1948 1549 1234 960 420 669.57 58,217,157 2374 1757 1354 1027 773 373 590.63 36,238,438

annotation, was to identify differentially expressed transcripts. For this purpose, RSEM (Li and Dewey, 2011) and bowtie2 (Langmead and Salzberg, 2012) were used to estimate expression values of each crude oil treatments (Control, 50 μg/L, 100 μg/L and 200 μg/L), based on alignment to Trinity assembly. Then, the differential expression analysis was performed with edgeR based on normalized read counts (Supplementary material 2). A total of 3196 transcripts were expressed differentially among four treatments with FDR < 0.05 and FC ≥ 4 (FDR: False Discovery Rate, FC: Fold Change). A subset of 1022 transcripts expressed differentially (FDR < 0.001 and FC ≥ 4, Fig. 7) between a couple of four libraries were selected for the following analysis. From this subset, the highest number of differentially expressed transcripts (up + down) corresponded to the treatments vs control comparisons. In this sense, the comparison between 100 μg/L vs control generated 496 up/down-regulated transcripts while the highest number of up/down-regulated transcripts between treatments was 305 (100 μg/ L vs 200 μg/L, Fig. 8A). Compared with control, the largest number of transcripts that differentially expressed exclusively was obtained in 100 μg/L (Fig. 8B). To know the overrepresented biological functions of transcripts expressed differentially, a gene ontology enrichment analysis was performed between control and three crude oil treatments. The enrichment analysis was performed with GOseq R package (Young et al., 2010). Compared to control, transcripts in 50 μg/L were enriched in the extracellular region, mitochondrial inner membrane and cytochrome-c oxidase activity. In 100 μg/L, the enriched functions were associated with the structural constituent of ribosome, innate immune response, complement activation, among others. Negative regulation of ERK1 and ERK2 cascade, negative regulation of fibroblast proliferation and transposition DNA-mediated were among the most significantly enriched functions in 200 μg/L treatment (Fig. 9). Interestingly in 200 μg/ L, a central component in energy metabolism, such as cytochrome c oxidase (subunit 1, 2 and 3), was down-regulated concerning the levels found in 50 and 100 μg/L (Supplementary material 3).

Fig. 4. (A) AHs and (B) PAHs levels (μg/g) on a dry weight basis detected in whole soft tissues of Crassostrea virginica exposed to different concentrations of crude oil mixture for 14 days (n = 3 pools/treatment, 15 oysters per pool, except to LS (20 oysters)). T0: 0 days; T1: 7 days; T2: 14 days; LS: Lagoon sampling. Data shown as mean ± SD.

Control

61,356 86,947 39.19

representing 35 functional categories reported to the second level of GO and related with 24,320 transcripts (28%). Cellular (9%), single-organism (9%) and metabolic (8%) processes were the main categories in Biological Process. In Molecular Functions, binding (59%) was the most abundant category, while cell (16%), cell part (16%) and organelle (15%) were for Cellular Components (Fig. 6). Other categories of interest in this study were biological regulation (8%), stimulus response (6%), immune system process (5%), among others.

3.6. Transcripts associated with metabolism of xenobiotics and immune system A wide range of transcripts related to phase I and II enzymes and transporters, contributing to metabolism of xenobiotics and antioxidants, were obtained from the annotation process. A total of 146 transcripts were associated with Cytochrome P450, including all families of this enzyme. However, up-regulated transcripts (200 μg/L vs

3.5. Differential expression and enrichment gene function The next goal, after C. virginica transcriptome assembly and 7

Comparative Biochemistry and Physiology, Part C 225 (2019) 108571

E.A. López-Landavery, et al.

Fig. 5. Classification of Cluster of Orthologous Group (COG) based on Crassostrea virginica transcriptomic response.

transcripts corresponded to GST family members. From this group, two transcripts were induced in response to crude oil exposure, GST theta 1 (FC = 6.27, 100 μg/L, FDR = 0.030) and GST mu 3 (FC = −156.16, 50 μg/L, FDR = 0.031). The most abundant members were GSTA and microsomal GST3. On the other hand, 76 transcripts coded for Sulfotransferases family members (SULT). For this enzyme, only was one transcript associated to SULT1B1 up-regulated in all crude oil

control) corresponded only to Cytochrome P450 4B1 and 3A24 with fold changes of 11.68 and 6.23, respectively. Conversely, Cytochrome P450 2B19 and 3A29 were the most down-regulated transcripts in 200 μg/L with fold changes of −38.17 and − 26.85, respectively. The most abundant family was CYP2. Phase II enzymes of xenobiotics metabolism were associated with Glutathione S-transferases and Sulfotransferases. A total of 62

Table 3 Top-five pathways per category obtained from KEGG database for the transcriptome of the eastern oyster Crassostrea virginica. Category (%)

Pathway

KO ID

% (inside category)

Metabolism (26.68)

Metabolic pathway Biosynthesis of secondary metabolites Biosynthesis of antibiotics Purine metabolism Microbial metabolism in diverse environments Ubiquitin-mediated proteolysis Spliceosome Protein processing in endoplasmic reticulum RNA transport Ribosome PI3K-Akt signaling pathway MAPK signaling pathway Ras signaling pathway Neuroactive ligand-receptor interaction Calcium signaling pathway Endocytosis Focal adhesion Lysosome Regulation of actin cytoskeleton Cell cycle NOD-like receptor signaling pathway Retrograded endocannabinoid signaling Insulin signaling pathway Platelet activation Chemokine signaling pathway Pathways in cancer Human papillomavirus infection HTLV-I infection Proteoglycans in cancer Epstein-Barr virus infection

01100 01110 01130 00230 01120 04120 03040 04141 03013 03010 04151 04010 04014 04080 04020 04144 04510 04145 04819 04110 04621 04723 04910 04611 04062 05200 05165 05166 05205 05169

23.93 7.13 4.43 3.90 3.60 11.49 11.25 10.51 8.80 8.07 6.51 5.37 4.88 4.88 4.23 7.89 6.86 6.17 6.17 5.83 2.90 2.67 2.44 2.32 2.32 5.43 4.68 3.76 3.59 3.43

Genetic information processing (8.18)

Environmental information processing (12.29)

Cellular processes (11.67)

Organismal systems (17.25)

Human diseases (23.93)

8

Comparative Biochemistry and Physiology, Part C 225 (2019) 108571

E.A. López-Landavery, et al.

Fig. 6. Gene Ontology (GO) analysis of Crassostrea virginica transcriptome reported to the second level of GO. (A) Biological Process, (B) Molecular Functions and (C) Cellular Components.

treatments when it was compared with control with a maximum fold change of 11.07 (50 μg/L). Superoxide dismutase (SOD [CueZn]), a protective enzyme against oxidative stress was identified from C. virginica transcriptome. A total of 13 transcripts were associated with this enzyme, and one of them was highly up-regulated with a maximum fold change of 8.48 in 50 μg/L when compared with control. Stress response enzymes, such as heat shock proteins (HSPs) were also identified. A total of 322 transcripts corresponded to HSP family members. Some transcripts of HSP90 and HSP70 family members were up-regulated in a dose-dependent way until 100 μg/L. Compared to control, the highest fold change for HSP70 was 22.80 while it was 41.53 for HSP90. Transcripts associated with the immune system were identified and expressed differentially from C. virginica transcriptome in response to crude oil exposure. Some of these transcripts coded for C1q complement containing proteins, c-binding like, scavenger receptors, fibrinogen-related proteins, collagen chain-like, polyadenylated binding, among others, which together are classified as pattern recognition receptors or PRRs. Although mixed expression results for C1q complement were obtained, in general, these transcripts were up-regulated in a dose-dependent way until 100 μg/L while in 200 μg/L expression they decreased to similar or slightly higher levels compared with control. 4. Discussion Crassostrea virginica is a non-model species, and to date, a reference genome fully cured has not been developed yet. This study obtained the digestive gland transcriptome in response to crude oil exposure, generating 86,947 transcripts with N50 of 960 bp and average length contigs of 669.57 bp (considering all transcripts) and 590.63 bp (considering the longest isoform per gene). These metrics about average length contigs agreed with those reported by Droege and Hill (2008). Other studies of RNA-Seq reported in marine invertebrates include Ruditapes philippinarum (Ghiselli et al., 2012), Octopus vulgaris (Zhang et al., 2012), C. virginica (Lüchmann et al., 2014) Jenny et al., 2016), Cristaria plicata (Patnaik et al., 2016), and Sinonovacula constricta (Wang et al., 2016), among others. All these studies, including ours, have reported a high recovery of reads and mean contig lengths higher than 350 bp. These characteristics and the low cost per reads should allow more significant sequencing efforts, improving genomic resources and functional annotation of non-model species, such as C. virginica. Annotation levels of C. virginica transcripts fluctuated from 5.95% (eggNOG) to 28.00% (GO). Furthermore, the annotation level from BLASTX was 12.59%. These low levels were likely associated with limited information about genomes and transcriptomes of the eastern oyster and related species. Another cause for this result might have been associated with the database that was used for annotation, which in our case was UniProt/SwissProt, except for GO. NR non-redundant is another database frequently used for annotation of massive sequencing data. SwissProt is the section of UniProt knowledgebase, and it contains publicly available manually annotated protein sequences from a broad spectrum of organisms (Boutet et al., 2016) while NR non-redundant is a more extensive database used for routine searches and includes UniProt/SwissProt, PIR and PRF databases (Schuler, 2001). Transcriptome studies reporting higher annotation levels using NR non-redundant database in species of the Crassostrea and other genus include C. virginica (23.14%, Zhang et al., 2014; 19.80%, McDowell et al., 2014), C. brasiliana (47.5%, Lüchmann et al., 2014) and winged pearl oyster Pteria penguin (41.47%, Li et al., 2017) while a similar pattern of level annotation (GO > SwissProt/UniProt > COG/KEGG) with this study was reported in Pearl oyster Pinctada fucata (Wang et al., 2016). 9

Comparative Biochemistry and Physiology, Part C 225 (2019) 108571

E.A. López-Landavery, et al.

Fig. 7. Heat-map for differential expression levels in Crassostrea virginica exposed to different concentrations of crude oil. The points correspond to the beginning (Ctrl, 0 μg/L, day 0) and end of the hydrocarbons exposure (50–200 μg/L, day 14). Expression values are represented as log2 (CPM + 1), with CPM representing counts per million.

activity. These categories were also present in the response of Sacostrea glomerata exposed to Pyrene and Fluoranthene (Ertl et al., 2016). Likewise, the number of metabolic pathways found in this study was higher than that reported in adults of Perna viridis exposed to different chemical and physical stressors (Leung et al., 2014), oyster C. brasiliana exposed to diesel, Phenanthrene and sewage (Lüchmann et al., 2014) and embryos of P. viridis exposed to Benzo(a)pyrene (Jiang et al., 2016). As a representative member of the marine invertebrates from the Gulf of Mexico and based on its life history (experiments dramatic environmental stress), C. virginica is an excellent candidate to study the effects of crude oil on marine invertebrates as shown in the functional annotation analysis in this study.

Regardless of the database used and considering the importance of annotation for downstream analyses and further studies, another critical aspect to consider in the future (with more complete databases) is that current unblasted transcripts might be good candidates for gene discovery. The GO and KEGG analysis from C. virginica transcripts supported the idea about diversity of the family genes involved in the different process, functions and metabolic pathways in response to biotic and abiotic stressors. Zhang et al. (2012) proposed this idea about stress response plasticity in members of Crassostrea genus. Among the primary biological processes and functions, this study found biological regulation, stimulus response, immune system process, binding and catalytic

10

Comparative Biochemistry and Physiology, Part C 225 (2019) 108571

E.A. López-Landavery, et al.

Fig. 8. Transcripts differentially expressed (TDE) in response to crude oil exposure in Crassostrea virginica. (A) Numbers of up- and down-regulated TDE in each crude oil treatments comparison. The x-axis shows six comparisons and y-axis represents the total number of TDE. (B) Venn diagrams of TDE (up-regulated) among crude oil treatments compared with control.

transcripts was higher than those down-regulated. Nonetheless, while comparing treatments only, the relationship was reversed (downregulated > up-regulated) and dose-dependent. The above might imply that hydrocarbons inhibited the expression of transcripts significantly in oysters, modulating their biological functions and metabolic pathways. This implication would be related to the gene ontology enrichment analysis, which showed that GO categories in control were similar to 50 μg/L and varied from 100 to 200 μg/L. Groups associated with 200 μg/L included negative regulation of ERK1 and ERK2 cascade, fat cell differentiation and fibroblast proliferation. ERK1 and ERK2 belonged to mitogen-activated family protein group, and they were

Previous studies have shown the adverse effects of the different oil crude components on mollusk physiology (Snyder, 2000; Frouin et al., 2007; Lopes et al., 2012). These effects had been correlated in a dosedependent manner but up to their tolerance limit. For example, the magnitude of ratio change affected between double-stranded DNA and total DNA in the ovary of Chlamys farreri was correlated with Benzo[a] pyrene concentrations (Tian et al., 2013). Furthermore, a more considerable histological damage in digestive gland was reported in C. brasiliana exposed to 100 μg/L of Phenanthrene compared with 1000 μg/L and control (dos Reis et al., 2015). In this study, while comparing treatments and control, the number of up-regulated 11

Comparative Biochemistry and Physiology, Part C 225 (2019) 108571

E.A. López-Landavery, et al.

Fig. 9. Gene enrichment analysis for transcripts differentially expressed in response to crude oil exposure in Crassostrea virginica at day 14.

Scavenger receptor, Acetylcholinesterase (AChE), among others. In this sense, biotransformation and excretion of PAHs are composed of three phases. Phase I includes enzymes of the Cytochrome P450 family that through processes of oxidation, reduction and hydrolysis, they act on hydrophobic compounds converting them more water soluble. In Phase II, enzymes like Glutathione S-Transferase and sulfotransferases link metabolites (from phase I) with the polar endogenous substrates, increasing their solubility and facilitating their excretion. In Phase III, specialized transport like multixenobiotic transporter system leads the elimination of conjugates from the cell (Amiard-Triquet and Amiard, 2013). Based on the differential expression analysis (FDR < 0.001 and FC ≥ 4), we got mixed results for Cytochrome P450 gene expression. Up-regulated transcripts corresponded to Cytochrome P450 4B1 and 3A24, both overexpressed in 200 μg/L. CYP3 and CYP4 enzyme families have been known for their role in pesticide metabolism (Karatolos et al., 2011). In bivalve mollusks, overexpression of CYP3 and CYP4 was reported in S. glomerata exposed to Pyrene and Fluoranthene (Ertl et al., 2016) and in Nodipecten nodosus exposed to 200 μg/L of Phenanthrene (Piazza et al., 2016). Moreover, in polychaete Nereis virens, CYP342A1 was overexpressed in response to crude oil and Benz(a)anthracene exposure (Rewitz et al., 2004). CYP342A1 has been related to CYP4 enzyme family and along with CYP4BB1 mediate Pyrene metabolism (Jørgensen et al., 2005; Rewitz et al., 2006). Conversely, Cytochrome P450 2B19 and 3A29 were the most down-regulated transcripts in 200 μg/L. CYP2 and CYP3 have been involved in defense against toxic compounds in bivalves (Miao et al., 2011), and down-regulated expression of CYP enzymes under oil components exposure had been reported in embryos of C. gigas (Nogueira et al., 2017). Although CYP1 has a primary role in PAHs metabolism in vertebrates (Petrulis et al., 2001; Brown-Peterson et al., 2015) and it had been recently identified in C. virginica (Jenny et al., 2016), in our database it was not differentially expressed under differential expression parameters. Based on our results, some members of CYP4 and CYP3 families might have been induced by crude oil exposure; therefore, they were involved in PAHs metabolism in C. virginica. During phase II, metabolites from phase I can be conjugated to endogenous substrates through the enzymes such as Glutathione Stransferases (GSTs) and Sulfotransferases (SULTs) (Bustamante et al., 2012; Lacroix et al., 2014). GST family members are involved in the conjugation of glutathione with xenobiotics containing an electrophilic structure and in the defense against oxidative damage of DNA and lipids (Van der Oost et al., 2003; Blanchette et al., 2007). Dealing with electrophilic compounds is necessary because they can react with macromolecules, such as DNA, RNA, and proteins that control cell

involved in functions such as cell survival and apoptosis (Chen et al., 2001). Specifically, negative regulations of ERK1 and ERK2 were associated with cell survival, and in this study, they might have been supported by PI3K-Akt signaling pathway (Phosphatidylinositol 3-Kinase – Protein Kinase B) abundance from KEGG database. Although PI3K-Akt signaling pathway functions as a cell survival signal to protect cell from apoptosis, and it has been proposed as a power system to understand survival signal in response to stress in C. gigas (Epelboin et al., 2016), more studies are necessary to assess the significance of alteration on specific gene expression in these metabolic pathways in response to crude oil exposure. In this sense, another important aspect from the GO enrichment to 200 μg/L was the non-presence of the cytochrome c oxidase category. It was supported for the differential expression analysis, where transcripts related with the three cytochrome c oxidase subunits (1, 2 and 3) were up-regulated only to 50 and 100 μg/ L. Cytochrome c oxidase is involved in oxidative phosphorylation, which is part of the energy metabolism through electron transport chain (ETC, Tymoczko et al., 2013). Interestingly, the transcript identified as putative ferric-chelate reductase 1 homolog was down-regulated only at 200 μg/L (FC = −8.29, FDR = 6.33E−04) and has been reported to belong to the cytochrome b protein family (Vargas et al., 2003). The FeeS species and members of the cytochrome b family are components of complex III, which catalyzes the transport of reducing equivalents from coenzyme Q to cytochrome c (Bhagavan and Ha, 2015). The expression pattern found in this study for cytochrome c oxidase and putative ferric-chelate reductase 1 homolog may indicate an alteration in the energetic metabolism to 200 μg/L. Recently, Kirby et al. (2019) suggested that acute oil exposure in young adults of Mahi-mahi (Coryphaena hippurus) inhibited mitochondrial function, specifically in complex III of ETC. Furthermore, although cytochrome c oxidase is associated with the complex IV, the flow of electrons is not only directional from complex I-II to complex IV but also an interrelationship is found between complex III and IV (Garret and Grisham, 2017). Based on the above and the enrichment analysis, it may be suggested that exposure to 200 μg/L affects the mitochondrial function of C. virginica negatively. Different biomarkers have been used to assess chemical contamination in aquatic media, including hydrocarbons (Chapman et al., 2013; Amiard-Triquet and Berthet, 2015), which have been classified in (a) defense or damage biomarkers (Davies and Vethaak, 2012) or (b) acute-term or long-term effect biomarkers (Martínez-Gómez et al., 2010). In this study, transcripts corresponding to both categories were expressed differentially, such as Cytochrome P450, Metallothionein, Glutathione S-transferase (GST), Superoxide dismutase (SOD), 12

Comparative Biochemistry and Physiology, Part C 225 (2019) 108571

E.A. López-Landavery, et al.

SOD in cellular protection against oxidative stress produced by crude oil metabolism. Heat shock proteins (HSPs), a family of plasma proteins, allow the detection of environmental stress and help to modulate the organism response through protein folding, translocation to cellular compartments and the degradation of misfolded proteins (Srivastava, 2002; Parcellier et al., 2003; Wang et al., 2003). These proteins are present in homeostasis condition in the cell and induced for thermal shock, anoxia, salinity, oxidative stress and xenobiotic components (CossuLeguille and Vasseur, 2013). In this study, different HSPs (HSP70 and HSP90 families) were up-regulated in a crude oil dose-dependent way. Interestingly, the highest expression levels were recorded in 50 and 100 μg/L treatments while the expression was down-regulated in 200 μg/L. Compared to control and 200 μg/L, the highest fold change for HSP70 was 22.80 while for HSP90, it was 41.53. Studies in marine invertebrates suggesting the overexpression of HSP in response to hydrocarbons exposure have been reported in C. gigas (Boutet et al., 2004), C. brasiliana (Lüchmann et al., 2014) and S. glomerata (Ertl et al., 2016). The decrease in the HSPs in 200 μg/L may be associated to the high concentration of crude oil and exposure time (14 days), reflected probably in the enrichment analysis profile (Fig. 6). In the case of HSP70, Köhler et al. (2001) reported that high concentrations of pollutants might inhibit the biomarker response and decrease the adaptive capacity strongly to deal with these pollutants. Although HSP70 has been suggested as a biomarker in aquatic toxicology (Cruz-Rodrıguez and Chu, 2002), more research is necessary to characterize family members and to evaluate the effect of exposure time (acute vs chronic). Invertebrates, as vertebrates, recognize and respond to pollutant presence, such as hydrocarbons. In this study, exposure of C. virginica to crude oil showed not only activation of physiological functions associated with oil detoxification, stress protein, and stress response but also activation of the immune system. It is well known that suppression activation of the immune system will depend on pollutant concentration and exposure time. Here, a group of transcripts coding for different components of the immune system and classified as PRRs were upregulated. PRRs recognize the conserved pathogen-associated molecular patterns, activating intracellular signaling pathways and triggering the production of antimicrobial effectors (Akira et al., 2006; Yang et al., 2014; Wang et al., 2018). Up-regulation of PRRs in a dosedependent way under crude oil exposure suggests that this species responds effectively until 100 μg/L. Also, down-regulation in 200 μg/L with respect to 50 and 100 μg/L would indicate that this concentration might suppress the immune system, making the oysters more susceptible to pathogens. It might have been related to 2.50% of mortality in 200 μg/L at the end of the experiment and the presence of Perkinsus marinus (no mortality in other treatments, data not shown). In this sense, although the presence of P. marinus was detected through sequence homology, the mortality level at 200 μg/L was low, and we suggest that it was mainly due to the level of hydrocarbons present in such treatment. This fact would be supported by the level of bioaccumulation that we found (Fig. 4B), the enrichment categories of functional annotation (Fig. 9) and the analysis based on KEGG (Table 3). In the enrichment analysis at 200 μg/L we find categories associated with negative regulation and cell proliferation; whereas in the analysis based on KEGG, the top three categories were metabolism (26.68%), human diseases (23.93%) and organismal systems (17.25%). In the human diseases category, the most abundant pathways were associated with cancer (5.43%). It has been reported that the presence of PAHs not only induces immunotoxicity, mutagenicity and carcinogenicity (Neff, 1979), but it has also been shown to compromise protein regulation, immune response and stress response systems (Lüchmann et al., 2012). PAHs components of crude oil were bioaccumulated in soft tissues of C. virginica in a dose-dependent way. This dose-dependent bioaccumulation process was associated with the high affinity between fat content and hydrocarbons, which are chemically hydrophobic and their distribution in the different tissues. It means that once absorbed by the

growth and other essential functions (Hampel et al., 2016). This study identified a wide range of GST family members, but in general, their expression was not induced significantly between treatments. Only were two transcripts expressed differentially with respect to control, GST theta 1 and GST mu 3. GST theta 1 was expressed in a dose-dependent manner up to 100 μg/L, which agreed with that reported in C. brasiliana and C. virginica (Lüchmann et al., 2014; Jenny et al., 2016). Also, down-regulation of GST has been reported in Chlamys farreri exposed to Benzo[a]pyrene (Cai et al., 2014). Other contrasting results about gene expression of the GST family have been reported in mussel, scallop and oysters exposed to crude oil compounds (Akcha et al., 2000; Jing-Jing et al., 2009; Ertl et al., 2016; Piazza et al., 2016; Zacchi et al., 2017). According to contrasting results, Boutet et al. (2004) showed variations in the expression of GST family members dependent on time exposure (acute vs chronic) on C. gigas exposed to hydrocarbon mixture (0.1%) and under laboratory conditions. Other factors that could influence the GST expression patterns reported in marine invertebrates might have been associated with the type of PAH (individual vs mixture vs crude oil), PAH concentration, tissues and exposure conditions (laboratory vs field). Based on our results, GST theta represents a potential biomarker in aquatic toxicology, and it agrees with a significant relationship between fold induction and higher PAH content in the digestive gland reported by Jenny et al. (2016) in the same species. Sulfotransferases catalyze reactions of sulfation and conjugation, which are important for metabolism of endogenous (steroid hormones and bile salts) and xenobiotic (aromatic amines and PAHs) compounds, especially those with hydroxyl groups (Kennedy and Tierney, 2013). In this study, different members of the sulfotransferase family were identified in C. virginica transcriptome; although they were generally not significantly induced for crude oil exposure, the tendency was down-regulated. However, only was a transcript up-regulated (SULT1B1) in all hydrocarbon treatments. Similar to our results, this enzyme (SULT1B1-like) was also up-regulated in N. nodosus exposed to 200 μg/L of Phenanthrene for 24 and 96 h (Piazza et al., 2016). In other studies, up-regulation of only one transcript associated with SULT1C1 was reported in C. brasiliana in response to Phenanthrene (Lüchmann et al., 2014) while up-regulation of SULT-like transcripts was stated in response to Phenanthrene and different salinity exposure (Zacchi et al., 2017). Contrary to the above, down-regulation of SULT1C1 was reported in C. gigas after 21 days of hydrocarbon exposure (Boutet et al., 2004). In vertebrates, SULT1B1 has been reported as one of the essential SULT enzymes involved in xenobiotic metabolism (Riches et al., 2009), and although little is known in marine invertebrates about its function, our results suggest its participation in detoxification of xenobiotics in C. virginica. Reactive oxygen species (ROS) are constantly produced in aerobic organisms as part of their basal metabolism, including cellular respiration or phagocytic activity (Valavanidis et al., 2006). Another source of ROS corresponds to the activity of oxidative enzymes during xenobiotic metabolism. Oxidative enzymes, such as Cytochrome P450 or dioxygenase, produce free radicals with the potential to generate hydrogen peroxide (Roméo and Giambérini, 2013). ROS are toxic and dangerous to the cell, causing lipid peroxidation and DNA damage, inducing the production of antioxidant enzymes, such as Superoxide dismutase (SOD), Catalase and Glutathione peroxidase (Noh et al., 2015; Wang et al., 2016). In this study, one SOD transcript was highly induced in response to crude oil exposure (log2FC = 8.48, 50 μg/L). Studies reporting up-regulation of SOD in oysters exposed to individual components of crude oil, mixture or only crude oil, include C. brasiliana, C. gigas, and C. virginica (Boutet et al., 2004; Lüchmann et al., 2014; Jenny et al., 2016; Zacchi et al., 2017). The activity of antioxidant enzymes might reflect the antioxidant capacity of an organism in response to particular stress, and it would be related with the activity of enzymes from xenobiotics biotransformation phase, such as Cytochrome P450 (Jenny et al., 2016). Our results showed up-regulation of Cytochrome P450 family members and SOD, suggesting a main role of 13

Comparative Biochemistry and Physiology, Part C 225 (2019) 108571

E.A. López-Landavery, et al.

Acknowledgments

organisms, its elimination is difficult and favors its retention in storage organs, such as the digestive gland, added to the fact that a low PAH metabolism capacity, reported in bivalve mollusks (Snyder, 2000; Frouin et al., 2007; Gauthier et al., 2014), causes a more efficient bioaccumulation process. Moreover, as a consequence of the inherent hydrophobicity of PAHs, which increases with molecular weight and ring numbers, their incorporation and solubility in fatty tissues depend mainly on physicochemical interaction (Freire and Labarta, 2003). For example, Prada-Ríos and Zambrano-Ortiz (2006) reported that in the bivalve Anadara tuberculosa a higher accumulation of Fluoranthene (high molecular weight) than that of Naphthalene (low molecular weight) had been observed, demonstrating the enormous capacity of bivalves to incorporate PAHs of high molecular weight. Studies reporting bioaccumulation of PAH crude oil components in oysters include C. gigas (Boutet et al., 2004), C. brasiliana (dos Reis et al., 2015; Piazza et al., 2016) and S. glomerata (Ertl et al., 2016). Based on the above and the sessile lifestyle of juvenile and adult oysters and other mollusks that favor pollutant concentration, these species are widely used as sentinel organisms in aquatic toxicology. When the detection of differentially expressed genes or transcripts is one of the main objectives in RNA-Seq studies, the mathematical averaged method is preferred over the biological averaged one mainly because of its higher statistical resolution (Schurch et al., 2016). However, when it is difficult to obtain replicates for RNA-Seq, either because of costs, the scale of the study or the availability of individuals, the biological averaged method represents a valid alternative to the mathematical averaged one (Biswas et al., 2013). In this sense, the biological averaged method can sufficiently control biological variation so that differences in gene expression can be reliably detected (Kendziorski et al., 2003; Kendziorski et al., 2005; Rajkumar et al., 2015). Moreover, to minimize the effect of pooling bias that refers to the differences between the measured values in the pool and their corresponding individual replicates, it is necessary to include a considerable number of samples in the pool with an optimum from 8 to 20 and use a detection method sensitive and specific like edgeR (Robinson et al., 2010; Biswas et al., 2013; Rajkumar et al., 2015). In this study, each of the pools was made up of 12 individuals, and edgeR was used to detect differences in transcripts expression. In studies of RNA-Seq based on pooled samples, using between 10 and 20 individuals per pool, differentially expressed genes have been correctly detected (Xu et al., 2012). Furthermore, Biswas et al. (2013) reported that concordance levels in the detection of differentially expressed genes using edgeR was > 95% for the top 4000 features and > 88% for the top 400 features between the biological and mathematical averaged methods. In conclusion, the transcriptomic response in this study for C. virginica showed the activation of metabolic pathways and induction of specific genes in response to crude oil exposure and PAHs bioaccumulation. These genes were associated with different stages of xenobiotic metabolism and detoxification, which were generally associated with different categories: phase I and phase II biotransformation enzymes, lipid metabolism, induction of enzymes with antioxidant capacity and chaperones. These results contribute to expanding the current genomic resources for this representative species of the Gulf of Mexico. These data will help to develop new studies in aquatic toxicology focused on knowledge in depth of metabolic pathways that together with proteomic approaches will allow to obtain a complete idea about the eastern oyster response to crude oil. Supplementary data to this article can be found online at https:// doi.org/10.1016/j.cbpc.2019.108571.

The authors are grateful to Manuel Roel and Braulio Soto for their assistance during the experimental phase and Varinka Saenz for map elaboration and salinity data variation of Laguna Morales; to Diana Fischer for editorial services in English. This research was funded by the National Council of Science and Technology of Mexico – Mexican Ministry of Energy and Hydrocarbon Trust, project 201441. This is a contribution of the Gulf of Mexico Research Consortium (CIGoM). References Aguilar, C.A., Montalvo, C., Rodríguez, L., Cerón, J.G., Cerón, R.M., 2012. American oyster (Crassostrea virginica) and sediments as a coastal zone pollution monitor by heavy metals. Int. J. Environ. Sci. Technol. 9 (4), 579–586. Akcha, F., Izuel, C., Venier, P., Budzinski, H., Burgeot, T., Narbonne, J.F., 2000. Enzymatic biomarker measurement and study of DNA adduct formation in benzo [a] pyrene-contaminated mussels, Mytilus galloprovincialis. Aquat. Toxicol. 49 (4), 269–287. Akira, S., Uematsu, S., Takeuchi, O., 2006. Pathogen recognition and innate immunity. Cell 124 (4), 783–801. Amiard-Triquet, C., Amiard, J.C., 2013. Introduction. In: Amiard-Triquet, C., Amiard, J.C., Rainbow, P.S. (Eds.), Ecological Biomarkers: Indicators of Ecotoxicological Effects. CRC Press, Boca Raton, FL, USA, pp. 1–14. Amiard-Triquet, C., Berthet, B., 2015. Individual biomarkers. In: Amiard-Triquet, C., Amiard, J.C., Mouneyrac, C. (Eds.), Aquatic Ecotoxicology: Advancing Tools for Dealing With Emerging Risks. Academic Press, USA, pp. 153–182. Anderson, K., Taylor, D.A., Thompson, E.L., Melwani, A.R., Nair, S.V., Raftos, D.A., 2015. Meta-analysis of studies using suppression subtractive hybridization and microarrays to investigate the effects of environmental stress on gene transcription in oysters. PLoS One 10 (3), e0118839. Benjamini, Y., Hochberg, Y., 1995. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. 57, 289–300. https://doi.org/ 10.2307/2346101. Bhagavan, N.V., Ha, C.E., 2015. Essentials of Medical Biochemistry: With Clinical Cases. Academic Press. Biswas, S., Agrawal, Y.N., Mucyn, T.S., Dangl, J.L., Jones, C.D., 2013. Biological averaging in RNA-seq. arXiv preprint (arXiv:1309.0670). Blackburn, M., Mazzacano, C.A.S., Fallon, C., Black, S.H., 2014. Oil in our Oceans: A Review of the Impacts of the Oil Spills on Marine Invertebrates. The Xerces Society for Invertebrate Conservation. Portland, OR, USA 152 pp. Blanchette, B., Feng, X., Singh, B.R., 2007. Marine glutathione S-transferases. Mar. Biotechnol. 9 (5), 513–542. Bolger, A.M., Lohse, M., Usadel, B., 2014. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30 (15), 2114–2120. Boutet, I., Tanguy, A., Moraga, D., 2004. Response of the Pacific oyster Crassostrea gigas to hydrocarbon contamination under experimental conditions. Gene 329, 147–157. Boutet, E., Lieberherr, D., Tognolli, M., Schneider, M., Bansal, P., Bridge, A.J., Poux, S., Bougueleret, L., Xenarios, I., 2016. UniProtKB/Swiss-Prot, the manually annotated section of the UniProt KnowledgeBase: How to use the entry view. In: Edwards, D. (Ed.), Plant Bioinformatics. Methods in Molecular Biology, vol. 1374. Humana Press, New York, USA, pp. 23–54. Brown-Peterson, N.J., Krasnec, M., Takeshita, R., Ryan, C.N., Griffitt, K.J., Lay, C., Mayer, G.D., Bayha, K.M., Hawkins, W.E., Morris, J., 2015. A multiple endpoint analysis of the effects of chronic exposure to sediment contaminated with Deepwater Horizon oil on juvenile southern flounder and their associated microbiomes. Aquat. Toxicol. 165, 197–209. Bustamante, P., Luna-Acosta, A., Clemens, S., Cassi, R., Thomas-Guyon, H., Warnau, M., 2012. Bioaccumulation and metabolization of 14C-pyrene by the Pacific oyster Crassostrea gigas exposed via seawater. Chemosphere 87 (8), 938–944. Caffrey, J.M., Hollibaugh, J.T., Mortazavi, B., 2016. Living oysters and their shells as sites of nitrification and denitrification. Mar. Pollut. Bull. 112 (1–2), 86–90. Cai, Y., Pan, L., Hu, F., Jin, Q., Liu, T., 2014. Deep sequencing-based transcriptome profiling analysis of Chlamys farreri exposed to benzo [a] pyrene. Gene 551 (2), 261–270. Chapman, P.M., Wang, F., Caeiro, S.S., 2013. Assessing and managing sediment contamination in transitional waters. Environ. Int. 55, 71–91. Chen, Z., Gibson, T.B., Robinson, F., Silvestro, L., Pearson, G., Xu, B.E., Wright, A., Vanderbilt, C., Cobb, M.H., 2001. MAP kinases. Chem. Rev. 101 (8), 2449–2476. Cossu-Leguille, C., Vasseur, P., 2013. Aquatic biomarkers. In: Férard, J.F., Blaise, C. (Eds.), Encyclopedia of Aquatic Ecotoxicology. Springer, Dordrecht, Netherlands, pp. 49–66. Cruz-Rodrıguez, L.A., Chu, F.L.E., 2002. Heat-shock protein (HSP70) response in the eastern oyster, Crassostrea virginica, exposed to PAHs sorbed to suspended artificial clay particles and to suspended field contaminated sediments. Aquat. Toxicol. 60 (3–4), 157–168. Davies, I.M., Vethaak, D., 2012. Integrated marine environmental monitoring of chemicals and their effects. ICES Cooperative Research Report 315. dos Reis, I.M., Mattos, J.J., Garcez, R.C., Zacchi, F.L., Miguelão, T., Flores-Nunes, F., Toledo-Silva, G., Sasaki, S.T., Taniguchi, S., Bícego, M.C., Cargnin-Ferreira, E., Bainy, A.C.D., 2015. Histological responses and localization of the cytochrome P450 (CYP2AU1) in Crassostrea brasiliana exposed to phenanthrene. Aquat. Toxicol. 169,

Declaration of Competing Interest The authors declare that they have no conflict of interest.

14

Comparative Biochemistry and Physiology, Part C 225 (2019) 108571

E.A. López-Landavery, et al.

(12), 4252–4257. Kennedy, C.J., Tierney, K.B., 2013. Xenobiotic protection/resistance mechanisms in organisms. In: Laws, E.A. (Ed.), Environmental Toxicology: Selected Entries From the Encyclopedia of Sustainability Science and Technology. Springer, New York, NY, pp. 689–721. Kennish, M. J. 1991. Ecology of Estuaries: Anthropogenic Effects, vol. 1. CRC Press, Boca Raton, FL, USA. 494 pp. Kirby, A.R., Cox, G.K., Nelson, D., Heuer, R.M., Stieglitz, J.D., Benetti, D.D., Grosell, M., Crossley II, D.A., 2019. Acute Crude Oil Exposure Alters Mitochondrial Function and ADP Affinity in Cardiac Muscle Fibers of Young Adult Mahi-mahi (Coryphaena hippurus). Comparative Biochemistry and Physiology Part C, Toxicology & Pharmacology. Köhler, H.R., Bartussek, C., Eckwert, H., Farian, K., Gränzer, S., Knigge, T., Kunz, N., 2001. The hepatic stress protein (hsp70) response to interacting abiotic parameters in fish exposed to various levels of pollution. J. Aquat. Ecosyst. Stress. Recover. 8 (3–4), 261–279. Lacroix, C., Coquillé, V., Guyomarch, J., Auffret, M., Moraga, D., 2014. A selection of reference genes and early-warning mRNA biomarkers for environmental monitoring using Mytilus spp. as sentinel species. Mar. Pollut. Bull. 86 (1–2), 304–313. Langmead, B., Salzberg, S.L., 2012. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9 (4), 357. Lech, J.J., Bendt, J.R., 1980. Relationship between biotransformation and the toxicity and fate of xenobiotic chemicals in fish. Environmental Health Perspective 34, 115–131. Leung, P.T., Ip, J.C., Mak, S.S., Qiu, J.W., Lam, P.K., Wong, C.K., Chan, L.L., Leung, K.M., 2014. De novo transcriptome analysis of Perna viridis highlights tissue-specific patterns for environmental studies. BMC Genomics 15 (1), 804. Li, B., Dewey, C.N., 2011. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC bioinformatics 12 (1), 323. Li, H., Liu, B., Huang, G., Fan, S., Zhang, B., Su, J., Yu, D., 2017. Characterization of transcriptome and identification of biomineralization genes in winged pearl oyster (Pteria penguin) mantle tissue. Comparative Biochemistry and Physiology Part D: Genomics and Proteomics 21, 67–76. Livingstone, D.R., 1998. The fate of organic xenobiotics in aquatic ecosystems: quantitative and qualitative differences in biotransformation by invertebrates and fish. Comparative Biochemistry and Physiology - A Molecular and Integrative Physiology 120 (1), 43–49. https://doi.org/10.1016/S1095-6433(98)10008-9. Lopes, B., Ferreira, A.M., Bebianno, M.J., 2012. Responses of CYP450 dependent system to aliphatic and aromatic hydrocarbons body burden in transplanted mussels from South coast of Portugal. Ecotoxicology 21 (3), 730–749. Lüchmann, K.H., Mattos, J.J., Siebert, M.N., Dorrington, T.S., Toledo-Silva, G., Stoco, P.H., Grisard, E.C., Bainy, A.C.D., 2012. Suppressive subtractive hybridization libraries prepared from the digestive gland of the oyster Crassostrea brasiliana exposed to a diesel fuel water-accommodated fraction. Environ. Toxicol. Chem. 31 (6), 1249–1253. Lüchmann, K.H., Dafre, A.L., Trevisan, R., Craft, J.A., Meng, X., Mattos, J.J., Zacchi, F.L., T.S., Dorrington, Schroeder, D.C., Bainy, A.C.D., 2014. A light in the darkness: New biotransformation genes, antioxidant parameters and tissue-specific responses in oysters exposed to phenanthrene. Aquat. Toxicol. 152, 324–334. https://doi.org/10. 1016/j.aquatox.2014.04.021. Martínez-Gómez, C., Vethaak, A.D., Hylland, K., Burgeot, T., Köhler, A., Lyons, B.P., Thain, J., Gubbins, M.J., Davies, I.M., 2010. A guide to toxicity assessment and monitoring effects at lower levels of biological organization following marine oil spills in European waters. ICES J. Mar. Sci. 67 (6), 1105–1118. McDowell, I.C., Nikapitiya, C., Aguiar, D., Lane, C.E., Istrail, S., Gomez-Chiarri, M., 2014. Transcriptome of American oysters, Crassostrea virginica, in response to bacterial challenge: insights into potential mechanisms of disease resistance. PLoS One 9 (8), e105097. Meador, J.P., 2003. Bioaccumulation of PAH's in marine invertebrates. In: Douben, P.E. (Ed.), PAHs: An Ecotoxicological Perspective. John Wiley & Sons, Chichester, UK, pp. 147–171. Men, B., He, M., Tan, L., Lin, C., Quan, X., 2009. Distributions of polycyclic aromatic hydrocarbons in the Daliao River Estuary of Liaodong Bay, Bohai Sea (China). Mar. Pollut. Bull. 58, 818–826. Miao, J., Pan, L., Liu, N., Xu, C., Zhang, L. 2011. Molecular cloning of CYP4 and GSTpi homologues in the scallop Chlamys farreri and its expression in response to benzo[a] pyrene exposure. Marine genomics,4(2), 99-108. Murphy, B., Lingman, S., Richter, B., Cralson, R., 2012. Simultaneous extraction of PAHs and PCBs from environmental samples using accelerated solvent extraction. Application note 1025. Thermo Fisher Scientific. Available from: http://www. dionex.com/en-us/webdocs/113838-AN1025-ASE-PCBs-PAHs-mussel-tissue-soilAN70253_E.pdf. Neff, J. M. 1979. Polycyclic Aromatic Hydrocarbons in the Aquatic Environment: Sources, Fates and Biological Effects. Applied Science, Barking, Essex, UK, 262 pp. Neff, J.M., 2002. Bioaccumulation in Marine Organisms: Effect of Contaminants From Oil Well Produced Water (First). Elsevier Ltd., Amsterdam. Neff, J.M., Boehm, P.D., Haensly, W.E., 1985. Petroleum contamination and biochemicalalterations in oysters (Crassostrea gigas) and plaice (Pleuronectes platessa) from bays impacted by the Amoco-Cadiz crude-oil spill. Mar. Environ. Res. 17 (2–4), 281–283. Neff, J.M., Ostazeski, S., Gardiner, W., Stejskal, I., 2000. Effects of weathering on the toxicity of three offshore Australian crude oils and a diesel fuel to marine animals. Environ. Toxicol. Chem. 19 (7), 1809–1821. Newell, R.I.E., Fisher, T.R., Holyoke, R.R., Cornwell, J.C., 2005. Influence of eastern oysters on nitrogen and phosphorus regeneration in Chesapeake Bay, USA. In: Dame, R.F., Olenin, S. (Eds.), The Comparative Roles of Suspension-feeders in Ecosystems. NATO Science Series IV: Earth and Environmental Series, Vol, vol. 47. Springer, Dordrecht, pp. 93–120.

79–89. Droege, M., Hill, B., 2008. The genome sequencer FLX system-longer reads, more applications, straight forward bioinformatics and more complete data sets. J. Biotechnol. 136, 3–10. Eisler, R., 1987. Polycyclic aromatic hydrocarbon hazards to fish, wildlife, and invertebrates: a synoptic review. U.S. Fish and Wildlife Service Biological Report 85 (1.11), 1–55. Enwere, R., 2009. Environmental Risk Management of Contamination of Marine Biota by Hydrocarbons Specifically Those Arising Following an Oil Spill. Ph.D. Thesis. The Robert Gordon University (257 pp). Epelboin, Y., Quintric, L., Guévélou, E., Boudry, P., Pichereau, V., Corporeau, C., 2016. The Kinome of Pacific oyster Crassostrea gigas, its expression during development and in response to environmental factors. PLoS One 11 (5), e0155435. Ertl, N.G., O'Connor, W.A., Brooks, P., Keats, M., Elizur, A., 2016. Combined exposure to pyrene and fluoranthene and their molecular effects on the Sydney rock oyster, Saccostrea glomerata. Aquat. Toxicol. 177, 136–145. Finch, B.E., Stefansson, E.S., Langdon, C.J., Pargee, S.M., Blunt, S.M., Gage, S.J., Stubblefield, W.A., 2016. Photo-enhanced toxicity of two weathered Macondo crude oils to early life stages of the eastern oyster (Crassostrea virginica). Mar. Pollut. Bull. 113 (1–2), 316–323. Fingas, M., 2001. The Basics of Oil Spill Cleanup. Lewis Publishershttps://doi.org/10. 1016/0025-326X(80)90262-3. 245 pp. Freire, J., Labarta, U., 2003. The Prestige: impacts on resources and marine ecosystems. In: Fernández-Latorre, S. (Ed.), The footprint of the Fuel. Essays on Prestige. Fundación, Coruña, pp. 104–135 (In Spanish). Frouin, H., Pellerin, J., Fournier, M., Pelletier, E., Richard, P., Pichaud, N., Rouleau, C., Garnerot, F., 2007. Physiological effects of polycyclic aromatic hydrocarbons on softshell clam Mya arenaria. Aquat. Toxicol. 82 (2), 120–134. Garret, R.H., Grisham, C.M., 2017. Biochemistry, Sixth edition. Cengage Learning, Boston, MA, USA 1330 pp. Gauthier, P.T., Norwood, W.P., Prepas, E.E., Pyle, G.G., 2014. Metal–PAH mixtures in the aquatic environment: a review of co-toxic mechanisms leading to more-than-additive outcomes. Aquat. Toxicol. 154, 253–269. Ghiselli, F., Milani, L., Chang, P.L., Hedgecok, D., Davis, J.P., Nuzhdin, S.V., Passamonti, M., 2012. De novo assembly of the Manila clam Ruditapes philippinarum transcriptome provides new insights into expression bias, mitochondrial doubly uniparental inheritance and sex determination. Molecular Biology Evolution 29 (2), 771–786. Glaser, J.A., Foerst, D.L., McKee, G.D., Quave, S.A., Budde, W.L., 1981. Trace analyses for wastewaters. Environ. Sci. Technol. 15, 1426–1435. Grabherr, M.G., Haas, B.J., Yassour, M., Levin, J.Z., Thompson, D.A., Amit, I., Adiconis, X., Fan, L., Raychowdhury, R., Zeng, Q., Chen, Z., Mauceli, E., Lindblad-Toh, K., Friedman, N., Regev, A., 2011. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29 (7), 644–652. Grundy, M.M., Ratcliffe, N.A., Moore, M.N., 1996. Immune inhibition in marine mussels by polycyclic aromatic hydrocarbons. Mar. Environ. Res. 42 (1–4), 187–190. Haas, B.J., Papanicolaou, A., Yassour, M., Grabherr, M., Blood, P.D., Bowden, J., Couger, M.B., Eccles, D., Li, B., Lieber, M., MacManes, M.D., Ott, M., Orvis, J., Pochet, N., Strozzi, F., Weeks, N., Westerman, R., William, T., Dewey, C.N., Henschel, R., LeDuc, R.D., Friedman, N., Regev, A., 2013. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat. Protoc. 8 (8), 1494. Hampel, M., Blasco, J., Díaz, M.L.M., 2016. Biomarkers and effects. In: Blasco, J., Chapman, P.M., Campana, O., Hampel, M. (Eds.), Marine Ecotoxicology: Current Knowledge and Future Issues. Academic Press 321 pp. Hu, B., Zhao, J., Lai, B., Qin, Y., Wang, H., Hu, G., 2016. LcGST4 is an anthocyaninrelated glutathione S-transferase gene in Litchi chinensis Sonn. Plant Cell Rep. 35 (4), 831–843. Jenny, M.J., Walton, W.C., Payton, S.L., Powers, J.M., Findlay, R.H., O'Shields, B., Diggins, B., Pinkerton, M., Porter, D., Crane, D.M., Tapley, J., Cunningham, C., 2016. Transcriptomic evaluation of the American oyster, Crassostrea virginica, deployed during the Deepwater Horizon oil spill: evidence of an active hydrocarbon response pathway. Mar. Environ. Res. 120, 166–181. Jeong, W.G., Cho, S.M., 2007. Long-term effect of polycyclic aromatic hydrocarbon on physiological metabolisms of the Pacific oyster, Crassostrea gigas. Aquaculture 265 (1–4), 343–350. Jiang, X., Qiu, L., Zhao, H., Song, Q., Zhou, H., Han, Q., Diao, X., 2016. Transcriptomic responses of Perna viridis embryo to benzo(a)pyrene exposure elucidated by RNA sequencing. Chemosphere 163, 125–132. Jing-Jing, M., Lu-Qing, P., Jing, L., Lin, Z., 2009. Effects of benzo [a] pyrene on DNA damage and histological alterations in gonad of scallop Chlamys farreri. Mar. Environ. Res. 67, 47–52. Jørgensen, A., Giessing, A., Rasmussen, L.J., Andersen, O., 2005. Biotransformation of the polycyclic aromatic hydrocarbon pyrene in the marine polychaete Nereis virens. Environ. Toxicol. Chem. 24 (11), 2796–2805. Karatolos, N., Pauchet, Y., Wilkinson, P., Chauhan, R., Denholm, I., Gorman, K., Nelson, D.R., Bass, C., French-Constant, R.H., Williamson, M.S., 2011. Pyrosequencing the transcriptome of the greenhouse whitefly, Trialeurodes vaporariorum reveals multiple transcripts encoding insecticide targets and detoxifying enzymes. BMC Genomics 12 (1), 56. Karnauskas, M., Schirripa, M.J., Kelble, C.R., Cook, G.S., Craig, J.K., 2013. Ecosystem status report for the Gulf of Mexico. NOAA Technical Memorandum NMFS-SEFSC 653, 52. Kendziorski, C.M., Zhang, Y., Lan, H., Attie, A.D., 2003. The efficiency of pooling mRNA in microarray experiments. Biostatistics 4 (3), 465–477. Kendziorski, C., Irizarry, R.A., Chen, K.S., Haag, J.D., Gould, M.N., 2005. On the utility ofpooling biological samples in microarray experiments. Proc. Natl. Acad. Sci. 102

15

Comparative Biochemistry and Physiology, Part C 225 (2019) 108571

E.A. López-Landavery, et al.

Park, MD, pp. 467–513. Sindermann, C. J. 2006. Coastal Pollution: Effects on Living Resources and Humans. CRC Press, Boca Raton, Florida, USA. 306 pp. Sistema de Venecia, 1958. Symposium on the classification of brackish waters, Venice April 8–14, 1958. Archives Oceanography and Limnology 11 (Suppl), 1–248. Snyder, M.J., 2000. Cytochrome P450 enzymes in aquatic invertebrates: recent advances and future directions. Aquat. Toxicol. 48 (4), 529–547. Srivastava, P., 2002. Roles of heat-shock proteins in innate and adaptive immunity. Nat. Rev. Immunol. 2 (3), 185. Tian, S., Pan, L., Sun, X., 2013. An investigation of endocrine disrupting effects and toxic mechanisms modulated by benzo [a] pyrene in female scallop Chlamys farreri. Aquat. Toxicol. 144, 162–171. Tolosa, I., de Mora, S.J., Fowler, S.W., Villeneuve, J., Bartocci, J., Cattini, C., 2005. Aliphatic and aromatic hydrocarbons in marine biota and coastal sediments from the Gulf and the Gulf of Oman. Mar. Pollut. Bull. 50, 1619–1633. Tymoczko, J. L., Berg, J. M., Stryer, L. 2013. Biochemistry: A Short Course. Second edition. Macmillan Education, New York, NY, USA. 862 pp. Valavanidis, A., Vlahogianni, T., Dassenakis, M., Scoullos, M., 2006. Molecular biomarkers of oxidative stress in aquatic organisms in relation to toxic environmental pollutants. Ecotoxicol. Environ. Saf. 64 (2), 178–189. Van der Oost, R., Beyer, J., Vermeulen, N.P.E., 2003. Fish bioaccumulation and biomarkers in environmental risk assessment: a review. Environ. Toxicol. Pharmacol. 13, 57–149. Vargas, J.D., Herpers, B., McKie, A.T., Gledhill, S., McDonnell, J., Van Den Heuvel, M., Davies, K.E., Ponting, C.P., 2003. Stromal cell-derived receptor 2 and cytochrome b561 are functional ferric reductases. Biochimica et Biophysica Acta (BBA)-Proteins and Proteomics 1651 (1–2), 116–123. Varney, R.L., Galindo-Sánchez, C.E., Cruz, P., Gaffney, P.M., 2009. Population genetics of the eastern oyster Crassostrea virginica (Gmelin, 1791) in the Gulf of Mexico. J. Shellfish Res. 28 (4), 855–864. Vignet, C., Larcher, T., Davail, B., Joassard, L., Le Menach, K., Guionnet, T., Lyphout, L., Ledevin, M., Goubeau, M., Budzinski, H., Bégout, M.L., 2016. Fish reproduction is disrupted upon lifelong exposure to environmental PAHs fractions revealing different modes of action. Toxics 4 (4), E26. Wang, S., Diller, K.R., Aggarwal, S.J., 2003. Kinetics study of endogenous heat shock protein 70 expression. J. Biomech. Eng. 125 (6), 794–797. Wang, Z., Shao, Y., Li, C., Zhang, W., Duan, X., Zhao, X., Qiu, Q., Jin, C., 2016. RNA-seq analysis revealed ROS-mediated related genes involved in cadmium detoxification in the razor clam Sinonovacula constricta. Fish & shellfish immunology 57, 350–361. Wang, L., Song, X., Song, L., 2018. The oyster immunity. Developmental & Comparative Immunology 80, 99–118. Widdows, J., Bakke, T., Bayne, B.L., Donkin, P., Livingstone, D.R., Lowe, D.M., Moore, M.N., Evans, S.V., Moore, S.L., 1982. Responses of Mytilus edulis on exposure to the water-accommodated fraction of North Sea oil. Mar. Biol. 67 (1), 15–31. https://doi. org/10.1007/BF00397090. Wrigth, D. A., Welburn, P. 2002. Environmental Toxicology, vol. 11. Cambridge University Press, New York, USA. 630 pp. Xu, J., Sun, J., Chen, J., Wang, L., Li, A., Helm, M., Dubovsky, S.L., Bacanu, S.A., Zhao, Z., Chen, X., 2012. RNA-Seq analysis implicates dysregulation of the immune system in schizophrenia. BMC Genomics 13 (8), S2. Yang, C., Wang, L., Zhang, H., Wang, L., Huang, M., Sun, Z., Sun, Y., Song, L., 2014. A new fibrinogen-related protein from Argopecten irradians (AiFREP-2) with broad recognition spectrum and bacteria agglutination activity. Fish & shellfish immunology 38 (1), 221–229. Young, M.D., Wakefield, M.J., Smyth, G.K., Oshlack, A., 2010. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 11 (2), R14. Zacchi, F.L., de Lima, D., Flores-Nunes, F., Mattos, J.J., Lüchmann, K.H., de Miranda Gomes, C.H.A., Bícego, M.C., Taniguchi, S., Sasaki, S.T., Bainy, A.C.D., 2017. Transcriptional changes in oysters Crassostrea brasiliana exposed to phenanthrene at different salinities. Aquat. Toxicol. 183, 94–103. Zhang, X., Mao, Y., Huang, Z., Qu, M., Chen, J., Ding, S., Hong, J., Sun, T., 2012. Transcriptome analysis of the Octopus vulgaris central nervous system. BMC Genomics 7 (6), e40320. Zhang, G., Fang, X., Guo, X., Li, L., Luo, R., Xu, F., Yang, P., Zhang, L., Wang, X., Qi, H., Xiong, Z., Que, H., Xie, Y., Holland, P., Paps, J., Zhu, Y., Wu, F., Chen, Y., Wang, Y., Peng, C., Meng, J., Yang, L., Liu, J., Wen, B., Zhang, N., Huang, Z., Zhu, Q., Feng, Y., Mount, A., Hedgecok, D., Xu, Z., Liu, Y., Domazet-Lošo, T., Du, Y., Sun, X., Zhang, S., Liu, B., Cheng, P., Jiang, Z., Li, J., Fan, D., Wang, W., Fu, W., Wang, T., Wang, B., Zhang, J., Peng, Z., Li, Y., Li, N., Wang, J., Chen, M., He, Y., Tan, F., Song, X., Zheng, Q., Huang, R., Yang, H., Du, X., Chen, L., Yang, M., Gaffney, P.M., Wang, S., Luo, L., She, Z., Ming, Y., Huang, W., Zhang, S., Huang, B., Zhang, Y., Qu, T., Ni, P., Miao, G., Wang, J., Wang, Q., Steinberg, C.E.W., Wang, H., Li, N., Qian, L., Zhang, G., Li, Y., Yang, H., Liu, X., Wang, J., Yin, Y., Wang, J., 2012. The oyster genome reveals stress adaptation and complexity of shell formation. Nature 490 (7418), 49. Zhang, L., Li, L., Zhu, Y., Zhang, G., Guo, X., 2014. Transcriptome analysis reveals a rich gene set related to innate immunity in the eastern oyster (Crassostrea virginica). Mar. Biotechnol. 16 (1), 17–33.

Nogueira, D.J., Mattos, J.J., Dybas, P.R., Flores-Nunes, F., Sasaki, S.T., Taniguchi, S., Schmidt, E.C., Bouzon, Z.L., Bícego, M.C., Melo, C.M.R., Toledo-Silva, G., Bainy, A.C.D., 2017. Effects of phenanthrene on early development of the Pacific oyster Crassostrea gigas (Thunberg, 1789). Aquat. Toxicol. 191, 50–61. Noh, E.M., Park, Y.J., Kim, J.M., Kim, M.S., Kim, H.R., Song, H.K., Hong, O.Y., So, H.S., Yang, S.H., Kim, J.S., Park, S.H., Youn, H.J., You, Y.O., Choi, K.B., Kwon, K.B., Lee, Y.R., 2015. Fisetin regulates TPA-induced breast cell invasion by suppressing matrix metalloproteinase-9 activation via the PKC/ROS/MAPK pathways. Eur. J. Pharmacol. 764, 79–86. Omar-Ali, A., Hohn, C., Allen, P.J., Rodriguez, J., Petrie-Hanson, L., 2015. Tissue PAH, blood cell and tissue changes following exposure to water accommodated fractions of crude oil in alligator gar, Atractosteus spatula. Mar. Environ. Res. 108, 33–44. Parcellier, A., Gurbuxani, S., Schmitt, E., Solary, E., Garrido, C., 2003. Heat shock proteins, cellular chaperones that modulate mitochondrial cell death pathways. Biochem. Biophys. Res. Commun. 304 (3), 505–512. Patnaik, B.B., Wang, T.H., Kang, S.W., Hwang, H.J., Park, S.Y., Park, E.B., Chung, J.M., Son, D.K., Kim, C., Kim, S., Lee, J.S., Han, Y.S., Park, H.S., Lee, J.S., 2016. Sequencing, de novo assembly, and annotation of the transcriptome of the endangered freshwater pearl bivalve, Cristaria plicata, provides novel insights into functional genes and marker discovery. PLoS One 11 (2), e0148622. Petrulis, J.R., Chen, G., Benn, S., LaMarre, J., Bunce, N.J., 2001. Application of the ethoxyresorufin-O-deethylase (EROD) assay to mixtures of halogenated aromatic compounds. Environ. Toxicol. 16 (2), 177–184. Piazza, R.S., Trevisan, R., Flores-Nunes, F., Toledo-Silva, G., Wendt, N., Mattos, J.J., Lima, D., Taniguchi, S., Sasaki, T.S., Mello, A.C.P., Zacchi, F.L., Serrano, M.A.S., Gomes, H.A.M., Bícego, M.C., de Almeida, E.A., Bainy, A.C.D., 2016. Exposure to phenanthrene and depuration: changes on gene transcription, enzymatic activity and lipid peroxidation in gill of scallops Nodipecten nodosus. Aquat. Toxicol. 177, 146–155. Prada-Ríos, J.E., Zambrano-Ortiz, M.M., 2006. Acute toxicity and bioaccumulation of two polycyclic aromatic hydrocarbons (naphthalene and fluoranthene) in a bivalve mollusc. Boletín Científico CCCP 13, 53–64 (In Spanish). Rajkumar, A.P., Qvist, P., Lazarus, R., Lescai, F., Ju, J., Nyegaard, M., Mors, O., Børglum, A.D., Li, Q., Christensen, J.H., 2015. Experimental validation of methods for differential gene expression analysis and sample pooling in RNA-Seq. BMC Genomics 16 (1), 548. Rewitz, K.F., Kjellerup, C., Jørgensen, A., Petersen, C., Andersen, O., 2004. Identification of two Nereis virens (Annelida: Polychaeta) cytochromes P450 and induction by xenobiotics. Comparative Biochemistry and Physiology Part C: Toxicology & Pharmacology 138 (1), 89–96. Rewitz, K.F., Styrishave, B., Løbner-Olesen, A., Andersen, O., 2006. Marine invertebrate cytochrome P450: emerging insights from vertebrate and insect analogies. Comparative Biochemistry and Physiology Part C: Toxicology & Pharmacology 143 (4), 363–381. Reynaud, S., Deschaux, P., 2006. The effects of polycyclic aromatic hydrocarbons on the immune system of fish: a review. Aquat. Toxicol. 77 (2), 229–238. https://doi.org/ 10.1016/j.aquatox.2005.10.018. Riches, Z., Stanley, E.L., Bloomer, J.C., Coughtrie, M.W., 2009. Quantitative evaluation of the expression and activity of five major sulfotransferases (SULTs) in human tissues: the SULT “pie”. Drug Metab. Dispos. 37 (11), 2255–2261. Robinson, M.D., McCarthy, D.J., Smyth, G.K., 2010. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26 (1), 139–140. Roméo, M., Giambérini, L., 2013. History of biomarkers. In: Amiard-Triquet, C., Amiard, J.C., Rainbow, P.S. (Eds.), Ecological Biomarkers: Indicators of Ecotoxicological Effects. CRC Press, Boca Raton, FL, USA, pp. 15–44. Roméo, M., Wirgin, I., 2011. Biotransformation of organics contaminants and the acquisition of resistance. In: Amiard-Triquet, C., Rainbow, P.S., Roméo, M. (Eds.), Tolerance to environmental contaminants. CRC Press, Boca Raton, FL, USA, pp. 175–208. SAGARPA, 2013. Statistical Yearbook of Aquaculture and Fisheries (Anuario estadístico de acuacultura y pesca). 385 pp. (In Spanish). Sambrook, J., Maccallum, P., Russel, D. 2001. Molecular Cloning: A Laboratory Manual, Third ed. Cold Springs Harbour Press, NY, USA. 2222 pp. Schäfer, S., Köhler, A., 2009. Gonadal lesions of female sea urchin (Psammechinus miliaris) after exposure to the polycyclic aromatic hydrocarbon phenanthrene. Mar. Environ. Res. 68 (3), 128–136. https://doi.org/10.1016/j.marenvres.2009.05.001. Schuler, G.D., 2001. Sequence alignment and database searching. In: Baxevanis, A.D., Ouellette, B.F. (Eds.), Bioinformatics: A Practical Guide to the Analysis of Genes and Proteins. vol. 43. John Wiley & Sons, pp. 187–214. Schurch, N.J., Schofield, P., Gierliński, M., Cole, C., Sherstnev, A., Singh, V., Wrobel, N., Gharbi, K., Simpson, G.G., Owen-Hughes, T., Blaxter, M., Barton, G., 2016. How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use? RNA 22 (6), 839–851. Sericano, J.L., Wade, T.L., Brooks, J.M., 1996. Accumulation and depuration of organic contaminants by the American oyster (Crassostrea virginica). Sci. Total Environ. 179, 149–160. Shumway, S., 1996. Natural environmental factors. In: Kennedy, V., Newell, R., Eble, A. (Eds.), The Eastern Oyster Crassostrea virginica. Maryland Sea Grant College, College

16