The chemical defensome of five model teleost fish

How an organism copes with chemicals is largely determined by the genes and proteins that collectively function to defend against, detoxify and eliminate chemical stressors. This integrative network includes receptors and transcription factors, biotransformation enzymes, transporters, antioxidants, and metal- and heat-responsive genes, and is collectively known as the chemical defensome. Teleost fish is the largest group of vertebrate species and can provide valuable insights into the evolution and functional diversity of defensome genes. We have previously shown that the xenosensing pregnane x receptor (pxr, nr1i2) is lost in many teleost species, including Atlantic cod (Gadus morhua) and three-spined stickleback (Gasterosteus aculeatus), but it is not known if compensatory mechanisms or signaling pathways have evolved in its absence. In this study, we compared the genes comprising the chemical defensome of five fish species that span the teleosteii evolutionary branch often used as model species in toxicological studies and environmental monitoring programs: zebrafish (Danio rerio), medaka (Oryzias latipes), Atlantic killifish (Fundulus heteroclitus), Atlantic cod, and three-spined stickleback. Genome mining revealed evolved differences in the number and composition of defensome genes that can have implication for how these species sense and respond to environmental pollutants, but we did not observe any candidates of compensatory mechanisms or pathways in cod and stickleback in the absence of pxr. The results indicate that knowledge regarding the diversity and function of the defensome will be important for toxicological testing and risk assessment studies.

www.nature.com/scientificreports/ subfunctionalizations. For example, we have previously shown that several losses of the pregnane x receptor (pxr, or nr1i2) have occurred independently across teleost evolution 9 . PXR is an important xenosensor and as a ligandactivated transcription factor one of the key regulators of the chemical defensome 10,11 . The importance of PXR in response to chemical stressors in vertebrates, raises questions of how some fish species cope without this gene. Thus, the objective of this study was to compare the chemical defensome of zebrafish (Danio rerio), medaka (Oryzias latipes), and Atlantic killifish (Fundulus heteroclitus), which are species that have retained a pxr gene, to Atlantic cod (Gadus morhua) and three-spined stickleback (Gasterosteus aculeatus), that have lost this gene by independent mechanisms. Genomes of zebrafish, medaka, killifish, and stickleback were chosen based on their widespread use as model species in both developmental and toxicological studies [12][13][14][15] . Atlantic cod is an ecologically and economically important species in the North Atlantic Ocean, and has commonly been used as a bioindicator species in environmental monitoring programs [16][17][18] . The genome of Atlantic cod was first published in 2011 19 , which has facilitated its increased use as model in toxicological studies [20][21][22][23][24] . Although these fish are all benthopelagic species, their natural habitats range from freshwater to marine environments, from tropical to temperate conditions, and reach sizes ranging from less than five cm (zebrafish and medaka) and 15 cm (killifish and stickleback) to 200 cm (Atlantic cod).
Our results include an overview of the full complement of putative chemical defensome genes in these five fish species. Furthermore, by using previously published transcriptomics data, we investigated the transcriptional expression of defensome genes in early development of zebrafish and stickleback, and their responses to the model polycyclic aromatic hydrocarbon contaminant, benzo(a)pyrene.
In conclusion, our study represents the first interspecies comparison of the full complement of chemical defensome genes in teleost model species. We found that although most defensome genes have been retained in the teleost genomes over millions of years, there are distinct differences between the species. Based on our results, we suggest a holistic approach to analyze omics datasets from toxicogenomic studies, where differences in the chemical defensome gene complement are taken into consideration.
Two main approaches were used to identify the genes related to the chemical defensome of the five fish species. First, we searched the current annotations in NCBI for Atlantic cod or ENSEMBL for zebrafish, stickleback, killifish, and medaka using gene names. For the well-annotated zebrafish genome, this approach successfully identified the genes that are part of the chemical defensome, as previously mapped by Stegeman, et al. 3 . Alternative gene sequences, e.g. splice variants, were removed when identifying the number of genes within the different gene families, but included in the RNA-Seq analyses.
However, only relying on annotations will not identify all defensome genes in the other less characterized fish genomes. Thus, secondly, we also performed hidden Markov model (HMM) searches using HMMER and Pfam profiles representing protein families that are part of the chemical defensome (available at FAIRDOM-Hub: https:// doi. org/ 10. 15490/ faird omhub.1. datafi le. 3956.1) in the genomes of the remaining four fish species. Putative orthologs of the retrieved protein sequences were identified using reciprocal best hit BLAST searches against the well-annotated zebrafish proteome. To capture any species-specific duplications in the fish genomes compared to the zebrafish reference genome, hits from one-way BLAST hits were also included. The identified peptide sequence IDs were subsequently converted to their related gene IDs using the BioMart tool on ENSEMBL (https://m. ensem bl. org/ bioma rt/ martv iew) and R package "mygene" (https:// doi. org/ doi: 10. 18129/ B9. bioc. mygene). Finally, the resulting gene lists were refined to contain only members of gene families and subfamilies related to the chemical defensome, using the same defensome gene lists as in the first approach (https:// doi. org/ 10. 15490/ faird omhub.1. datafi le. 3957.1).
Transcription of chemical defensome genes in early development. No embryos were directly used in this study. RNA-Seq datasets of early developmental stages of zebrafish (expression values in Transcripts Per Million from ArrayExpress: E-ERAD-475) and stickleback (sequencing reads from NCBI BioProject: PRJNA395155) were previously published by White, et al. 25 and Kaitetzidou, et al. 26 , respectively.
For stickleback, embryos were sampled at early morula, late morula, mid-gastrula, early organogenesis, and 24 h post hatching (hph). The sequencing data was processed and analyzed following the automatic pipeline RASflow 27 . Briefly, the reads were mapped to the stickleback genome downloaded from ENSEMBL of version release-100. HISAT2 28 was used as aligner and featureCounts 29 was used to count the reads. The library sizes were normalized using Trimmed Mean of M values (TMM) 30 and the Counts Per Million (CPM) were calculated using R package edgeR 31 . The source code and relevant files can be found on GitHub: https:// github. com/ zhxia okang/ fishD efens ome/ tree/ main/ devel opmen talSt ages/ stick leback/ RASfl ow. www.nature.com/scientificreports/ The zebrafish dataset included 18 time points from one cell to five days post fertilization (dpf). In order to best compare to the available stickleback developmental data, we chose to include the following time points in this study: Cleavage_2 cell (early morula), blastula_1k cell (late morula), mid-gastrula, segmentation_1-4-somites (early organogenesis), and larval_protruding_mouth (24 hph).
The temporal clustering was performed on both species. The R package "tscR" 32 , which clusters time series data using slope distance was applied here, to account for the gene expression variation pattern over time, instead of the absolute gene expression values.
Exposure response of defensome genes in zebrafish early development. RNA-Seq datasets of zebrafish exposed to benzo(a)pyrene (B(a)P) (gene counts from NCBI GEO: GSE64198) were previously published by Fang, et al. 33 . Briefly, the datasets are results from the following in vivo experiments: Adult parental zebrafish were waterborne-exposed to 50 μg/L B(a)P or 0.1 mL/L ethanol vehicle control, and their eggs collected from days 7 to 11. The eggs were further raised in control conditions or continuously exposed to 42.0 ± 1.9 μg/L B(a)P until 3.3 and 96 (hours post fertilization, hpf), where mid-blastula state and complete organogenesis is reached, respectively. At 3.3 or 96 hpf, embryos (50/pool) or larvae (15/pool) were pooled for RNA isolation, giving three replicate groups of control and exposed at 3.3 hpf and two replicate groups of control and exposed at 96 hpf.
The RNA-Seq reads were mapped to the zebrafish genome and the genes were quantified in the original publication by Fang, et al. 33 . In our study, the gene counts were then normalized followed by differential expression analysis using edgeR 31 and the chemical defensome genes were then identified based on our defensome gene list of zebrafish.

Results and discussion
Chemical defensome genes present in model fish genomes.. The full complement of chemical defensome genes in zebrafish, killifish, medaka, stickleback and Atlantic cod are available at the FAIRDOMHub (https:// doi. org/ 10. 15490/ faird omhub.1. datafi le. 3958.1). In short, genome analyses of the selected fish species show that the number of chemical defensome genes range from 446 in stickleback to 510 in zebrafish ( Fig. 1). Although the number of putative homologous genes in each subfamily varies, we found that all gene subfamilies of the chemical defensome is represented in each species, except for the absence of pxr in stickleback and Atlantic cod. While we believe that all major protein families that are related to the chemical defensome are included, we cannot preclude that there are individual genes that are not identified due to the genome quality and level of annotation.

Soluble receptors and transcription factors.
Stress-activated transcription factors serve as important first responders to many chemicals, and in turn regulate the transcription of other parts of the chemical defensome. Nuclear receptors (NRs) are a superfamily of structurally similar, ligand-activated transcription factors, where members of subfamilies NR1A, B, C, H, and I (such as retinoid acid receptors, peroxisomal proliferatoractivated receptors, and liver X receptor), NR2A and B (hepatocyte nuclear factors and retinoid x receptors), and NR3 (such as estrogen receptors and androgen receptor) are involved in the chemical defense [34][35][36][37] . All NR subfamilies were found in the five fish genomes, except for the nr1i2 gene (Supplementary Table 2, available at FAIRDOMHub: https:// doi. org/ 10. 15490/ faird omhub.1. docum ent. 925.3). We did not find an expansion of any of the groups that indicated an evolved compensational receptor in the absence of pxr.
NR1I2, or pregnane x receptor (PXR) is considered an important xenosensor responsible for the transcription of many genes involved in the biotransformation of xenobiotics 10,38 . We have previously shown that loss of the pxr gene has occurred multiple times in teleost fish evolution 9 , including in Atlantic cod and stickleback. Interestingly, our searches did not reveal a pxr gene in the ENSEMBL genome assembly of Japanese medaka HdRr, which is considered the reference medaka strain 39,40 . In contrast, a pxr gene is annotated in the ENSEMBL genomes of the closely related medaka strains HSOK and HNI. Our previous study identified a sequence similar to zebrafish Pxr in the MEDAKA1 (ENSEMBL release 93) genome 9 , and a partial coding sequence (cds) of pxr is cloned from medaka genome 41 . However, the specific strain of these resources is not disclosed.
To assess the possible absence of pxr in the medaka HdRr strain, we also performed synteny analysis. In vertebrate species, including fish, pxr is flanked by the genes maats1 and gsk3b. These genes are also annotated in medaka HdRr, but the specific gene region has a very low %GC and low sequence quality. Thus, we suspect that the absence of pxr in the Japanese medaka HdRr genome is due to a sequencing or assembly error. However, until more evidence of its absence can be presented, we chose to include the medaka pxr gene (UniProt ID A8DD90_ORYLA) in our resulting list of medaka chemical defensome genes.
Other important transcription factors are the basic helix-loop-helix Per-Arnt-Sim (bHLH/PAS) proteins and the oxidative stress-activated transcription factors. bHLH/PAS proteins are involved in circadian rhythms (such as clock and arntl), hypoxia response (such as hif1a and ncoa), development (such as sim), and the aryl hydrocarbon receptor pathway (ahr and arnt) (as reviewed by Gu, et al. 42 and Kewley, et al. 43 ), while oxidative stress-activated transcription factors respond to changes in the cellular redox status and promote transcription of antioxidant enzymes 44,45 . The latter protein family includes nuclear factor erythroid-derived 2 (NFE2), NFE2-like (NFE2L, also known as NRFs) 1, 2, and 3, BACH, the dimerization partners small-Mafs (MafF, MafG and MafK), and the inhibitor Kelch-like-ECH-associated protein 1 (KEAP1) 46 . Putative orthologs for all these subfamilies were identified in all five fish genomes.
Biotransformation enzymes. In the first phase of xenobiotic biotransformation, a set of enzymes modifies substrates to more hydrophilic and reactive products. The most important gene family of oxygenases is the www.nature.com/scientificreports/ cytochrome P450 enzymes (CYPs, EC 1.14. -.-), a large superfamily of heme-proteins that initiate the biotransformation of numerous xenobiotic compounds through their monooxygenase activity 47 . The subfamilies considered to be involved in xenobiotic transformation is Cyp1, Cyp2, Cyp3, and Cyp4, and genes of these families were found in all fish species. Interestingly, our results show that the zebrafish genome holds almost twice as many genes belonging to the cyp2 subfamily compared to the other four fish species. This is in line with previous findings comparing the CYPome in zebrafish to that in human and cod 48,49 , but the reason for the cyp2 gene bloom in zebrafish remains unknown. The number of genes in each subfamily differs slightly from previous mappings of the CYPome of zebrafish and cod 48,49 (Supplementary Table 3 1.1), and prostaglandin-endoperoxide synthases (PTGS, also known as cyclooxygenases, EC 1.14.99.1). Of these, aldh represented the largest family in our study, with number of putative gene orthologs ranging from 19 to 22 genes. In comparison, fmo had only one gene in zebrafish and four in killifish and cod.
Furthermore, reductases modify chemicals by reducing the number of electrons. Reductases include aldo-keto reductases (AKRs, EC 1.1.1), hydroxysteroid dehydrogenases (HSDs, EC 1.1.1), epoxide hydrolases (EPHXs, EC 3.3.2.9 and EC 3.3.2.10), and the NAD(P)H:quinone oxidoreductases (NQOs, EC 1.6.5.2). Interestingly, the number of putative orthologous genes in the nqo reductase families varied greatly between the fish species, ranging from one in zebrafish to ten in stickleback. A phylogenetic analysis of the evolutionary relationship of the sequences (Fig. 2), shows that all fish species have a Nqo1 annotated gene. In addition, medaka, killifish, stickleback and cod have three to nine other closely related genes. The endogenous functions of the different nqo genes found in fish, and thus the consequences of their putative evolutionary gain in teleost fish, remains unknown and should be studied further.
In the second phase of biotransformation, endogenous polar molecules are covalently attached to xenobiotic compounds by transferases, thus facilitating the generation of more water-soluble products that can be excreted from the cells and the organism. An important class of such conjugating enzymes are glutathione S-transferases (GSTs, EC 2.5.1.18), which are divided into three superfamilies: the cytosolic GSTs (divided in six subfamilies designated alpha through zeta), the mitochondrial GST (GST kappa) and the membrane-associated GST (designated MAPEG) 52,53 . Other classes of conjugating enzymes in vertebrates include cytosolic sulfotransferases Our searches identified a gstp gene in zebrafish and cod genomes, but not in medaka, killifish and stickleback. Gstp is previously identified as the major Gst isoenzyme in livers of marine salmonid species 54 . Although the specificity of GSTP is not fully understood, its activity seems related to oxidative stress 55 . Furthermore, the total number of gst encoding genes was substantially higher in zebrafish (19 genes) compared to the other fish species (9 in stickleback, 10 in medaka, and 13 in killifish and in cod).
Similarly, the number of ugt encoding genes is considerably higher in zebrafish compared to the other fish genomes. Whereas 31 ugt genes were found annotated in zebrafish, we only identified 15 in killifish, 11 genes in medaka, 17 in stickleback, and 16 in cod. All Ugt subfamilies (ugt1, -2, -5, and -8) are represented in the different fish genomes but include a varying number of homologs. Previous publications indicate that the number of zebrafish ugt encoding genes are as high as 45, with the ugt5 subfamily only existing in teleost and amphibian species 56 . One study has found cooperation of NQO1 and UGT in detoxification of vitamin K3 in HEK293 cell line 57 . However, it is not known if there is any correlation between the high number of ugt genes and low number of nqo genes in zebrafish, relative to the other fish species.
Transporter proteins. Energy-dependent efflux transport of compounds across both extra-and intracellular membranes is facilitated by ATP-binding cassette (ABC) transporters. In humans, the ABCs are organized into seven subfamilies, named ABC A through G, where proteins of B (also known as MDR1), C and G are known to be involved in multidrug resistance (MDR) 58 . Our results showed that the number of abc genes in each subfamily were quite similar between the fish, except for Atlantic cod that showed a lower number of abcb and abcc compared to the other species (Supplementary Table 4, available at FAIRDOMHub: https:// doi. org/ 10. 15490/ faird omhub.1. docum ent. 925.3). Abcc1, also known as multidrug resistance-associated protein (mrp1) were identified in all species, whereas abcb1, also known as multidrug resistance protein (mdr1), were only observed in killifish. Neither was abcb5, which is suggested as a close ortholog to abcb1 in zebrafish, found in medaka, stickleback and cod, which is in line with previous results 59 . A separate group of proteins is called the Solute Carrier (SLC) 'superfamily, ' which consists of diverse non-homologous groups of ion and metal transporting membrane proteins that facilitate passive transport 60 . Relevant solute carrier proteins include the drug transporting SLC22 and SLC47, the zinc transporting SLC30 and SLC39, the copper transporting SLC31, and the organic anion transporting SLCO 61 .
We found that the number of abcb, abcc, abcg, slc, and slco genes was similar between the fish species, with zebrafish holding a slightly higher number of homologs. MDR and P-glycoproteins have been relatively understudied in fish [62][63][64] . A clade of abch transporters related to abcg is found in some fish species, including zebrafish, but not in Japanese medaka, stickleback and cod 65 . However, as the endogenous function of these genes are not determined, we have not included them specifically into this study.
Together with xenobiotic metabolizing enzymes, induction of genes and enzymatic activity involved in antioxidant defense has long been recognized as a gold standard in the biomarker approach to environmental studies 69 . Putative orthologs for all antioxidant genes were clearly identified in the five fish genomes examined.
Heat-responsive genes. Heat-responsive genes represent the largest functional group of genes in the chemical defensome and act in response to a wide range of endogenous and exogenous stressors, such as temperatureshock and heavy metal exposure 70 . In response to stressors, heat shock factors (HSFs) regulate transcription of heat shock proteins (HSPs) 71 . HSPs are divided into families based on their molecular size, and each subfamily has various cellular tasks, including cytoskeleton modulation, protein folding, and chaperone functioning 72,73 .
Not much is known about the heat shock protein expression in fish 70 . We found that all fish genomes hold putative orthologs of hsf, heat shock binding proteins (hsbp) and hsp. The number of putative orthologs of hsp and dnaj (formerly known as hsp40) is high in all fish genomes, with the highest number in killifish with 40 and 60 genes, respectively.
In our study, we found putative orthologs for all gene families, except a mt encoding gene in stickleback or cod genome assemblies. Metallothioneins are cysteine-rich, low molecular weight proteins, and can thus be lost due to low-quality sequence and subsequent assembly. In discrepancy with the genome data, Mts are previously described in both Atlantic cod (Hylland et al. 1994) and stickleback (Uren Webster 2017), and these protein IDs were included into our overview. Similarly, only one or two mt genes were found in zebrafish, killifish, and medaka, and this low number is in line with previous findings on metallothioneins in fish 75 . www.nature.com/scientificreports/ Expression of defensome genes during early development of fish. The developmental stage at which a chemical exposure event occurs greatly impacts the effect on fish. In general, chemical exposures during early developmental stages of fish cause the most adverse and detrimental effects. Based on data from the ECETOC Aquatic Toxicity database, fish larvae are more sensitive to substances than embryos and juveniles 76 . However, it is not known how the sensitivity is correlated to the expression of the chemical defensome. As examples in this study, we mapped the expression of the full complement of chemical defensome genes during early development using previously published transcriptomics data from zebrafish 25 and stickleback 26 (Fig. 3, relevant data available on FAIRDOMHub: https:// doi. org/ 10. 15490/ faird omhub.1. assay. 1379.1).
Our results showed that many defensome genes are not expressed until after hatching in both species. The delayed genes belonged to all functional categories but were especially prominent where there are several paralogs within the same gene subfamily, for example the transporters and the transferases (Fig. 3). In contrast, the heat shock proteins hspa8 and hsp90ab1 were highly transcribed at all stages in both zebrafish and stickleback (Fig. 3).
Hspa8 is a constitutively expressed member of the Hsp70 subfamily, which is previously known as important in rodent embryogenesis 77 . The role of Hsp90ab1 in development is less known 78 .
Using temporal clustering analysis (relevant data available on FAIRDOMHub: https:// doi. org/ 10. 15490/ faird omhub.1. assay. 1379.2), we further investigated whether the absence of pxr in stickleback lead to a shift in target gene expression during development. Our results did not find clustering patterns of genes related to pxr, or any other transcription factor, in neither fish. Cyp3a65, a well known zebrafish Pxr target gene 79 , were transcribed at the late morula stage in stickleback and at the mid gastrula stage in zebrafish. This difference could possibly be explained by the difference in synteny of cyp3a65 in zebrafish and cyp3a65 in a number of fish species 49 . In a previous study, zebrafish cyp3a65 appear to have bimodal expression 49 , but this was not observed in the set of developmental stages included in our study. For the glutathione transferases, we found substantial variations in expression between different genes, which is also previously shown at the protein expression and enzyme activity level 80,81 . Thus, although it has been demonstrated that genes regulated by common transcription factors tend to be located spatially close in the genome sequence and thus facilitate a concerted gene expression 82 , we did not observe such patterns in our datasets. Exposure response of the defensome genes. Next, we studied the transcriptional effect of waterborne exposure to 42.0 ± 1.9 μg/L benzo(a)pyrene (BaP), a well-known Ahr agonist, on the chemical defensome genes at embryonic and larval stages of zebrafish (relevant data available on FAIRDOMHub: https:// doi. org/ 10. 15490/ faird omhub.1. datafi le. 3961.1). In the original publication of the RNA-Seq data, Fang, et al. 33 showed that the number of differentially expressed genes increased from 8 at the embryonic stage (3.3. hpf) to 1153 at the larvae stage (96 hpf), where the affected pathways included the Ahr detoxification pathway. Focusing on the chemical defensome genes, we show that following exposure at the larval stage, but not embryonic stage, we found a trend of clustered regulation of functionally grouped genes (Fig. 4). In general, transcription factors were downregulated, whereas biotransformation enzymes were upregulated. However, the BaP xenosensor, ahr2, and the oxidative stress-responsive transcription factor nfe2l1a, were both upregulated (1.4-and 2.4-fold, respectively). The crosstalk between these transcription factors following exposure to chemical stressors is previously studied in zebrafish 83,84 . Furthermore, we showed that at the embryonic stage the exposure led to a strong upregulation of cyp2aa9 (14-fold) and an increased transcription of single genes such as ahr2, nfe2, aanat1, and hspb1 (Fig. 4a). These effects were not identified using the differential expression analysis in the original study but can indicate low level changes in the chemical defensome signaling network. Cyp1a, which is an established biomarker of exposure to BaP and other polycyclic aromatic hydrocarbons (PAH) in fish [85][86][87] , was not induced at this stage. However, as described in the original study 33 , we found a strong induction of cyp1a (50-fold) at the larvae developmental stage (Fig. 4b). Induction of zebrafish cyp1a is previously shown from 24 hpf following exposure to the Ahr model-agonist TCDD 88 .

Summary and perspectives
The chemical defensome is essential for detoxification and subsequent clearance of xenobiotic compounds, and the composition of the defensome can determine the toxicological responses to many chemicals. Based on the previous annotations of the chemical defensome in other species, our results showed that the number of chemical defensome genes ranged from 446 in three-spined stickleback to 510 in zebrafish due to a varying number of gene homologs in the evolutionarily conserved modules. Of the five fish included in this study, zebrafish has the highest number of gene homologs in most gene families, with the interesting exception of the nqo reductases where medaka, killifish, cod, and especially stickleback, had retained a higher number of homologs compared to only one in zebrafish.
We have previously shown that the stress-activated receptor pxr gene has been lost in stickleback and cod, but is retained in zebrafish, Atlantic killifish and medaka 9 . The main objective of this study was to explore whether other compensatory genes or signaling pathways had evolved in its absence. However, based on our results, we did not observe such compensatory mechanisms. Still, there could be other variables linked to the absence of pxr that fell outside the scope of this study, such as sequence divergences, regulatory elements and expression variation.
Furthermore, we analyzed the transcriptional levels of the defensome genes in early development of zebrafish and stickleback. Importantly, the full complement of defensome genes was not transcribed until after hatching. Comparing the transcriptional effects of waterborne exposure to BaP, a well-known Ahr agonist, on chemical defensome genes at two developmental stages of zebrafish, we found a trend for clustering of gene functional groups at the larvae, but not embryonic, stage.
This study presents characterization of the chemical defensome in five different fish species and at different developmental stages as a way of illustrating and understanding genome-based interspecies and stage-dependent  www.nature.com/scientificreports/ differences in sensitivity and response to chemical stressors. As this study has focused on the presence or absence of chemical defensome genes, and the variation across model teleosts, the role of intraspecies, strain-dependent variants in defensome genes were not included. Several studies have identified defensome gene variants linked to pollution tolerance in fish populations, e.g. in the Ahr pathway [89][90][91][92] . In the study by Lille-Langøy, et al. 93 , we previously showed that single-nucleotide polymorphisms (SNPs) in the zebrafish pxr gene affect ligand activation patterns. Thus, sequence variation and variation in expression and inducibility also play important roles in chemical sensitivity differences between species, and SNP and copy number variation contributes to sensitivity differences within species. Traditionally, studying single molecular biomarkers of exposure has proven very useful in toxicological studies 69,94 . Now, the recent advances in omics technologies enable a more holistic view of toxicological responses, including gene set enrichment analysis and pathway analysis approaches [95][96][97][98] . However, these analyses can be challenging when working with less studied and annotated species, such as marine teleosts. As seen from our results, studying the full gene complement of the chemical defense system can identify trends of grouped responses that can provide a better understanding of the overall orchestrated effects to chemical stressors, e.g. applied in genome-scale metabolic reconstructions 99 . Such insights will be highly useful in chemical toxicity testing and environmental risk assessment.  ) and (b) zebrafish larvae (96 hpf) following exposure to benzo(a)pyrene. The transcription is shown as log2 fold change between exposed and control group at each timepoint. The genes are grouped into their functional categories in the chemical defensome and the name of some genes are indicated for clarity.