Metabolism disruption analysis of zebrafish larvae in response to BPA and BPA analogs based on RNA-Seq technique

Ecotoxicology and Environmental Safety 174 (2019) 181–188

Contents lists available at ScienceDirect

Ecotoxicology and Environmental Safety journal homepage: www.elsevier.com/locate/ecoenv

Metabolism disruption analysis of zebrafish larvae in response to BPA and BPA analogs based on RNA-Seq technique

T

⁎

Wenhui Qiua,b,1, , Shuai Liuc,1, Feng Yanga,b, Peiyao Dongd, Ming Yangc, Minghung Wonga, ⁎ Chunmiao Zhenga,b, a

Guangdong Provincial Key Laboratory of Soil and Groundwater Pollution Control, School of Environmental Science and Engineering, Southern University of Science and Technology, Shenzhen 518055, China b State Environmental Protection Key Laboratory of Integrated Surface Water-Groundwater Pollution Control, School of Environmental Science and Engineering, Southern University of Science and Technology, Shenzhen 518055, China c School of Environmental and Chemical Engineering, Shanghai University, Shanghai 200444, China d Institute of Water Sciences, College of Engineering, Peking University, Peking 100871, China

A R T I C LE I N FO

A B S T R A C T

Keywords: Bisphenol A Bisphenol F Bisphenol S Zebrafish embryos RNA sequencing Metabolism

Bisphenol A (BPA) is an environmentally ubiquitous chemical widely used in industry and is known to have adverse effects on organisms. Given the negative effect, BPA-free products have been developed with BPA analogs such as bisphenol F (BPF) and bisphenol S (BPS); however, these analogs are proving to exhibit toxicity similar to that of BPA. In the present study, we aimed to identify and compare the underlying mechanisms of toxicity of BPA, BPF, and BPS at the transcriptional level by conducting global transcriptome sequencing (RNASeq) on zebrafish embryos. RNA-seq results showed that the expression levels of 285, 191, and 246 genes were significantly changed in zebrafish larvae after embryos were treated for 120 h with 100 μg/L BPA, BPF, and BPS, respectively. Among the genes exhibiting altered expression, a substantial number were common to two or three exposure groups, suggesting consistent toxicity between the three bisphenols. We further validated the expression levels of 19 differentially expressed genes by qRT-PCR, using sequencing RNA and the RNA samples after treatment by 0.01, 1, and 100 μg/L bisphenols under identical condition, the results were similar to RNASeq. Moreover, functional enrichment analysis indicated that metabolism was the main pathway which disrupted in zebrafish larvae by bisphenols treatment. Protein–protein interaction network analysis indicated that six DEGs (ces, cda, dpyd, upp1, upp2, and cmpk2) interact together in the drug metabolism of zebrafish. In summary, our study revealed changes in the transcription of genes upon bisphenols treatment in zebrafish larvae for the first time, indicating that BPF and BPS may cause adverse effects similar to BPA via their involvement in various biological processes, providing a solid foundation for further research on the toxicology of BPA analogs.

1. Introduction Endocrine-disrupting chemicals (EDCs) are a class of exogenous chemicals that can interfere with hormone function by altering homeostasis or other signaling components of endocrine system (Braun, 2017; Zoeller et al., 2012). As a typical EDC, bisphenol A [2,2-bis(4hydroxyphenyl)propane, BPA], is ubiquitous throughout the ecological environment, at the highest levels of μg/L or μg/g (Corrales et al., 2015), because of its extensive use, random discharge and high leachability (Honeycutt et al., 2017; Huang et al., 2012), for instance, the concentration of BPA spanned from 0.317 to 1275 ng/L in water

treatment plants (Muhamad et al., 2016), and more than 1000 µg/L in landfill leachate and sewage treatment effluents (Crain et al., 2007; Kang et al., 2007). Previous studies have described the adverse effects of BPA on the cardiovascular, metabolic, nervous, and immune systems of organisms (Anbarasu, 2017; Belcher et al., 2015; Burgos-Aceves et al., 2016; Valentino et al., 2016). Therefore, legislation restricting the use of BPA has prompted the development of new alternatives to replace it in industry, such as bisphenol F (4,4′-dihydroxydiphenylmethane, BPF) and bisphenol S (4,4′-sulfonyldiphenol, BPS), which also widespread in the environment as emerging organic contaminants (EOCs) nowadays.

⁎

Corresponding authors at: Guangdong Provincial Key Laboratory of Soil and Groundwater Pollution Control, School of Environmental Science and Engineering, Southern University of Science and Technology, Shenzhen 518055, China. E-mail addresses: [email protected] (W. Qiu), [email protected] (C. Zheng). 1 These authors contributed equally. https://doi.org/10.1016/j.ecoenv.2019.01.126 Received 31 October 2018; Received in revised form 29 December 2018; Accepted 17 January 2019 0147-6513/ © 2019 Published by Elsevier Inc.

Ecotoxicology and Environmental Safety 174 (2019) 181–188

W. Qiu, et al.

2.2. Zebrafish embryos

BPF are structurally similar to BPA (Table S1), it is commonly used to produce lacquers, varnishes, liners, adhesives plastics, and water pipes (Cabaton et al., 2009). Likewise, BPS, another BPA structural analog (Table S1), is also used in baby bottles, thermal receipt paper, food packaging, and personal care products (Qiu et al., 2018c). Similar to BPA, BPF and BPS can also leach from products (Bittner et al., 2014; Zheng et al., 2016) and disperse throughout the environment, such that the occurrences of BPF and BPS have been reported in environment and body fluid. One study (Huang et al., 2018) have summarized the concentrations of bisphenols in the surface water and sediments of some lakes and rivers in China, Japan, and Korea, revealing that the highest concentration of BPA, BPF, and BPS in water were up to 7.50, 2.85 and 65.60 μg/L, and 13.37, 9.65, and 1.97 μg/g in sediments, respectively. Moreover, in indoor dust samples from twelve countries, the mean concentration of BPA, BPF, and BPS were detected at 1, 1, and 0.22 μg/ g (Wang et al., 2015), and also highest at 37.7, 212, and 21 ng/mL in human urine (Liao et al., 2012; Zhou et al., 2014). These results suggest that BPF and BPS are already widely distributed throughout the environment. Furthermore, in humans, using the urine of people living near electronic waste recycling facilities, Zhang et al. (2016) have revealed that the dismantling process of electric waste can cause BPs pollution and exposure, and also proved the BPA and BPS concentration were positively related to the oxidative stress of humans. To sum up, the increasing usage of BPF and BPS and their occurrence in the environment may lead them to acquire pollution status similar to that of BPA, generating speculation regarding their toxic to humans. A substantial number of studies have reported that BPF and BPS possess toxicity similar to that of BPA, which can cause disruption to processes mediated by estrogen receptors (ERs), G protein-coupled estrogen receptors (GPERs), thyroid hormone receptors (THRs), and peroxisome proliferator-activated receptor gamma (PPARγ) in various organisms (Huang et al., 2016; Macczak et al., 2017; Skledar et al., 2016; Vinas and Watson, 2013). Chen et al. (2016) noted that BPS and BPF showed moderate to slightly acute toxicity and weak estrogenic activity in exposed organisms. The endocrine system of zebrafish was also affected by disruption of the hypothalamic-pituitary-gonadal (HPG) or hypothalamic-pituitary-thyroid (HPT) axis following BPF (Yang et al., 2017) and BPS (Zhang et al., 2017) exposure, respectively. Our previous in vitro and in vivo studies also showed that potential toxicity on fish which induced by BPF and BPS was similar with BPA treatment (Qiu et al., 2018a, 2018b, 2018c); moreover, these three bisphenols can alter the immunity of offspring zebrafish to the same degree (Dong et al., 2018). Although a number of studies reported the adverse effects of BPA, BPF, and BPS on zebrafish, changes in the transcriptome profile of exposed zebrafish remained unclear. To understand the underlying toxicity mechanisms of BPA and BPA analogs in zebrafish, global transcriptome sequencing (RNA-Seq) was performed using the embryos of zebrafish (Danio rerio) after exposure to BPA, BPF, and BPS in the present study. The mRNA expression profile and disturbances in related biological processes and pathways were analyzed in the exposed zebrafish larvae. The present study aimed to study the similarity in toxicity of BPA, BPF, and BPS on zebrafish; the results will prove valuable for the further study of BPA analogs.

The adults zebrafish were fed and fertilized zebrafish embryos were obtained according to our previous study (Qiu et al., 2018b) and provided in supplementary information. Those which had reached the blastula stage were considered to have developed normally and were selected for exposure experiments at four hours post-fertilization (hpf). The embryos were kept in a 28 °C incubator in 1 × E3 medium (5 mM NaCl, 0.17 mM KCl, 0.33 mM CaCl2, 0.3 mM MgSO4). In the experimental stage, four hpf zebrafish embryos were randomly distributed in 90 mm glass diameter petri dishes (100 embryos per dish) immediately prior to treatment with chemical exposure solutions (50 mL). 2.3. Experimental design US Environmental Protection Agency has established an oral reference dose of BPA was 100 μg/L in drinking water (Willhite et al., 2008), thus, in the present study, in order to comprehensively study on the mode of action of BPA analogs in zebrafish, the exposure concentration of BPA, BPF, and BPS was set at 100 μg/L. The embryos were exposed to 100 μg/L BPA, BPF, and BPS in E3 medium, using solutions prepared by stepwise dilution of the stock solutions (10 g/L). A solvent control group (0.005% DMSO) incubated in the embryo medium was prepared in parallel with the exposure groups, all groups were set up with three biological parallels. Both control and bisphenol exposure solutions were replaced completely every 24 h; the concentrations of BPA, BPF, and BPS in the exposure solutions after 24 h treatment were detected by LC-MS/MS (Agilent 1260 Infinity, Santa Clara, CA, USA) according to our previous study (Qiu et al., 2016). No significant changes were detected between nominal concentrations and the concentrations detected at 24 h (Table S2). After 120 h exposure at 124 hpf (Qiu et al., 2016, 2018b), the zebrafish larvae in each dish were collected, then immersed in liquid nitrogen to be snap-frozen and stored at − 80 °C for subsequent RNA isolation. The results of hatching rates, survival rates, and body length at 124 hpf were similar to our previous study (Qiu et al., 2018b). Concentration-dependent experiments were performed under identical conditions with the bisphenol exposure concentrations at 0.01, 1, and 100 μg/L. 2.4. RNA isolation, cDNA library construction, and sequencing Total RNA of larvae in DMSO-, BPA-, BPF-, and BPS-treatment groups were extracted by TRNzol Total RNA Extraction Reagents (Tiangen Biotech Co., Ltd., Beijing, China) following the manufacturer's instructions, the RNAs of three replicate in each treatment group were mixed as one sequencing sample. After quality assessment of RNA by NanoDrop™ (Thermo Fisher Scientific, Waltham, MA, USA), Qubit (Invitrogen, Carlsbad, CA, USA), and Agilent 2100 (Agilent, Santa Clara, CA, USA), high quality RNA samples (28 S:18 S=2.0–2.2, RIN > 9.0) were used to construct libraries with a commercial kit (Mega Genomics, Beijing, China). The libraries were then sequenced with 150 cycle paired-end sequencing (300 cycles total) using the Illumina HiSeq PE150 platform (Illumina, San Diego, CA, USA) at Mega Genomics (Beijing, China).

2. Materials and methods

2.5. Bioinformatics analysis of sequencing results

2.1. Chemicals

RNA-Seq experiments were performed using Illumina HiSeq PE150 platform. In total, 28.64 Gb raw data were obtained from DMSO, BPA, BPF, and BPS groups. The RNA-Seq reads were pre-processed using a custom Perl script, based on a filtration algorithm (Garcia et al., 2012) to remove unreliable reads and trim low-quality sections. The processed reads were mapped to reference zebrafish genome using TopHat2 (Kim et al., 2013). All reads were then assembled using Cufflinks (Trapnell et al., 2010) and compared to the reference genome with Cuffcompare to predict new genes; gene annotations were obtained from nr (Deng

BPA (CAS number 80-05-7, 99 +%), BPF (CAS Number 620-92-8, 98 +%), and BPS (CAS Number 80-09-1, 98 +%) were purchased from Sigma-Aldrich (St. Louis, MO, USA). These compounds were dissolved in dimethyl sulfoxide (DMSO, Sigma-Aldrich) to obtain stock solutions of 10 g/L each, all of which were stored at 4 °C until use. All other chemicals used were of analytical grade and purchased from SigmaAldrich. 182

Ecotoxicology and Environmental Safety 174 (2019) 181–188

W. Qiu, et al.

et al., 2006), Swiss-Prot (Apweiler et al., 2004), GO (Ashburner et al., 2000), COG (Tatusov et al., 2000), KOG (Koonin et al., 2004), Pfam (Finn et al., 2014), eggNOG (Huerta-Cepas et al., 2016) and KEGG (Kanehisa et al., 2004) databases by blast (Altschul et al., 1997). The FPKM (Fragments Per Kilobase of transcript per Million mapped reads) values of sequenced genes were counted with the method described by Trapnell et al. (2010). Differential expression analysis was performed using the EBSeq package (Leng et al., 2013). Genes with |log2(fold change)| ≥ 1 and false discovery rate (FDR) < 0.01 were considered significant for differentially expressed genes (DEGs). Finally, functional enrichment analysis of DEGs was carried out using WebGestalt (Wang et al., 2013). The significantly enriched KEGG (Kyoto Encyclopedia of Genes and Genomes) and PANTHER pathways were analyzed using KOBAS 3.0 (Xie et al., 2011). Moreover, in order to compare the 162 gene sets in zebrafish from the KEGG database, the gene sets from the DMSO samples vs ALL bisphenol-treated samples (average fold-change values in BPA, BPF, and BPS) were evaluated using the Java GSEA implementation. A gene set was considered enriched when p-value < 0.05 and FDR < 25%. Protein-protein interaction (PPI) network analysis of DEGs in the three bisphenol groups was then performed using the STRING database (Szklarczyk et al., 2017). 2.6. Validation of transcriptome analysis using quantitative real time PCR (qRT-PCR) To validate the transcriptome analysis of DEGs, the remaining RNA samples from the four treatment groups (DMSO, BPA, BPF, and BPS) and that of in the concentration-dependent experiment were used for qRT-PCR, cDNA templates were synthesized from 1 μg of each sample using a QuantScript RT Kit (Tiangen Biotech Co., Ltd., Beijing, China). The method of qRT-PCR was provided in supplementary information. Rpl13a (ribosomal protein L13a) was used as an internal reference gene for each tested sample according to a previous study (Fan et al., 2010). Specific primers for the selected target genes were synthesized by Invitrogen (USA) (Table S3). Ct-based relative quantification with efficiency correction normalized to rpl13a was calculated using the 2-△△Ct method.

Fig. 1. A, The numbers and lengths of new genes annotated to COG, KEGG, GO, Swissport, KOG, Pfam, eggNOG, and nr databases (L=length); B, The predicted fish species to which the new annotated genes mapped in nr database.

(data not shown). 3.2. Identification of new mRNA genes in zebrafish RNA samples All the sequenced reads were assembled using Cufflinks and compared to the reference genome with Cuffcompare; 26,473 genes were obtained in four zebrafish RNA samples from the annotation to several databases (COG, KEGG, GO, Swissport, KOG, Pfam, eggNOG, and nr), which include 1467 novel genes. Furthermore, all 1467 new genes were compared to the functional annotation databases listed above, with 218, 383, 450, 659, 828, 890, 1312, and 1457 in COG, KEGG, GO, Swissport, KOG, Pfam, eggNOG, and nr, respectively (Fig. 1A) More importantly, the lengths of the majority of annotated genes were more than 300 bp (Fig. 1A). In the nr database annotation results, the quantity of animal species that the new genes annotated was counted and presented in Fig. 1B, with 81.87% annotated in Danio rerio, and only 17.45% in other fish species. This proved the accuracy of the sequencing data in the present study. The GO items and KEGG pathway analysis showed that these new genes have different potential functions in the regulation of various biological processes in zebrafish (Table S5 and Table S6). These results could enrich the annotation file of the zebrafish genome and provide fundamental data for the study of these previously unknown genes.

2.7. Statistical analysis The qRT-PCR data are presented as mean ± standard error of the mean (SEM). The correlations between the RNA-Seq and qRT-PCR results were analyzed by the linear-fitting method on OriginPro 8.0 (OriginLab, Northampton, MA, USA) and by Pearson correlation on SPSS Statistics 20 (SPSS Inc., Chicago, IL, USA). The figures were drawn by OriginPro 8.0, GraphPad Prism 5 (La Jolla, CA, USA), MeV 4.6.0, and Cytsocape 3.6. Statistical analyses were performed by one-way analysis of variance (ANOVA), followed by the least significant difference (LSD) in homogeneity test or Tamhane's in inhomogeneity test. A value of p < 0.05 was used to indicate statistically significant differences with an asterisk (*, n = 3). 3. Results 3.1. RNA sequencing results Zebrafish larvae RNA of DMSO and three bisphenol-exposure groups were subjected to RNA sequencing, each treatment group have one sequencing RNA which mixed from its three replicate samples. After quality control, 27.82 Gb clean data including 46.5, 47.4, 45.5, and 46.2 million clean reads were obtained for DMSO, BPA, BPF, and BPS groups, respectively (Table S4), with 96% of bases generating Qscores ≥ 20 and GC contents ranging from 46.7% to 47.26%. The mapped ratio of clean reads in four groups ranged from 77.42% to 77.85% when referencing with the zebrafish genome database. The saturation test showed that effective data were acquired for all four sample groups

3.3. Identification of DEGs between DMSO group and BPA, BPF, and BPS groups The FPKM values of all detected genes in the DMSO, BPA, BPF, and BPS groups are shown in Fig. S1. In these four treatment groups, the FPKM values of the majority of genes were around 10. In order to evaluate the adverse effects of bisphenol exposure on zebrafish larvae, the DEGs were analyzed, with the BPA, BPF, and BPS exposure groups 183

Ecotoxicology and Environmental Safety 174 (2019) 181–188

W. Qiu, et al.

Fig. 2. The number of DEGs in three bisphenol treatment groups compared to the DMSO control group (A); Hierarchical cluster analysis of DEG expression profiles in three bisphenol treatment groups, B (BPA), C (BPF), and D (BPS). The colors in the heat map represent the expression changes (log2 fold change) in the bisphenol treatment groups relative to the DMSO group, red represents upregulation and green represents downregulation; E, Venn diagram of DEGs between three bisphenol groups compared to the DMSO group. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

assessed against the DMSO group. Following the screening standard, 285, 191, and 246 DEGs were obtained in comparisons between BPA, BPF, and BPS vs DMSO groups, respectively; among these DEGs, there were 53 upregulated and 232 downregulated genes, 34 upregulated and 157 downregulated genes, as well as 43 upregulated and 203 downregulated genes, respectively (Fig. 2A). The expression profile of genes in each DEG set (BPA, BPF, and BPS) in all the three bisphenol treatment groups are evaluated and the hierarchical cluster analysis are shown in Fig. 2B (BPA), C (BPF), and D (BPS). Heat map results reveal that the expression changes in DEGs for each bisphenol group were not all in the same direction in other two bisphenol groups, when compared to the DMSO group. This suggested different toxic effects occurred in zebrafish larvae after treatment with three different bisphenol analogs, although overlapping genes were present in the three DEG sets, which from the results of Venn diagram (Fig. 2E). As shown in Fig. 2E, 20 of the same DEGs were obtained in all three groups of zebrafish treated with different bisphenol analogs, the majority of these 20 genes were both down- or up-regulated in all three treatment groups, except for three: jun dimerization protein 2 (jdp2) up in BPA and BPF groups, while down in BPS group; LOC103910878 and interleukin-1 beta (il1β) up in BPA group, while down in BPF and BPS groups (Table S7). Besides, BPA and BPF dysregulated 41 of same genes, corresponding to 22 genes in comparison with the BPA and BPS treatment groups and 41 genes in comparison with the BPF and BPS treatment groups (Table S7), in which, still have little genes did not in the same change directions. Based on these results, it can be speculated that the three different bisphenols caused both similar and different toxic effects in each of the three treatment groups. Fig. 3. A, qRT-PCR results of 19 DEGs, significant differences versus DMSO control are indicated as *p < 0.05 (ANOVA, LSD's, or Tamhane's test, n = 3); B, Correlation analysis between qRT-PCR and RNA-Seq results, indicating a significantly positive correlation [linear-fitting method (R2) and * *p < 0.01 by Pearson correlation]; C, Hierarchical cluster analysis of 11 DEG expression profiles after treatment by different concentrations of three bisphenols compared to DMSO group (n = 3).

3.4. Validation of DEG expression profiles in zebrafish larvae after bisphenol exposure by qRT-PCR To confirm the results of RNA-Seq, 19 DEGs in the RNA-Seq results were selected to verify their expression patterns by qRT-PCR, using cDNA from the remaining RNA samples in DMSO, BPA, BPF, and BPS treatment groups. The fold changes of these 19 genes (two upregulated and 17 downregulated) reported in the RNA-Seq data are shown in Fig. S2. All of these genes were significantly changed in at least one bisphenol exposure group, and in the same direction of change in the other groups, when compared with the DMSO group. The qRT-PCR results (Fig. 3A) revealed the expression changes of 19 selected genes were similar to the RNA-Seq results, although the magnitude of fold change was more significant in RNA-Seq. In the qRT-PCR results, the expression levels of two genes [interferon-induced very large GTPase 1

(gvinp1) and protein-glutamine gamma-glutamyltransferase K/transglutaminase 1, K polypeptide (tgm1)] were significantly induced by all three bisphenols, except for only a slight increase of gvinp1 in the BPA group. In the other 17 genes, expression levels of 12 genes [three-finger protein 5 (LOC100003647); protein kinase C eta, b (prkch); tribbles pseudokinase 1 (trib1); regulator of G-protein signaling 2 (rgs2); FK506 binding protein 5 (fkbp5); argonaute RISC catalytic component 3b 184

Ecotoxicology and Environmental Safety 174 (2019) 181–188

W. Qiu, et al.

pathway in diabetic complications,’ related to diabetes. To sum up, these results suggest that three different bisphenols (BPA, BPF, and BPS) can severely damage the metabolic system of zebrafish larvae. Six genes in BPA DEGs set were further explored, all of which were enriched in the ‘Metabolic pathways,’ ‘Drug metabolism-other enzymes,’ and ‘Pyrimidine metabolism’ pathways, including five downregulated genes [cytidine deaminase (cda); dihydropyrimidine dehydrogenase (dpyd); uridine phosphorylase 1 (upp1); uridine phosphorylase 2 (upp2); and thymidine phosphorylase (tymp)] and one upregulated gene [cytidine monophosphate (UMP-CMP) kinase 2, mitochondrial (cmpk2)] (Fig. 4A). Of even greater interest, PPI network analysis showed that five enzymes (CDA, DPYD, UPP1, UPP2, and CMPK2) encoded by the genes mentioned above (cda, dpyd, upp1, upp2, and cmpk2) interacted with another protein, CES2 [encoded by a downregulated gene, carboxylesterase 2 (ces2)] involved only in the ‘Drug metabolism-other enzymes’ pathway (Fig. 4C), to participate in related physiological functions in zebrafish larvae. In Fig. 4D, qRT-PCR results showed relatively similar expression changes in the levels of cmpk2, cda, upp1, and upp2 expression when compared with RNA-Seq results, although no concentration-response relationships were found after 0.01, 1, and 100 μg/L bisphenol treatments.

(ago3b); epidermal growth factor receptor kinase substrate 8-like protein 1 (eps8l1); immunoglobulin superfamily (igsf); potassium channel subfamily K member 1 (kcnk1); complement component 7b (c7b); tetraspanin35 (tspan35); and solute carrier family 16, member 9a (slc16a9a)] in zebrafish larvae were significantly reduced after BPA, BPF, and BPS treatment. However, the expression levels of the five remaining genes [Epstein-Barr virus induced gene 3 (ebi3); chemokine (C-X-C motif) ligand 18b (cxcl18b); retinol dehydrogenase 10b (rdh10b); serine palmitoyltransferase, long chain base subunit 2b (sptlc2b); and protein kinase C delta b (prkcdb)] were not all significantly decreased. For instance, ebi3 and prkcdb expression levels were increased in BPS-exposed zebrafish larvae. The correlations of expression fold changes in the 19 DEGs between RNA-Seq and qRT-PCR results were analyzed by the linear-fitting (R2) and Pearson correlation methods, with R2 (0.691, 0.647, and 0.499) and Pearson correlation coefficients (0.842, 0.817, and 0.726) in BPA, BPF, and BPS groups, which illustrated in Fig. 3B. Results indicated significantly positive correlations for all three different bisphenols, demonstrating that the transcriptional changes determined by qRT-PCR were consistent with results obtained by RNA-Seq, thereby establishing the validity of the RNA-Seq data. In order to study the concentration-response relationship for 11 selected genes, zebrafish embryos were then exposed to three different concentrations of bisphenol (0.01, 1, and 100 μg/L) under conditions identical to those set in the RNA-seq experiment (Fig. 3C and Table S8). In which, two genes (gvinp1 and tgm1) were upregulated in all bisphenol-exposed groups with strong concentration-response relationship. In the other nine genes, expression levels for all of them were reduced after BPA, BPF, and BPS treatments, except for slight increases in some relatively low concentration exposure groups. Thus, the RNAseq results were again validated, using different exposure samples, which suggests that BPA, BPF, and BPS indeed have adverse effects on zebrafish embryos/larvae by altering the mRNA expression levels of genes related to diverse biological processes.

4. Discussion In the present study, RNA-Seq was used to investigate the similar effects among BPA, BPF, and BPS on zebrafish larvae for the first time, and our results revealed that exposure to these three bisphenols can alter the expression of several genes in zebrafish larvae. Obviously, larvae is more sensitive when exposed to BPA than to BPF or BPS in same concentration, however, the negative effect is relatively similar among these three bisphenols treatment, which involved in three important pathways: ‘Apoptosis signaling pathway’, ‘AGE-RAGE signaling pathway in diabetic complications’, and ‘Metabolic pathways’, and this may proved that BPA-substitutes were not pretty safe to organisms. PANTHER pathway analysis indicated that ‘Apoptosis signaling pathway’ was significantly affected by bisphenol exposure, with dysregulation of c-Fos and hsp70. As a proto-oncogene, c-Fos may be induced by numerous stimuli to interact with c-JUN to regulate cell proliferation and differentiation (Filardo et al., 2002; Shaulian and Karin, 2002). BPA has been proven to induce the expression of c-Fos by GPER, leading to cell apoptosis (Pupo et al., 2012; Wang et al., 2017). According to the RNA-seq results in the present study, c-Fos was also significantly increased in BPA- and BPF-treated groups (column chart in Fig. 4A), which in line with previous studies. HSPs (heat shock proteins) exhibit strong cytoprotective effects, maintaining proper protein folding, disaggregation, and refolding of misfolded proteins, as well as targeting damaged proteins for degradation (Kong et al., 2011; Tutar and Tutar, 2010), hsp70 also inhibits cell apoptosis by interacting with or regulating the expression of bcl-2 (Kim et al., 2002). BPA, BPF, and BPS significantly reduced the expression of hsp70 in zebrafish larvae in the present study (column chart in Fig. 4A), suggested that cell apoptosis and DNA damage occurred in zebrafish larvae after bisphenol exposure, which well complemented the similar report that BPA and BPS could induce the apoptosis of zebrafish larvae or red common carp macrophages (Qiu et al., 2018c; Wu et al., 2017). Another enriched KEGG pathway ‘AGE-RAGE signaling pathway in diabetic complications’ is related to advanced glycation end products (AGE) and receptor of AGE (RAGE), which plays an important role in the pathogenesis of diabetes (Moser et al., 2007; Ramasamy et al., 2011), all three bisphenols changed at least one gene's mRNA level (Fig. 4A), indicating a potential disruption, suggesting that bisphenol could also be involved in diabetes via this pathway. Most interesting, various metabolism system of zebrafish larvae were severely affected by bisphenol treatment (Fig. 4), in the past researches, it has been reported that BPA induced metabolic disruption in organisms through complex interactions with ERs, THRs, and

3.5. Functional enrichment analysis of DEGs in bisphenol treatment groups To further investigate the underlying effects of bisphenol treatments on the DEGs identified in zebrafish larvae, GO categories of DEGs were analyzed with WebGestalt. Among all DEGs in three bisphenol treatment groups, 9, 6, and 4 different GO items were significantly enriched (Table S9); all GO items were different in three bisphenol treatment groups. Moreover, the expression profiles of annotated genes in bisphenol-exposed groups were counted. As shown in Table S9, the expression changes of enriched genes in the BPA exposed group were almost the same as those in the BPF, while distinct from the BPS. These findings suggest that the toxic effects of BPA and BPF on zebrafish were more similar than BPS, though the common dysregulation of several genes was detected in all three bisphenol treatment groups. The KEGG pathway enrichment results were analyzed using KOBAS 3.0 (Fig. 4A). According to the results, 48 different DEGs were significantly enriched in 10 KEGG pathways and one PANTHER pathway (p < 0.05). The only enriched PANTHER pathway, ‘Apoptosis signaling pathway,’ is related to two genes, fosab (c-Fos) in BPA and BPF and heat shock protein 70 (hsp70) in all three bisphenol groups. In the KEGG enrichment results, ‘Metabolism’ was the main enriched KEGG pathway in all DEGs of the three bisphenol treatment groups, which includes nine pathways: ‘Metabolic pathways,’ ‘Drug metabolism-other enzymes,’ ‘Pyrimidine metabolism,’ ‘Arachidonic acid metabolism,’ ‘Glycine, serine and threonine metabolism,’ ‘Glyoxylate and dicarboxylate metabolism,’ ‘Tryptophan metabolism,’ ‘Steroid biosynthesis,’ and ‘Fatty acid biosynthesis.’ Meanwhile, from the GSEA results (Fig. 4B and Fig. S3), the DEG sets related to ‘Metabolic pathways’ were significantly enriched in the KEGG pathway, similar to the results obtained from KOBAS, which proved the consistency and accuracy of our analysis. The remaining pathway in the KEGG analysis was the ‘AGE-RAGE signaling 185

Ecotoxicology and Environmental Safety 174 (2019) 181–188

W. Qiu, et al.

Fig. 4. A, KEGG pathway analysis of all DEGs in BPA, BPF, and BPS treatment groups, the different colors in the circle represent different DEG sets; B, enrichment pathway results of DEGs by GSEA; C, PPI network of selected proteins; D, gene expression levels of cmpk2, cdab, upp1, and upp2, after treatment by BPA, BPF, and BPS at 0.01, 1, and 100 μg/L; significant differences versus DMSO control are indicated as *p < 0.05 (ANOVA, LSD's, or Tamhane's test). The four column charts in A and C are the mRNA expression changes of four genes in RNA-seq results, respectively, * represent the significant changes of gene mRNA level.

metabolism, reflecting the adverse effects on metabolic activity that occurred in zebrafish larvae. Carboxylesterase is one type of detoxifying enzyme that enhances the degradation of xenobiotics in organisms which has proven to be an appropriate biomarker in polluted environments (Jeon et al., 2016; Wu et al., 2014). In human fetal liver,

glucocorticoid receptors (GRs), which result in adipogenesis, weight gain, diabetes, or insulin level changes (Alonso-Magdalena et al., 2010; Casals-Casas and Desvergne, 2011; Schug et al., 2011). In the present study, the PPI network analysis showed the interactions of six enzymes (CES2, CDA, UPP1, UPP2, DPYD, and CMPK2) related to drug

186

Ecotoxicology and Environmental Safety 174 (2019) 181–188

W. Qiu, et al.

Conflict of interest

high BPA decreased the expression of ces2 and ces5a (Nahar et al., 2014). Similarly, BPA, BPF, and BPS also reduced ces2 mRNA expression levels in zebrafish larvae (column chart in Fig. 4A). Reduced levels of CES2 suggest that bisphenol exposure may affect the normal detoxification functions in zebrafish larvae. The following four enzymes (CDA, UPP1, UPP2, and DPYD) all have the ability to metabolize purines or pyrimidines (Rengaraj et al., 2013). As a liver enzyme, CDA catalyzes the inactivation of cytidine and dCyd to uridine and deoxyuridine, respectively (Galmarini et al., 2001), and resistance to external nucleosidic analog drugs. In cancer cells, an increase of CDA could accelerate the metabolism of therapeutic drugs to weaken their curative effects (Galmarini et al., 2001; Serdjebi et al., 2015), while a decrease of CDA may reduce the metabolic activity of external drugs or other organic substances in normal organisms or cells, thereby inducing toxicity. After CDA-related processes, uridine and deoxyuridine were catalyzed to uracil by UPP1/UPP2 and TYMP; meanwhile, DPYD is an important enzyme in the process by which uracil changes to β-alanine and thymine changes to 3-aminoisobutanoate (Galmarini et al., 2001). All the processes described above are part of the catabolic pathway of pyrimidine metabolism; for instance, the synthesis and metabolism of the anti-cancer drug 5-fluorouracil require the presence of CDA, UPP1, UPP2, TYMP, and DPYD (Lam et al., 2016). CMPK is one type of monophosphate kinase which participates in the de novo and salvage pathways of pyrimidine metabolism. Xu et al. reported a UMP-CMPK2 enzyme that could phosphorylate dUMP, dCMP, UMP and CMP, with a preference for the deoxynucleotides (Xu et al., 2008). Therefore, alteration of these six genes suggested the disruption of metabolism system of zebrafish larvae. Moreover, our GSEA analysis (Fig. 4B and Fig. S3) further confirmed the results that bisphenol exposures indeed reduce the metabolic ability of zebrafish larvae and affect their normal development when respond to drugs, which is consistent with previous report (Fenichel et al., 2013; Hélièstoussaint et al., 2014; Sharpe and Drake, 2010).

The authors declare no competing financial interest. Appendix A. Supporting information Supplementary data associated with this article can be found in the online version at doi:10.1016/j.ecoenv.2019.01.126. References Alonso-Magdalena, P., Ropero, A.B., Soriano, S., Quesada, I., Nadal, A., 2010. BisphenolA: a new diabetogenic factor? Horm.-Int. J. Endocrinol. 9 (2), 118–126. Altschul, S.F., Madden, T.L., Schaffer, A.A., Zhang, J.H., Zhang, Z., Miller, W., Lipman, D.J., 1997. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 25 (17), 3389–3402. Anbarasu, C.R., 2017. Effects of Bisphenol A Exposure on Central Nervous System Development. Apweiler, R., Bairoch, A., Wu, C.H., Barker, W.C., Boeckmann, B., Ferro, S., Gasteiger, E., Huang, H.Z., Lopez, R., Magrane, M., Martin, M.J., Natale, D.A., O'Donovan, C., Redaschi, N., Yeh, L.S.L., 2004. UniProt: the universal protein knowledgebase. Nucleic Acids Res. 32, D115–D119. Ashburner, M., Ball, C.A., Blake, J.A., Botstein, D., Butler, H., Cherry, J.M., Davis, A.P., Dolinski, K., Dwight, S.S., Eppig, J.T., 2000. Gene ontology: tool for the unification of biology. Nat. Genet. 25 (1), 25. Belcher, S.M., Gear, R.B., Kendig, E.L., 2015. Bisphenol A alters autonomic tone and extracellular matrix structure and induces sex-specific effects on cardiovascular function in male and female CD-1 mice. Endocrinology 156 (3), 882–895. Bittner, G.D., Yang, C.Z., Stoner, M.A., 2014. Estrogenic chemicals often leach from BPAfree plastic products that are replacements for BPA-containing polycarbonate products. Environ. Health 13 (1), 41. Braun, J.M., 2017. Early-life exposure to EDCs: role in childhood obesity and neurodevelopment. Nat. Rev. Endocrinol. 13 (3), 161–173. Burgos-Aceves, M.A., Cohen, A., Smith, Y., Faggio, C., 2016. Estrogen regulation of gene expression in the teleost fish immune system. Fish Shellfish Immunol. 58, 42–49. Cabaton, N., Dumont, C., Severin, I., Perdu, E., Zalko, D., Cherkaoui-Malki, M., Chagnon, M.-C., 2009. Genotoxic and endocrine activities of bis(hydroxyphenyl) methane (bisphenol F) and its derivatives in the HepG2 cell line. Toxicology 255 (1–2), 15–24. Casals-Casas, C., Desvergne, B., 2011. Endocrine disruptors: from endocrine to metabolic disruption. Annu. Rev. Physiol. 73, 135–162. Chen, D., Kannan, K., Tan, H.L., Zheng, Z.G., Feng, Y.L., Wu, Y., Widelka, M., 2016. Bisphenol analogues other than BPA: environmental occurrence, human exposure, and toxicity - a review. Environ. Sci. Technol. 50 (11), 5438–5453. Corrales, J., Kristofco, L.A., Steele, W.B., Yates, B.S., Breed, C.S., Williams, E.S., Brooks, B.W., 2015. Global assessment of bisphenol A in the environment: review and analysis of its occurrence and bioaccumulation. Dose Response 13 (3) (1559325815598308). Crain, D.A., Eriksen, M., Iguchi, T., Jobling, S., Laufer, H., Leblanc, G.A., Guillette, L.J., 2007. An ecological assessment of bisphenol-A: evidence from comparative biology. Reprod. Toxicol. 24 (2), 225–239. Deng, Y., Li, J., Wu, S., Zhu, Y., Chen, Y., He, F., 2006. Integrated nr database in protein annotation system and its localization. Comput. Eng. 32 (5), 71–74. Dong, X., Zhang, Z., Meng, S., Pan, C., Yang, M., Wu, X., Yang, L., Xu, H., 2018. Parental exposure to bisphenol A and its analogs influences zebrafish offspring immunity. Sci. Total Environ. 610–611, 291–297. Fan, C.Y., Cowden, J., Simmons, S.O., Padilla, S., Ramabhadran, R., 2010. Gene expression changes in developing zebrafish as potential markers for rapid developmental neurotoxicity screening. Neurotoxicol. Teratol. 32 (1), 91–98. Fenichel, P., Chevalier, N., Bruckerdavis, F., 2013. Bisphenol A: an endocrine and metabolic disruptor. Ann. Endocrinol. 74 (3), 211–220. Filardo, E.J., Quinn, J.A., Frackelton Jr., A.R., Bland, K.I., 2002. Estrogen action via the G protein-coupled receptor, GPR30: stimulation of adenylyl cyclase and cAMP-mediated attenuation of the epidermal growth factor receptor-to-MAPK signaling axis. Mol. Endocrinol. 16 (1), 70–84. Finn, R.D., Bateman, A., Clements, J., Coggill, P., Eberhardt, R.Y., Eddy, S.R., Heger, A., Hetherington, K., Holm, L., Mistry, J., Sonnhammer, E.L.L., Tate, J., Punta, M., 2014. Pfam: the protein families database. Nucleic Acids Res. 42 (D1), D222–D230. Galmarini, C.M., Mackey, J.R., Dumontet, C., 2001. Nucleoside analogues: mechanisms of drug resistance and reversal strategies. Leukemia 15 (6), 875–890. Garcia, T.I., Shen, Y.J., Crawford, D., Oleksiak, M.F., Whitehead, A., Walter, R.B., 2012. RNA-Seq reveals complex genetic response to deepwater horizon oil release in Fundulus grandis. BMC Genom. 13 (1), 474. Hélièstoussaint, C., Peyre, L., Costanzo, C., Chagnon, M.C., Rahmani, R., 2014. Is bisphenol S a safe substitute for bisphenol A in terms of metabolic function? An in vitro study. Toxicol. Appl. Pharmacol. 280 (2), 224–235. Honeycutt, J.A., Nguyen, J.Q.T., Kentner, A.C., Brenhouse, H.C., 2017. Effects of water bottle materials and filtration on bisphenol A content in laboratory animal drinking water. J. Am. Assoc. Lab. Anim. Sci. 56 (3), 269–272. Huang, C., Wu, L.H., Liu, G.Q., Shi, L., Guo, Y., 2018. Occurrence and ecological risk assessment of eight endocrine-disrupting chemicals in urban river water and sediments of south china. Arch. Environ. Contam. Toxicol. 1–12. Huang, G.M., Tian, X.F., Fang, X.D., Ji, F.J., 2016. Waterborne exposure to bisphenol F

5. Conclusion The RNA-Seq results reported in this study reveal that mRNA level of 285, 191, and 246 genes in zebrafish embryos/larvae exposed to bisphenols (BPA, BPF, and BPS) was significantly altered. Based on the bioinformatic analysis results, multiple physiological processes and signaling pathways in zebrafish were affected by bisphenols treatment, especially the metabolism systems. Furthermore, the expression levels of metabolism system-related DEGs (ces2, cda, dpyd, upp1, upp2, and cmpk2) in zebrafish larvae were also changed in groups treatment by different concentrations of bisphenol. Our findings reveal that BPAsubstitutes BPF and BPS may exert comparative similar toxicities on non-target aquatic organisms when compared with BPA. Therefore, concerns about the safety of BPA analogs warrant further study. Acknowledgment This work was sponsored by National Key R&D Program of China (2018YFC0406504), National Natural Science Foundation of China (Grants 21707064), State Key Laboratory of Marine Environmental Science, the MEL Visiting Fellowship Program, State Environmental Protection Key Laboratory of Integrated Surface Water-Groundwater Pollution Control, Guangdong Provincial Key Laboratory of Soil and Groundwater Pollution Control (No. 2017B030301012), Shenzhen Municipal Science and Technology Innovation Committee through project Shenzhen Key Laboratory of Soil and Groundwater Pollution Control (No. ZDSY20150831141712549), Shenzhen Science and Technology Innovation Committee (KQTD2016022619584022, 201803023000902), Southern University of Science and Technology (Grant No. G01296001), and the Leading Talents of Guangdong Province Program (Chunmiao Zheng). 187

Ecotoxicology and Environmental Safety 174 (2019) 181–188

W. Qiu, et al.

chemicals and disease susceptibility. J. Steroid Biochem. Mol. Biol. 127 (3–5), 204–215. Serdjebi, C., Milano, G., Ciccolini, J., 2015. Role of cytidine deaminase in toxicity and efficacy of nucleosidic analogs. Expert Opin. Drug Metab. Toxicol. 11 (5), 665–672. Sharpe, R.M., Drake, A.J., 2010. Bisphenol A and metabolic syndrome. Endocrinology 151 (6), 2404–2407. Shaulian, E., Karin, M., 2002. AP-1 as a regulator of cell life and death. Nat. Cell Biol. 4 (5), E131–E136. Skledar, D.G., Schmidt, J., Fic, A., Klopcic, I., Trontelj, J., Dolenc, M.S., Finel, M., Masic, L.P., 2016. Influence of metabolism on endocrine activities of bisphenol S. Chemosphere 157, 152–159. Szklarczyk, D., Morris, J.H., Cook, H., Kuhn, M., Wyder, S., Simonovic, M., Santos, A., Doncheva, N.T., Roth, A., Bork, P., Jensen, L.J., von Mering, C., 2017. The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res. 45, D362–D368. Tatusov, R.L., Galperin, M.Y., Natale, D.A., Koonin, E.V., 2000. The COG database: a tool for genome-scale analysis of protein functions and evolution. Nucleic Acids Res. 28 (1), 33–36. Trapnell, C., Williams, B.A., Pertea, G., Mortazavi, A., Kwan, G., van Baren, M.J., Salzberg, S.L., Wold, B.J., Pachter, L., 2010. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28 (5), 511–515. Tutar, L., Tutar, Y., 2010. Heat shock proteins; an overview. Curr. Pharm. Biotechnol. 11 (2), 216–222. Valentino, R., D’Esposito, V., Ariemma, F., Cimmino, I., Beguinot, F., Formisano, P., 2016. Bisphenol A environmental exposure and the detrimental effects on human metabolic health: is it necessary to revise the risk assessment in vulnerable population? J. Endocrinol. Invest. 39 (3), 259–263. Vinas, R., Watson, C.S., 2013. Mixtures of xenoestrogens disrupt estradiol-induced nongenomic signaling and downstream functions in pituitary cells. Environ. Health 12 (1), 26. Wang, C.L., Zhang, J.X., Li, Q., Zhang, T.B., Deng, Z.S., Lian, J., Jia, D.H., Li, R., Zheng, T., Ding, X.J., Yang, F., Ma, C., Wang, R., Zhang, W.X., Wen, J.G., 2017. Low concentration of BPA induces mice spermatocytes apoptosis via GPR30. Oncotarget 8 (30), 49005–49015. Wang, J., Duncan, D., Shi, Z., Zhang, B., 2013. WEB-based gene set analysis toolkit (WebGestalt): update 2013. Nucleic Acids Res. 41 (W1), W77–W83. Wang, W., Abualnaja, K.O., Asimakopoulos, A.G., Covaci, A., Gevao, B., JohnsonRestrepo, B., Kumosani, T.A., Malarvannan, G., Minh, T.B., Moon, H.B., Nakata, H., Sinha, R.K., Kannan, K., 2015. A comparative assessment of human exposure to tetrabromobisphenol A and eight bisphenols including bisphenol A via indoor dust ingestion in twelve countries. Environ. Int. 83, 183–191. Willhite, C.C., Ball, G.L., McLellan, C.J., 2008. Derivation of a bisphenol A oral reference dose (RfD) and drinking-water equivalent concentration. J. Toxicol. Environ. Health B. Crit. Rev. 11 (2), 69–146. Wu, H.H., Gao, C., Guo, Y.P., Zhang, Y.P., Zhang, J.Z., Ma, E.B., 2014. Acute toxicity and sublethal effects of fipronil on detoxification enzymes in juvenile zebrafish (Danio rerio). Pestic. Biochem. Physiol. 115, 9–14. Wu, M.H., Pan, C.Y., Chen, Z., Jiang, L.H., Lei, P.H., Yang, M., 2017. Bioconcentration pattern and induced apoptosis of bisphenol A in zebrafish embryos at environmentally relevant concentrations. Environ. Sci. Pollut. Res. 24 (7), 6611–6621. Xie, C., Mao, X.Z., Huang, J.J., Ding, Y., Wu, J.M., Dong, S., Kong, L., Gao, G., Li, C.Y., Wei, L.P., 2011. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 39, W316–W322. Xu, Y.J., Johansson, M., Karlsson, A., 2008. Human UMP-CMP kinase 2, a novel nucleoside monophosphate kinase localized in mitochondria. J. Biol. Chem. 283 (3), 1563–1571. Yang, Q., Yang, X.H., Liu, J.N., Ren, W.J., Chen, Y.W., Shen, S.B., 2017. Effects of BPF on steroid hormone homeostasis and gene expression in the hypothalamic-pituitarygonadal axis of zebrafish. Environ. Sci. Pollut. Res. 24 (26), 21311–21322. Zhang, D.H., Zhou, E.X., Yang, Z.L., 2017. Waterborne exposure to BPS causes thyroid endocrine disruption in zebrafish larvae. Plos One 12 (5), e0176927. Zhang, T., Xue, J., Gao, C., Qiu, R.L., Li, Y.X., Huang, M.Z., Kannan, K., 2016. Urinary concentrations of bisphenols and their association with biomarkers of oxidative stress in people living near e-waste recycling facilities in China. Environ. Sci. Technol. 50, 4045–4053. Zheng, S., Shi, J.C., Hu, J.Y., Hu, W.X., Zhang, J., Shao, B., 2016. Chlorination of bisphenol F and the estrogenic and peroxisome proliferator-activated receptor gamma effects of its disinfection byproducts. Water Res. 107, 1–10. Zhou, X., Kramer, J.P., Calafat, A.M., Ye, X., 2014. Automated on-line column-switching high performance liquid chromatography isotope dilution tandem mass spectrometry method for the quantification of bisphenol A, bisphenol F, bisphenol S, and 11 other phenols in urine. J. Chromatogr. B: Anal. Technol. Biomed. Life Sci. 944, 152–156. Zoeller, R.T., Brown, T.R., Doan, L.L., Gore, A.C., Skakkebaek, N.E., Soto, A.M., Woodruff, T.J., Saal, F.S.V., 2012. Endocrine-disrupting chemicals and public health protection: a statement of principles from the endocrine society. Endocrinology 153 (9), 4097–4110.

causes thyroid endocrine disruption in zebrafish larvae. Chemosphere 147, 188–194. Huang, Y.Q., Wong, C.K.C., Zheng, J.S., Bouwman, H., Barra, R., Wahlstrom, B., Neretin, L., Wong, M.H., 2012. Bisphenol A (BPA) in China: a review of sources, environmental levels, and potential human health impacts. Environ. Int. 42, 91–99. Huerta-Cepas, J., Szklarczyk, D., Forslund, K., Cook, H., Heller, D., Walter, M.C., Rattei, T., Mende, D.R., Sunagawa, S., Kuhn, M., Jensen, L.J., von Mering, C., Bork, P., 2016. eggNOG 4.5: a hierarchical orthology framework with improved functional annotations for eukaryotic, prokaryotic and viral sequences. Nucleic Acids Res. 44 (D1), D286–D293. Jeon, H.J., Lee, Y.H., Kim, M.J., Choi, S.D., Park, B.J., Lee, S.E., 2016. Integrated biomarkers induced by chlorpyrifos in two different life stages of zebrafish (Danio rerio) for environmental risk assessment. Environ. Toxicol. Pharmacol. 43, 166–174. Kanehisa, M., Goto, S., Kawashima, S., Okuno, Y., Hattori, M., 2004. The KEGG resource for deciphering the genome. Nucleic Acids Res. 32, D277–D280. Kang, J.H., Aasi, D., Katayama, Y., 2007. Bisphenol A in the aquatic environment and its endocrine-disruptive effects on aquatic organisms. Crit. Rev. Toxicol. 37 (7), 607–625. Kim, D., Pertea, G., Trapnell, C., Pimentel, H., Kelley, R., Salzberg, S.L., 2013. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 14 (4), R36. Kim, R., Tanabe, K., Uchida, Y., Emi, M., Inoue, H., Toge, T., 2002. Current status of the molecular mechanisms of anticancer drug-induced apoptosis - The contribution of molecular-level analysis to cancer chemotherapy. Cancer Chemother. Pharmacol. 50 (5), 343–352. Kong, X.C., Zhang, D., Qian, C., Liu, G.T., Bao, X.Q., 2011. FLZ, a novel HSP27 and HSP70 inducer, protects SH-SY5Y cells from apoptosis caused by MPP. Brain Res. 1383, 99–107. Koonin, E.V., Fedorova, N.D., Jackson, J.D., Jacobs, A.R., Krylov, D.M., Makarova, K.S., Mazumder, R., Mekhedov, S.L., Nikolskaya, A.N., Rao, B.S., Rogozin, I.B., Smirnov, S., Sorokin, A.V., Sverdlov, A.V., Vasudevan, S., Wolf, Y.I., Yin, J.J., Natale, D.A., 2004. A comprehensive evolutionary classification of proteins encoded in complete eukaryotic genomes. Genome Biol. 5 (2), R7. Lam, S.W., Guchelaar, H.J., Boven, E., 2016. The role of pharmacogenetics in capecitabine efficacy and toxicity. Cancer Treat. Rev. 50, 9–22. Leng, N., Dawson, J.A., Thomson, J.A., Ruotti, V., Rissman, A.I., Smits, B.M.G., Haag, J.D., Gould, M.N., Stewart, R.M., Kendziorski, C., 2013. EBSeq: an empirical Bayes hierarchical model for inference in RNA-seq experiments. Bioinformatics 29 (8), 1035–1043. Liao, C., Liu, F., Alomirah, H., Loi, V.D., Mohd, M.A., Moon, H.B., Nakata, H., Kannan, K., 2012. Bisphenol S in urine from the United States and seven Asian countries: occurrence and human exposures. Environ. Sci. Technol. 46 (12), 6860–6866. Macczak, A., Cyrkler, M., Bukowska, B., Michalowicz, J., 2017. Bisphenol A, bisphenol S, bisphenol F and bisphenol AF induce different oxidative stress and damage in human red blood cells (in vitro study). Toxicol. In Vitro 41, 143–149. Moser, B., Desai, D.D., Downie, M.P., Chen, Y., Yan, S.F., Herold, K., Schmidt, A.M., Clynes, R., 2007. Receptor for advanced glycation end products expression on T cells contributes to antigen-specific cellular expansion in vivo. J. Immunol. 179 (12), 8051–8058. Muhamad, M.S., Salim, M.R., Lau, W.J., Yusop, Z., 2016. A review on bisphenol A occurrences, health effects and treatment process via membrane technology for drinking water. Environ. Sci. Pollut. Res. 23 (12), 11549–11567. Nahar, M.S., Kim, J.H., Sartor, M.A., Dolinoy, D.C., 2014. Bisphenol A-associated alterations in the expression and epigenetic regulation of genes encoding xenobiotic metabolizing enzymes in human fetal liver. Environ. Mol. Mutagen. 55 (3), 184–195. Pupo, M., Pisano, A., Lappano, R., Santolla, M.F., De Francesco, E.M., Abonante, S., Rosano, C., Maggiolini, M., 2012. Bisphenol A induces gene expression changes and proliferative effects through GPER in breast cancer cells and cancer-associated fibroblasts. Environ. Health Perspect. 120 (8), 1177–1182. Qiu, W., Zhan, H., Tian, Y., Zhang, T., He, X., Luo, S., Xu, H., Zheng, C., 2018a. The in vivo action of chronic bisphenol F showing potential immune disturbance in juvenile common carp (Cyprinus carpio). Chemosphere 205, 506–513. Qiu, W.H., Shao, H.Y., Lei, P.H., Zheng, C.M., Qiu, C.X., Yang, M., Zheng, Y., 2018b. Immunotoxicity of bisphenol S and F are similar to that of bisphenol A during zebrafish early development. Chemosphere 194, 1–8. Qiu, W.H., Yang, M., Liu, S., Lei, P.H., Hu, L., Chen, B., Wu, M.H., Wang, K.J., 2018c. Toxic effects of bisphenol S showing immunomodulation in fish macrophages. Environ. Sci. Technol. 52 (2), 831–838. Qiu, W.H., Zhao, Y.L., Yang, M., Farajzadeh, M., Pan, C.Y., Wayne, N.L., 2016. Actions of bisphenol A and bisphenol S on the reproductive neuroendocrine system during early development in zebrafish. Endocrinology 157 (2), 636–647. Ramasamy, R., Yan, S.F., Schmidt, A.M., 2011. Receptor for AGE (RAGE): signaling mechanisms in the pathogenesis of diabetes and its complications. Ann. N. Y. Acad. Sci. 1243 (1), 88–102. Rengaraj, D., Lee, B.R., Jang, H.J., Kim, Y.M., Han, J.Y., 2013. Comparative metabolic pathway analysis with special reference to nucleotide metabolism-related genes in chicken primordial germ cells. Theriogenology 79 (1), 28–39. Schug, T.T., Janesick, A., Blumberg, B., Heindel, J.J., 2011. Endocrine disrupting

188