Comparative analysis of distinctive transcriptome profiles with biochemical evidence in bisphenol S- and benzo[a]pyrene-exposed liver tissues of the olive flounder Paralichthys olivaceus

Flounder is a promising model species for environmental monitoring of coastal regions. To assess the usefulness of liver transcriptome profiling, juvenile olive flounder Paralichthys olivaceus were exposed to two pollutants, bisphenol S (BPS) and benzo[a]pyrene (BaP), which have different chemical characteristics and have distinct modes of metabolic action in teleost. Six hours after intraperitoneal injection with BPS (50 mg/kg bw) or BaP (20 mg/kg bw), liver transcriptomes were analyzed using the Illumina Hiseq 3000 platform. Interestingly, the transcriptome was highly sensitive and was distinctively expressed in response to each chemical. The primary effect of BPS was significantly increased transcription of egg process and vitellogenesis related genes, including vitellogenins (vtg1, vtg2), zona pellucida sperm-binding proteins (zp3, zp4), and estrogen receptors (erα, erβ), with increases in plasma 17β-estradiol (E2) and vitellogenin (VTG) concentrations. Following BaP treatment, detoxification- and biotransformation-related genes such as cyp1a1 and UDP-glucuronosyltransferase (ugt1a1) were significantly increased, with an increase in EROD activity. In both transcriptomes, mRNA expression of genes involved in antioxidant defense systems was increased, while genes involved in innate immunity were decreased upon BPS or BaP exposure with a decrease in complement activity. This study provides useful insight into the chemical-specific hepatic transcriptional response of P. olivaceus and suggests a basis for further studies examining biomarker application of liver transcriptomes for environmental pollution.


Introduction
Flatfish are bottom-dwelling fishes that spend most of their life cycle on the bottom of estuaries and coastal regions. Thus, flatfish are acutely or chronically exposed to a wide range of environmental pollutants derived from agricultural and industrial wastes, municipal sewage, and human activity [1]. Due to their benthic habitats and wide-ranging geographical distributions, flatfish could be potential indicators of water and sediment quality, using their molecular, biochemical, and physiological sensitivity. Numerous biomarkers have been continuously developed to understand the potential effects of environmental pollutants on flatfish in their role as sentinels [2][3][4][5]. The capacity of a biomarker to respond may be prolonged or restricted by complex molecular pathways with combinations of certain signal cascades, as numerous aquatic pollutants are released in different forms with distinct modes of action in waterbodies. Particularly, biochemical markers, including molecular responses, have been identified as powerful and cost-effective approaches to obtain information on the status of the environment and the effects of pollution on resident species, including flatfish [5][6][7][8][9]. However, the accessibility of genomic resources (e.g., whole genome and transcriptome) is limited in most marine fish, even though gene expression profiling has been widely applied to predict the potential toxicity of various chemical compounds and to elucidate the underlying molecular mechanisms.
The olive flounder Paralichthys olivaceus (also known as the Japanese flounder or bastard halibut) used in this study is one of the most cosmopolitan aquaculture species in east and south Asia, including South Korea, Japan, and China. P. olivaceus is desirable for commercial culture and was subsequently introduced in many Asia countries for aquaculture production. In 2011, the catch of wild P. olivaceus was approximately 4,600 metric tons and its aquaculture production was over 40,000 metric tons [10]. One of the most important steps in developing a model species for environmental monitoring is the availability of genomic information. Genomic platforms have been successfully applied to understand diverse molecular and physiological characteristics of P. olivaceus [11][12][13][14]. The recently published whole genome sequence of P. olivaceus provides a unique advantage for genomic application among flatfish [15]. Despite these advantages, application of genomic information of P. olivaceus in marine environmental monitoring has received little attention to date.
The primary objective of the present study was to test whether the transcriptome profiling of the olive flounder could be usefully applied for the monitoring of marine pollutants. Because the sensitivity of signal pathways is very different for chemicals, two distinct chemical pollutants were employed in this study. Bisphenol S (BPS) is a structural analog of bisphenol A (BPA) and is used as an alternative for BPA due to increasing concerns over the endocrine disrupting potential of BPA [16]. However, many in vivo and in vitro studies have suggested that BPS is as hormonally active as BPA with endocrine disrupting effects [17]. Benzo[a]pyrene (BaP) is a polycyclic aromatic hydrocarbon (PAH) and a ubiquitous marine pollutant derived from numerous human activities [18]. Adverse effects of BPA (e.g., carcinogenicity, genotoxicity, immunotoxicity, and neurotoxicity) and its biotransformation mechanism have been extensively reviewed [19,20]. Biotransformation of xenobiotics and drug metabolism (e.g., oxidation, reduction, hydrolysis) occurs by the action of one or more enzymes in the liver tissue [21,22]. In particular, both chemicals have consistently been found up to microgram level in the sediments and waterbodies of the coastal regions of south Asia, which are known to be among the olive flounders' habitats [23][24][25][26]. Therefore, we comparatively analyzed BPS-and BaP-exposed liver tissues of P. olivaceus to consider its hepatic transcriptome as a promising tool for pollution monitoring of coastal regions.

Total RNA extraction and library construction
The entire procedure for Illumina RNA sequencing (RNA seq) was performed in Macrogen Inc. (Seoul, South Korea). A set of RNA samples for each chemical, comprised of three biological replicates, was not pooled. Each liver tissue was frozen in liquid nitrogen and homogenized with a glass pestle. Total RNA was extracted from individual olive flounder (two fish from control and three fish from BPS and BaP, respectively) using TRIzol TM Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. DNA digestion was performed using DNase I (Sigma Aldrich, St Louis, MO, USA). Prior to Illumina RNA sequencing, total RNA quantity and quality were checked by analyzing the ratios A230/260 and A260/280 using a NanoDrop 1 2000 Spectrophotometer (Thermo Scientific, Wilmington, DE, USA) and an Agilant 2100 Bioanalyzer (Agilent, Böblingen, Germany) respectively. Samples satisfying the criteria of RNA integrity > 7.5 were used for library preparation.
Library construction was performed using Tureseq™ RNA sample prep Kit (Illumina, San Diego, CA, USA) at Macrogen Inc. (Seoul, South Korea) [27]. Briefly, the poly A-containing mRNA was purified using oligo dT-coated magnetic beads. The first strand of cDNA was synthesized using random hexamers and subsequently double-stranded cDNA was synthesized and purified [28]. After adaptor ligation, proper DNA fragments were selected and enriched to create the final cDNA library template. All libraries were sequenced by Illumina Hiseq 3000 with read length of 101 bp.

De novo assembly, differentially expressed genes (DEGs) identification, and gene ontology (GO) analysis
(nr) protein database to evaluate sequence similarity to genes of other species using an E-value cutoff of 1e-06.
The abundance of each transcript was analyzed by direct alignment of each read to the control group assembly using the Bowtie 2 alignment tool (http://bowtie-bio.sourceforge.net/ bowtie2) with the default parameters and RSEM method (RNA-Seq by Expectation-Maximization; http://deweylab.biostat.wisc.edu/rsem). The quantified number of mapped reads was normalized as fragments per kilobase of transcripts per million fragments mapped (FPKM) value. The false discovery rate (FDR) was used to predict the P value threshold for statistical analysis. The FPKM value of each sample was compared using the R package provided by Trinity platform. When the FDR value of a contig was less than 0.05 and its FPKM value was more than a two-fold change, the contig was considered to contain a differentially expressed gene (DEG). Hierarchical clustering analysis was employed in order to construct a heat map using MeV v.7.4 (Dana-Farber Cancer Institute, Boston, MA, USA) software. Gene Ontology (GO) analysis of the contigs was performed using the GOstats program (http://www. bioconductor.org/packages/3.3/bioc/html/GOstats) and Fisher's Exact Test (P < 0.05), as implemented in the sequence annotation tool Blast2GO [30].

Data deposition
All the raw data are available in the Sequencing Read Archive (SRA) of NCBI under the Bio-Project number PRJNA421418, with the accession number SRP126397 (SRX3459054-SRX3459061).

Quantitative PCR validation
Total RNAs were reverse-transcribed into cDNA using a first-strand cDNA synthesis kit (Invitrogen, Carlsbad, CA, USA). The optical density at 260 nm and 280 nm (OD260/280) was approximately measured as 1.9, and the OD260/230 was between 1.8 and 1.9. The quantitative reverse transcription polymerase chain reaction (qRT-PCR) was performed using a two-step procedure. The β-actin gene was used as a reference control for real-time qRT-PCR. To quantify the mRNA expression level, we used the comparative CT method (2 -ΔΔCt method) with Roto-Gene Q (Qiagen, Hilden, Germany) according to the manufacturer's instruction. mRNA expression for each gene was measured in triplicate. The specific primers for the qRT-PCR analysis were designed based on RNA-SEQ library using the ABI PRISM Primer Express software (Applied Biosystems, Darmstadt, Germany).

EROD activity
The 7-ethoxyresorufin-O-deethylase (EROD) activity assay was conducted according to previous studies with slight modifications (e.g. buffer volume, employed basic instruments) [31]. The liver tissues of individual fish were homogenized in 0.25 M sucrose and centrifuged at 9,000 × g for 20 min at 4˚C. To sample the microsomal fraction, the supernatant was transferred and centrifuged at 105,000 × g for 60 min at 4˚C. Approximately 35 μg of microsomal fraction was transferred into the EROD buffer (0.1 M NaPO4, pH 7.6) with the addition of 7-ethoxyresorufin (ER) solution (0.4 mM ER in DMSO) and NADPH solution (10 mM NADPH in distilled water). The background subtraction was prepared with the same mixed buffer with the absence of the microsomal fraction. The resourufin production was measured with 530 nm excitation and 595 nm emission filters with a Varioskan Flash fluorometer (Thermo Fisher Scientific, Tewksbury, MA, USA). Results were represented as μmol min/mg protein.

Plasma analysis
In our preliminary study, no significant changes in plasma 17β-estradiol (E2) or vitellogenin (VTG) were observed at 6 h after the injections of both chemicals. Thus, we extended our analysis to 24 h after the injection.
The E2 levels in the plasma were quantified by a standard ethyl ether extraction method using a commercial Estradiol ELISA kit (Cayman Chemical Company, Michigan, USA). Approximately 1 ml of blood was taken from the caudal vein of individual anesthetized olive flounders (n = 5) with a heparin syringe. The whole blood sample was centrifuged at 3,000 × g for 5 min and the supernatant was transferred for E2 analysis. The extract was centrifuged at 2,000 × g for 30 min and the ether phase was removed. The samples were dried using a centrifugal evaporator. The intra-assay coefficients of variance for both the fish were < 10%. The E2 concentration was measured with a Varioskan Flash spectrophotometer (Thermo Fisher Scientific, Tewksbury, MA, USA) at 420 nm, and was expressed as pg/mL.
The concentration of VTG in the plasma samples was determined using a commercially available enzyme-linked immunosorbent assay (ELISA) kit for VTG (MyBioSource Inc, San Diego, USA), following the manufacturer's instructions. The intra-assay coefficients of variance for both the fish were < 12%. The VTG concentration was measured with a Varioskan Flash spectrophotometer (Thermo Fisher Scientific, Tewksbury, MA, USA) at 420 nm, and was expressed as μg/mL.

Haemolytic complement activity
The activity of the alternative complement pathway hemolytic activity (ACH) was conducted according to the classical method established in fish, with slight modification (e.g. buffer volume, employed basic instruments) [32]. Since no significant change in complement activity was observed at 6 h after the injections of both chemicals, we extended our analysis to 24 h after the injection. Sheep red blood cells (SRBC; 1.5 × 10 6 cells; National Institute of Toxicological Research, South Korea) were used as the targets of the assay. Reconstitution of the haemolytic activity of olive flounder serum was measured by incubation of the SRBC (5 μl) in 6% olive flounder serum (25 μl) in 10 mM phenol red-free Hank's buffer with 5 mM Mg 2+ and 0.15 NaCl (pH 7.3) in a 96-well plate. The plate was incubated for 90 min at 20˚C with gentle shaking. Haemolysis was measured spectrophotometrically with a Varioskan Flash spectrophotometer (Thermo Fisher Scientific, Tewksbury, MA, USA) at 414 nm. Complete (100%) and no haemolysis (0%) were measured by adding the washed SRBC (25 μl) to distilled water (100 μl) and phenol red-free Hank's buffer, respectively. The ACH was calculated as the reciprocal of the serum dilution, causing 50% lysis of SRBC (ACH 50 ; U/ml) based on the value of Y/ 1 − Y against the reciprocal of the serum dilutions on a log-log scaled graph.

Statistical analysis
Data are presented as mean ± standard deviation (S.D.). We performed a one-way analysis of variance (ANOVA) followed by Duncan's test to identify differences between the experimental groups. All analyses were performed using SPSS ver. 17.0 (SPSS Inc., Chicago, IL).

Sequencing and assembly
To compare the molecular response to different chemical characteristics in the olive flounder, we conducted an RNA-seq based investigation of transcriptional changes in liver tissues exposed to BPS or BaP for 6 h. Three libraries for each chemical and two libraries for the control group were subjected to Illumina RNA seq. The sequenced liver cDNA libraries contained a large number of raw reads, at 62 to 83 million reads per library (S1 Table). After the trimming process, 62-82 million reads were retained from each library.
The de novo assembly generated 124,125 transcripts for the liver tissue of P. olivaceus ( Table 1). The lengths of total transcripts ranged from 201 to 17,656 bp with an average length of 984 bp. The overall GC ratio of total transcripts was 47% and the N50 value of those transcripts was 2,110 bp. By employing a bioinformatics platform, 38,909 unigenes were filtered out from the total transcripts (S2 Table). Gene annotation of whole unigenes was performed by BLASTx analysis with the reciprocal BLAST best-hit method (E-value < 1 × 10-10) using NCBI nr protein database. Of unigenes, 28,919 (74%) and 29,979 (77%) transcripts were annotated with GO and KO identity, respectively.

mRNA expression profiling
Each library of the three groups (i.e. control, BPS exposure, BaP exposure) was individually aligned to the assembled reference transcriptome in order to compare differentially expressed transcripts. Overall mRNA expressions analyzed by both hierarchical clustering and multidimensional scaling (MDS) plot revealed that those of all three groups separated clearly (Fig 1,  S1 Fig). No outlier library was observed in each group, and relatively similar expression profiles were found with hierarchical clustering (Fig 1A). To investigate whether clear differences existed between the control and treated groups, we performed a principle component analysis (PCA), to provide a general overview ( Fig 1B). All members involved in each group closely clustered together. While the two control libraries had a relatively higher data variance regarding component two, the scatter was very low regarding the more indicative component one (39.9%). Taken together, both analyses revealed similar relative expression levels for all members analyzed in each group.
A total of 917 unique genes showed significantly different mRNA expression between control and BPS-exposed fish with the criteria > ± 2-fold change and P < 0.05 (603 transcripts upregulated; 314 transcripts downregulated), while 1,628 transcripts showed significant differences in the BaP-exposed fish compared to the control (993 transcripts upregulated; 635 transcripts downregulated) ( Fig 1C).
Of upregulated transcripts, 305 mRNAs of BPS-and 695 mRNAs of BaP-exposed fish were chemical-preferentially expressed, whereas 298 mRNAs were commonly upregulated in both chemical exposures ( Fig 1D). In the case of downregulated transcripts, 219 mRNAs were common to both chemicals, while 95 and 416 mRNAs were unique to BPS and BaP, respectively.

Fig 1. Comparison of transcriptional expression patterns of whole libraries and analysis of chemical-specifically expressed genes. (A)
Transcriptional pattern analysis of each library (i.e., two control, three BPS-exposed, and three BaP-exposed liver tissues) by employing heat map and hierarchical clustering; (B) PCA plot analysis of transcriptional profile of each library. Each sample is depicted with a different color. (C) Number of statistically significant transcripts (i.e., over 2 fold; P < 0.05) in P. olivaceus liver tissues exposed to BPS or BaP; (D) the number of uniquely or commonly up-or downregulated transcripts in the P. olivaceus liver tissues exposed to BPS or BaP. Detailed list of the commonly modulated genes is included in S5 and S6 Tables. https://doi.org/10.1371/journal.pone.0196425.g001 Chemical-specific transcriptome profiles of olive flounder

Analysis of gene ontology (GO) and KEGG pathway
Global functionality of the BPS-and BaP-exposed liver transcriptomes compared to the control group was analyzed by both GO analysis and KEGG pathway mapping. Transcripts were categorized according to their functions within each GO class, such as biological process, cellular process, or molecular function (Fig 2, S3 Table). Overall GO composition was analyzed and Comparison of Gene Ontology (GO) in terms of (A) biological process, (B) cellular components, and (C) molecular function that were enriched in the BPS-or BaP-exposed liver tissues of P. olivaceus. Composition of each GO term is represented as a percentage. For detailed information see in S3 Table. the majority was similar among transcriptomes. In the biological process class, the numbers with the GO term "Metabolic process" were increased in the BPS-(2,154) and BaP-exposed liver (2,297) compared that in the control (1,949). Similarly, counts for the GO term "Response to stimulus" were higher in the BPS-(461) and BaP-exposed liver (520) compared to the number in the control (401). Interestingly, the numbers of most GO terms involved in the molecular function class were higher in the BPS-and BaP-exposed liver transcriptomes than those of the control. However, the number of genes categorized in each function of the BPS-exposed liver transcriptome was higher than that from BaP for most GO terms of the three major categories.
To test the distinct sensitivity and further usefulness of hepatic transcriptome profiling between the BPS and BaP exposures, all the assembled reads were assigned on the functional classification of KEGG (Table 2 and S4 Table). In both transcriptomes, a high number of transcripts had annotations related to "Metabolic pathways" and "Biosynthesis of secondary metabolites". Overall KEGG classification showed similar distributions and compositions between transcriptomes. To analyze whether the gene composition matching certain pathways is chemical-specific, we selected a pathway which showed a distinct difference in gene number, and compared each gene's involvement. Although the pathway, entitled "Metabolic pathways", showed the largest matched gene number, we were not able to adopt it, as numerous pathways and sub-pathways are involved in the KEGG classification. Among the classifications, the "Protein processing in endoplasmic reticulum" pathway showed different matched genes, as 34 genes of BaP transcriptome were matched to the pathway, while only 12 genes from BPS transcriptome were incorporated in the pathway (Fig 3A).
Since the pathway of "Metabolism of benzo[a]pyrene by cytochrome P450" is registered in the KEGG database, we additionally checked the number and composition of matched genes to the pathway. As shown in Fig 3B, CYP1A-involved metabolism genes were highly expressed in the BaP-exposed liver transcriptome. In addition, significant increases in EROD activity were observed at 5 and 50 mg/kg bw of BaP injection into the liver tissue, while no significant change was measured in the BPS-exposed olive flounder (Fig 3C). The signaling, BaP-mediated CYP1A expression, is not uncommon in teleosts, but comparison of the mated gene clearly showed that the flounder hepatic transcriptome is differentially sensitive to specific chemicals; this profiling will be useful for monitoring marine contaminants.

Comparison of the BPS-and BaP-exposed liver transcriptomes
A comparative analysis of the P. olivaceus liver transcriptomes was performed to identify transcripts conserved in the liver tissue and those unique to each chemical exposure. Table 3 and  pyrene metabolism in "Metabolism of xenobiotics by cytochrome P450" in the BPS-or BaP-exposed liver tissues of P. olivaceus. The orange color denotes the gene mapped from the BPS-exposed liver transcriptome to each KEGG pathway. The blue-green color denotes the gene mapped from the BaP-exposed liver transcriptomes to each KEGG pathway. Notably, no gene was solely mapped from the BPS-exposed liver transcriptome to a KEGG pathway. (C) The hepatic EROD activity in the liver tissue of BPS-or BaP-injected olive flounder. Data are presented as the mean ± standard deviation (S.D.). Significant differences from the control value are indicated by an asterisk ( Ã ) on the data bar (P < 0.05). https://doi.org/10.1371/journal.pone.0196425.g003 Chemical-specific transcriptome profiles of olive flounder Table 4 show the 20 most highly up-and downregulated genes in the BPS-or BaP-exposed liver tissue, respectively. Overall, 4 of the 20 most highly upregulated genes in the BPS-exposed liver transcriptome were related to egg process and vitellogenesis, such as vitellogenin 1 (vtga), 2 (vtgb), and zona pellucida sperm-binding protein (zp3 and zp4) ( Table 3). Since the transcription of vtg genes and synthesis of the VTG protein are strongly associated with estrogen Chemical-specific transcriptome profiles of olive flounder level, we further analyzed estrogen-related transcriptional and physiological responses. In the increased genes, we identified that the estrogen receptor (ER) family such as erα and erβ genes were significantly increased compared to their expression in BaP-exposed liver tissue (Fig 4A).
To confirm the potential correlation between upregulation of the vtg and er genes, we measured the plasma E2 and VTG concentrations in the serum samples after BPS or BaP injection. Chemical-specific transcriptome profiles of olive flounder As a result, significant increases in plasma E2 were observed after 24 h in the BPS-injected olive flounder (P < 0.05), while no significant change was measured in the BaP treatment Chemical-specific transcriptome profiles of olive flounder (p > 0.05) (Fig 4B). A significant elevation in the plasma VTG was also observed only after 24 h in the BPS-injected olive flounder (P < 0.05) (Fig 4C).
In the case of the BaP-exposed liver, genes involved in drug/xenobiotic metabolism were highly upregulated, as the most highly expressed gene was cyp1a1, with 176 a fold change, followed by UDP-glucuronosyltransferase (ugt1a1) ( Table 4). We observed that 10 genes related to drug/xenobiotic metabolism and the antioxidant defense system (i.e. UDP-glucuronosyltransferase, ugt2A1; glutathione S-transferase, gst; microsomal GST, mgst; GST zeta, gstz; and catalase, cat) were commonly upregulated by BPS and BaP exposures (S5 Table). The hepatic enzymatic activity of the total GST protein was significantly elevated by the injections of both chemical (P < 0.05) (Fig 5A). Additionally, the CAT activity showed significant increases in its response to both chemicals (P < 0.05) (Fig 5B), while only BaP increased the SOD activity at 6 h after the injection (P < 0.05) (Fig 5C).
Although numerous genes associated with diverse pathways were downregulated by BPS or BaP exposure, mRNA expressions of several genes involved in the innate immunity including complement genes were decreased by both BPS and BaP exposure (Fig 4D, S6 Table). To further assess the potential correlation between the decrease in mRNA expression and protein activity, we measured the haemolytic complement activity of olive flounder serum. The plasma haemolytic complement activity of BPS-or BaP-exposed olive flounder was decreased, but only reached a statistically significant level after 24 h (P < 0.05) (Fig 4E). In this experiment, the DMSO also decreased the activity after 24 h (P < 0.05, observed in the plasma of BaPinjected olive flounder).

qPCR validation
To verify differentially expressed genes for each chemical, mRNA expressions of highly represented genes were analyzed by real-time qRT-PCR (Fig 6). The details on primer information are shown in S7 Table. Quantitative expressions of 15 randomly selected genes among the highly up-or downregulated genes in BPS-or BaP-exposed liver tissue showed the high reproducibility of the biological replicates. Most genes showed similar mRNA expression patterns between the two platforms, RNA-seq and qPCR, except for calr3 and ch in the BPS-injected sample, and vtg1 and zp4 in the BaP-injected sample.

Discussion
Overall comparison of GO compositions and their ratio among transcriptomes indicates that both chemicals have the potential to induce or inhibit transcriptional expression, and the modulated genes may be associated with cellular homeostasis in liver tissue. This result further implies that diverse pathways are active in the liver tissue of P. olivaceus, some of which have chemical-specific roles. KEGG analysis revealed that a high number of transcripts were associated with general metabolism and metabolite processes in both transcriptomes. A possible explanation for this could be that the fish liver tissue employs numerous enzymes and functional systems to maintain homeostasis by converting molecules (e.g., drug/xenobiotics metabolism) and regulating energy balance (e.g., glycogenesis, gluconeogenesis) [34]. Both chemicals used in this study are detoxified and excreted by metabolism in the fish liver [35][36][37][38][39]. BPS was biotransformed into a variety of gluco-and sulfo-conjugated metabolites in zebrafish models (primary hepatocytes, ZFL, and ZELH-zfER cell lines) [40]. BaP-diones accounted for the majority of metabolites in the liver tissue of two species of Ictalurid catfish [41]. In fact, increase or decrease in the transcriptome after exposure to exogenous chemicals is not surprising as this phenomenon has been continuously observed in diverse aquatic animals, including fish. However, it is notable that genes mapped to a KEGG pathway by comparative transcriptome analysis showed chemical-specific patterns. The results clearly suggest that transcriptome profiling and transcriptional sensitivity in the P. olivaceus liver tissue is promising as a strong biomarker with distinct characteristics on exposure to certain chemicals.
Vitellogenin (Vtg), as the precursor of egg yolk protein, is important for oocyte development and embryogenesis in oviparous oogenesis, as Vtg is secreted into the bloodstream and incorporated into the growing oocytes [42]. Estrogen and other hormones mainly control the synthesis of Vtg protein, and its extra-hepatic expressions have been detected in teleost in addition to its major production in the liver tissue [43]. Induction of Vtg is known to be a strong biomarker for estrogenic compounds, such as endocrine disrupting chemicals (EDCs) in fish [44][45][46]. Induction of Vtg has also been observed in immature fish (e.g., larvae and juveniles) [47], as shown in the present study. Although BPS has been considered as an alternative of BPA, its estrogenic potential and adverse effects (e.g., endocrine-disrupting effects) on the physiological functions have been consistently reported in vertebrates, including fish [17,48,49]. Developmental exposure to BPS significantly increased plasma Vtg levels in both male and female zebrafish [50]. In fish, vitellogenin production and vitellogenesis are mainly mediated through the induction of ERs [51,52]. In fact, BPS can directly bind to the ERs in a concentration-, tissue-, or species-specific manner [53][54][55][56]. Based on the transcriptional increases of the ER family and plasma E2 level observed in this species, we would conclude that BPS injection triggered vitellogenic responses by exerting estrogenic effects in P. olivaceus.
Zona pellucida proteins (ZPs) share a conserved ZP domain and are constituents of the egg chorion under estrogenic control in fish [57]. Although the potential function of ZPs in the recognition and interaction between egg and sperm has been rarely studied in fish, it is notable that xenoestrogen treatments strongly induced genes associated with fish zonagenesis [42,57]. In addition, mRNA expressions of ZPs (i.e., zp3, zpb, zpc) were significantly increased by the estrogen mimic octylphenol and nonylphenol exposure in the Atlantic salmon smolts and cichlid fish [58,59]. Thus, transcriptional expression of P. olivaceus ZPs can be assertive molecular indicators for xenoestrogens.
Besides the well-characterized carcinogenicity and genotoxicity of BaP, previous studies have suggested that BaP is biotransformed to epoxides and quinones by drug/xenobiotic metabolism via activation of the "AhR gene battery" [18,33,60,61]. In P. olivaceus, BaP injection strongly increased mRNA expression of the members of "AhR gene battery" in a mechanistic manner through phase I (e.g. cyp1a1) to phase II (e.g. ugt2A1, gsts) drug/xenobiotic metabolism with inductions of the EROD and GST activity. The AhR-mediated process was originally characterized as a regulator of PAHs such as BaP beyond the well-characterized TCDD-induced toxicity [62]. In addition to the classical AhR-mediated mechanism, induction of the cyp1a gene at both transcriptional and translational levels has been highlighted as one of the major biomarkers for BaP exposure in numerous teleost [63][64][65][66][67]. Several fish also showed BaP-triggered inductions of members of the UGT [68,69] or GST gene families [41,67,70,71]. Some GST proteins are strongly involved in the AhR-mediated process [33,61]. Taken together, the overall response of the "AhR gene battery" is likely to be sensitive to BaP injection in P. olivaceus. Although our results are not sufficient to prove actual biotransformation because enzymatic activity and chemical structure of BaP metabolites were not analyzed, a significant elevation of the mRNA expressions with induction of the EROD activity indicates their potential as strong biomarkers and early signals of BaP exposure in P. olivaceus.
The mRNA expression and enzymatic activity of genes associated with the antioxidant defense system was increased in the liver transcriptome of P. olivaceus after exposure to both chemicals studied. Over the last ten years, the prooxidant role of BPA has been extensively studied, including BPA-induced intracellular free radicals [e.g. reactive oxygen species (ROS)], oxidative stress, and subsequent DNA damage in numerous in vitro models and vertebrates, including fish [72][73][74]. Although there still needs to be rigorous study to understand BPSinduced oxidative stress, recent research suggests that BPS, as the main BPA substituent, can also promote the generation of ROS and alter antioxidant balance related to oxidative stress in vertebrates [75][76][77]. Taken together, the results of the present study suggest that BPS-triggered induction of oxidative stress may produce potential pleiotropic effects such as DNA damage, developmental toxicity, reproductive toxicity, and disease beyond its endocrine disrupting properties in fish.
Typically, PAHs, including BaP, are known to act as strong mutagens and/or carcinogens via oxidative stress induction. BaP can generate large amounts of intracellular ROS during the biotransformation process and the modulation of the antioxidant defense system [78]. In detail, its metabolites (e.g. BaP-3-phenol, BaP-quinone) are believed to directly bind to intracellular micro and macromolecules with the formation of ROS [79]. In teleost, studies on transcriptional activation and enzymatic modulation of the antioxidant defense system have consistently highlighted their biomarker potential for BaP exposure prior to the phenotypic and physiological toxic effects [67,71,[80][81][82][83]. Thus, we consider that BaP-triggered oxidative stress may be significant in the liver tissues of P. olivaceus. Although the example of KEGG classification, "Protein processing in endoplasmic reticulum", was used for testing different transcriptional sensitivities, the pathway also suggested biological evidence of BaP-triggered oxidative stress, and the subsequent induction of the endoplasmic reticulum-stress pathway, as shown in in vitro systems [84,85]. Since there is very limited information available on ubiquitin ligase (E3) complexes and ubiquitination in teleosts, we were not able to estimate its role in the BaP-specific transcriptional expression observed in this species. Activated AhR induced the expressions of estrogen receptor (ER) target genes, which promoted ER proteolysis through assembling the E3 complex [86]. The BaP-triggered AhR pathway could be associated with the ubiquitination of certain proteins damaged by BaP in the liver tissues of P. olivaceus.
Another interesting aspect revealed by transcriptome profiling is that mRNA expression of genes involved in the innate immune system were reduced by both BPS and BaP injection in this species. Similarly to the decreased mRNA expressions of several complement genes, the haemolytic complement activity was also significantly decreased by both chemicals. The innate immunity of vertebrates has been considered as one of the targets of BPA [87]. BPA modulates various aspects of the immune system, including innate and adaptive immunity, by altering B cell functions, T cell subsets, dendritic cell, and macrophage biology in aquatic animals [87,88]. In teleost, both acute and chronic exposure to BPA significantly affect the immune system such as inhibition of transcriptional expressions of genes related to the immune response [e.g., IFNγ, IL1β, IL10, Mx, TNFα, CC-chemokine, CXCL-clc, the Toll-like receptors (TLRs) signaling pathway (TLR3, TRIF, MyD88, SARM, IRAK4, and TRAF6), and NF-κB-associated immune genes], serum immune parameters [e.g. levels of immunoglobulin M (IgM), C-reaction Protein (CRP), and complement component 3] [87][88][89][90][91]. Although BPS assays on the fish immune system have not been extensively undertaken to date, a recent study suggested that immune defense is diminished by exposing parental generations to environmentally relevant concentrations of BPS in zebrafish [92]. Suppressive effects of BaP on immune systems have been extensively reviewed in teleost [93][94][95]. In particular, the liver tissues contribute to immune surveillance by controlling circulating antigens and innate immune cells [96,97]. Taken together, modulation of immune transcriptomes would be involved in dysregulation of innate immunity, as the physical and chemical barriers of the innate immune system can control numerous waterborne pathogens in fish.
The aim of this study was to elucidate the potential advantages, particularly in monitoring of coastal regions and in transcriptional profiling of the P. olivaceus liver tissue. Finally, we propose a schematic diagram for transcriptional responses involved in the spectrum of flatfish (i.e., chemical-specific) or common induction and inhibition of possible mechanisms in the BPS-or BaP-exposed liver tissues of P. olivaceus with biochemical evidence (Fig 7). Our results show that the application of transcriptome profiling is a promising testing method for the Chemical-specific transcriptome profiles of olive flounder study of exposure to different environmental pollutants, as each transcriptome can respond to certain chemicals via both specific and common response mechanisms.
Supporting information S1 Fig. Comparison of transcriptional expression patterns of whole libraries. A) Hierarchical clustering analysis of each library (i.e. two control liver tissues, three BPS-exposed liver tissues, and three BaP-exposed liver tissues) by employing heat map and hierarchical clustering. B) Overall transcriptional profile of each library. Similarity is depicted with different colors. (TIF) S1