Small CYP complement and no detectable CYP1-like gene

BLAST Searches using CYP1 family protein sequences from various vertebrates and the invertebrates did not result in CYP- like sequences in O. dioica geneome. We then searched the genome for all CYP sequences and unambiguously identified 23 unique genes (Additional file 1, Figure S1, Additional file 2). Thus, Oikopleura seems to have the smallest CYP gene complement of animal genomes sequenced thus far [27], http://drnelson.uthsc.edu/CytochromeP450.html. Phylogenetic analysis shows that none of the Oikopleura CYPs cluster with CYP1 family proteins from other organisms (Figure 1). A larger set of proteins and species including all the 23 Oikopleura CYPs and many CYP1 protein sequences from various organisms [28] were also used to construct a phylogenetic tree of similar topology to the one shown in Figure 1, but none of the Oikopleura CYPs cluster with CYP1 sequences (Additional file 1, Figure S1). Thus, CYP1 genes are either lost or have diverged beyond recognition in Oikopleura. In contrast, CYP1-like genes are present in the genomes of Ciona and sea urchin [2, 28]. Oikopleura has CYP3-like and CYP4-like genes but the 13 CYP2-like sequences do not cluster well with CYP2 family genes from other organisms (Additional file 1, Figure S1). The absence of CYP1-like genes appears consistent with the parallel absence of AhR signaling pathway (see below). Interestingly, our expression analysis showed that BaP and Clo induced CYP3-like-2 (CBY21750) and CYP2-like-11 (CBY24358) genes, respectively (Additional files 3 and 4).

Lack of detection of AhR signaling pathway genes

In mammals, bHLH-PAS domain proteins sense xenobiotics (AhR), oxygen (hypoxia-inducible factors/HIFs) and light (circadian locomoter output cycles protein kaput/clock) [18]. BLAST search in the Oikopleura genome and EST collections did not reveal a candidate ortholog for AhR. HMM searches for PAS domain proteins in predicted Oikopleura protein database and BLAST searches using various bHLH-PAS domain protein sequences of vertebrates and invertebrates identified three bHLH-PAS domain proteins in Oikopleura. One PAS domain protein identified (accession: CBY20245) is a potassium voltage-gated channel family like protein (not a member of bHLH-PAS domain proteins) and was not analyzed further. Phylogenetic comparison of the three Oikopleura bHLH-PAS domain proteins with representative homologs from mouse, Ciona and sea urchin showed none of them cluster with AhR (Figure 2). The failure to detect AhR and its repressor AhRR genes, as well as the downstream CYP1 family genes (Table 1) suggests the absence of an AhR mediated xenobiotic biotransformation signaling pathway in Oikopleura. In contrast, AhR orthologs are present in other invertebrates including the tunicate C. intestinalis[1, 2, 20, 28].

Oikopleura NR1I-like and NR1H-like nuclear receptors

The absence of detected AhR sequence prompted us to inspect the Oikopleura genome for candidate xenobiotic sensing nuclear receptors. There are about 40 candidate nuclear receptors in Oikopleura (our unpublished data), suggesting significant lineage-specific expansion of the superfamily. SXR, CAR and VDR (vitamin D receptor) belong to NR1I subfamily of nuclear receptor genes. SXR and CAR serve as xenobiotic sensors that activate the transcription of some CYP genes in vertebrates [11, 12]. In vertebrates, the NR1H subfamily of nuclear receptors liver-X-receptor (LXR) and farnesoid-X-receptor (FXR) are activated by sterols and bile acids, respectively [29]. LXR and FXR, also found in invertebrates [2], are potential targets of xenobiotic compounds. Genes encoding NR1I-like and NR1H-like nuclear receptors appear to be multiplied in the Oikopleura genome, even though the precise orthology relationships could not be clarified [10]. There are 6 NR1I and 10 NR1H (LXR/FXR-like) genes in the genome (Table 1) [10], and among these, there could be candidate xenobiotic receptors. No strong candidate for PPAR (NR1C) was found in Oikopleura (our unpublished data).

Nrf2 mediated oxidative stress response pathway genes in Oikopleura

We further explored the Oikopleura genome for transcription factors involved in oxidative stress response, particularly the Nrf2-complex and found homologs of Nrf2, kelch-like ECH-associated protein 1 (Keap1), and small Maf proteins (Table 1, Additional file 1, Figure S2). Phylogenetic comparison of the O. dioica proteins with representative proteins from other organisms suggested a candidate Nrf2-homolog gene for Oikopleura (accession: CBY09280) although it clusters with both Nrf2 and Bach proteins (Additional file 1, Figure S2). Nrf2 is a member of cap 'n' collar-basic leucine-zipper (cnc-bZIP) proteins that heterodimerize with small Maf family proteins and activate transcription from anti-oxidant response elements (ARE) in promoters of anti-oxidant defense genes [30]. Nrf2 is normally bound to its suppressor Keap-1 and is inactive in the cytosol, but it dissociates under oxidative stress from Keap-1 and activates the transcription of target genes [23]. Mammalian Bach 1 and Bach 2 proteins are Nrf2-related factors that also heterodimerize with Maf proteins and mediate oxidative stress signaling [31]. Thus, although uncertainties remain due to the divergence of Oikopleura sequences, components of Nrf2 signaling appear to be present.

Other groups of transcription factors involved in the oxidative stress response include the heat shock factors (HSFs) [32]. In Oikopleura, there are two or three HSF5-like genes (not shown) and one gene very similar to both HSF1 and HSF4 (here referred to as HSF1/4) (CBY25079). The HSF1/4 and many HSP genes were strongly induced by BaP treatment (Additional file 1, Table S1, Additional file 3), and HSF binding sites could be identified in the promoter regions of many of the BaP induced HSP genes (not shown).

Genes differentially regulated by BaP

BaP significantly changed the expression level of 762 annotated genes (336 genes up-regulated and 426 genes down-regulated) in O. dioica (Additional file 3). The list of affected genes includes several genes that are involved in drug metabolism and disposition. Xenobiotic Phase I and II biotransformation enzyme genes induced by BaP include genes encoding a cytochrome P450 family 3 (CYP3-like-2, CBY21750), carboxylesterase, glutathione S-transferase alpha (GSTα), thiosulfate sulfurtransferase and thioredoxin. Genes encoding Phase III transport proteins, ABC transporters and solute carrier proteins were also modulated. In BaP-treated animals, a co-regulation of genes in oxidative stress pathways was demonstrated by up-regulation of GSTα and the two subunits of glutamate-cysteine ligase (GCL), a rate limiting enzyme in glutathione synthesis (Additional file 1, Figure S3, Additional file 3). Another example of co- regulation is the simultaneous up-regulation of many heat-shock protein (HSP) genes and their transcription factor HSF1/4 gene [32, 33] (Additional file 1, Table S1). HSP genes were among the most strongly induced genes (Additional file 1, Table S1, Additional file 3). Some putative HSP proteins of Oikopleura showed only low sequence similarity to mouse proteins (BLASTP E-value ≥ 0.001) and therefore were not included in the list of genes used for pathway analysis (see Methods).

Pathways affected by BaP

In order to take advantage of the rich annotations of mammalian model organisms such as mouse, all predicted proteins of genes differentially regulated were annotated using BLAST searches against M. musculus proteome (BLASTP, e-value < 0.001). The best mouse hit for each Oikopeleura gene was then used in functional analyses using DAVID [34, 35] and various tools in MetaCore (GeneGo) [36]. It is therefore important to interpret the pathway analyses data with caution since some pathways present in mammalian models may not be conserved or relevant in Oikopleura. For example, enriched disease pathways are not relevant in Oikopleura, although they may give insights into mechanisms of the underlying expression changes of the homolog genes in Oikopleura.

MetaCore enrichment analyses showed most enriched BaP affected pathways are involved in apoptosis, oxidative stress response, immune response, protein folding and cell adhesion (Table 3A and 3B). The top scored pathway map (Apoptosis and survival_Endoplasmic reticulum stress response) is shown in Figure 3 with the BaP affected genes indicated. In MetaCore Interactome analysis option, networks can be built to see interconnectedness within a data set. Genes in the input list can have physical and functional interactions (such as binding, catalytic, phosphorylation, transcriptional regulation) that are used by the algorithm to build sub-networks ranked by relative enrichment with the uploaded genes and relative saturation of networks with canonical pathways [36]. A detailed legend for MetaCore pathway maps and network objects is shown in Additional file 5. The network building tool in MetaCore was used to create networks using the Transcription Regulation algorithm, and resulted in a list of ranked transcription factor centered networks (Additional file 6). One of the BaP-induced key transcription factors, JunD is a functional component of the AP-1 transcription factor complex. AP-1 is one of the top ranked transcription factors (Additional file 6) and many oxidative stress response genes including GCLm, GCLc and thioredoxin that are up-regulated by BaP are in the AP-1 network (Figure 4).

Rank	A. Pathway Maps	pValue
1	Apoptosis and survival_ER stress response pathway	6.5E-06
2	NRF2 regulation of oxidative stress response	7.7E-06
3	Role of alpha-6/beta-4 integrins in carcinoma progression	1.5E-05
4	Mechanisms of CFTR activation by S-nitrosoglutathione	1.8E-05
5	Immune response_IL-12 signaling pathway	1.8E-05
6	Oxidative stress_Role of ASK1 under oxidative stress	1.9E-05
7	Development_TGF-beta-dependent induction of EMT via MAPK	2.1E-05
8	Development_Glucocorticoid receptor signaling	2.4E-05
9	Cytoskeleton remodeling_TGF, WNT and cytoskeletal remodeling	2.4E-05
10	Immune response_IL-7 signaling in T lymphocytes	4.1E-05
	B. Networks
1	Protein folding_Response to unfolded proteins	2.9E-08
2	Cell adhesion_Platelet-endothelium-leucocyte interactions	9.0E-06
3	Protein folding_ER and cytoplasm	3.8E-05
4	Development_Blood vessel morphogenesis	4.7E-05
5	Reproduction_Spermatogenesis, motility and copulation	5.1E-05
6	Cell adhesion_Cell-matrix interactions	5.2E-05
7	Signal transduction_ESR1-nuclear pathway	6.3E-05
8	Development_Regulation of angiogenesis	2.8E-04
9	Proteolysis_Connective tissue degradation	5.3E-04
10	Signal transduction_NOTCH signaling	6.0E-04

To explore interactions within the set of BaP modulated genes, networks were generated in MetaCore using Analyze algorithm with default settings. In the top ranking networks appear genes that are also part of the enriched pathways and processes (Table 3A and 3B, Additional file 7). For example many stress response genes including those involved in the GO Biological Process "endoplasmic reticulum unfolded protein response" are represented in the third ranking network (Additional file 1, Figure S4).

Comparison of BaP induced responses between Oikopleura and human cells

In order to assess conservation of responses to BaP in mammals and Oikopleura, we compared pathways significantly affected by BaP in this experiment and in a gene set from a microarray experiment using human macrophages [37]. Many oxidative stress and apoptosis related pathways were affected in both Oikopleura and human cells (Additional file 1, Figure S5 A and B) suggesting conservation of responses to genotoxic carcinogenic compounds. However, consistent with our genome survey results that showed absence of AhR signaling genes (Figure 2), the AhR signaling pathway was not enriched in Oikopleura, but among the top enriched pathways in human cells (Additional file 1, Figure S5B). Other major differences in the top 10 pathways enriched are related to cell-cycle regulation, which were among the top over-represented in human cells but not in Oikopleura (Additional file 1, Figure S5B). This could be because of species specific differences in cell-cycle control mechanisms, which in Oikopleura involves endocycling, and not necessarily mitosis [38]. Generation of transcription factors centred over-connected hubs in MetaCore (Transcription Regulation algorithm) resulted in a largely similar list of transcription factors for Oikopleura and human BaP-regulated gene sets (Additional file 6, not shown). Not surprisingly, AhR was in the human list of transcription factor networks (not shown), but not in the Oikopleura list, suggesting lack of modulated AhR target genes in the latter. Another notable difference is that BaP appears to provoke stronger endoplasmic reticulum (ER) stress response in Oikopleura as suggested by significant over-representation of Apoptosis and survival_Endoplasmic reticulum stress response pathway in Oikopleura but not human cells (Additional file 1, Figure S5B). It is not clear why Endoplasmic reticulum stress response pathway is responding more in Oikopleura. This pathway can be adaptive, promoting survival or it can lead to induction of apoptosis [39]. However, Apoptosis genes respond equally in both Oikopleura and the human macrophage cells (Additional file 1, Figure S5A). ER stress response can also explain the significant enrichment of Cystic fibrosis disease as one of the processes affected by BaP in Oikopleura but not human cells (Additional file 1, Figure S5A). Some pathways and processes such as human diseases are enriched here because the pathway analyses were done on mouse homologs of differentially regulated Oikopleura genes, but the underlying changes in gene expression or cellular processes can be relevant also in Oikopleura. For example, inspection of the modulated genes that led to the enrichment of Cystic fibrosis disease reveals that they are involved in ER stress response (e.g. HSPs). ER stress response is thought to contribute to mechanisms leading to the development of lung fibrosis in mice [40].

Genes and pathways affected by Clo

Treatment of O. dioica with Clo resulted in differential regulation of 630 genes, with 166 genes up-regulated and 464 genes down-regulated (Additional file 4). Genes involved in the xenobiotic defense system in vertebrates were part of the list, including biotransformation enzyme genes such as CYP2-like-11 gene (CBY24358), conjugating enzyme genes such as GSTα, and some ABC transporters and solute carrier proteins (Additional file 4). Superoxide dismutase and many ABC transporters were also down-regulated by Clo.

The list of genes differentially regulated by Clo was used in various functional analyses using KEGG and MetaCore (GeneGo) as described above for BaP gene set. Clo affected mainly muscle contraction and energy metabolism pathways (Tables 4 and 5, Figure 5). Among the significantly enriched pathways, ECM-receptor interaction and Focal adhesion are also affected by BaP (Table 2). Clo had orchestrated effects on main energy pathways such as Glycolysis/Gluconeogenesis, Oxidative Phosphorylation and Citrate cycle (TCA cycle) (Table 4). While most pathways affected by Clo are presumably conserved between Oikopleura and mouse, some such as the disease pathways (Table 4) may be less relevant for Oikopleura. Remarkably, nearly all the genes affected by Clo in the enriched muscle component and energy pathways were down-regulated (Figure 5, not shown), suggesting functional suppression of energy metabolism and motility. This appears consistent with the reduced motility of the animals observed for the highest concentration of Clo (not shown).

Term	Count	PValue	Benjamini	FDR
ECM-receptor interaction	12	8.28E-05	1.15E-02	0.098
Valine, leucine and isoleucine degradation	9	1.15E-04	7.99E-03	0.135
Cardiac muscle contraction	11	2.31E-04	1.07E-02	0.273
Arginine and proline metabolism	9	3.19E-04	1.11E-02	0.376
Glycolysis/Gluconeogenesis	10	3.74E-04	1.04E-02	0.441
Hypertrophic cardiomyopathy (HCM)	11	4.27E-04	9.92E-03	0.503
Focal adhesion	17	8.40E-04	1.67E-02	0.987
Dilated cardiomyopathy	11	8.86E-04	1.54E-02	1.040
Alzheimer's disease	16	9.98E-04	1.54E-02	1.171
Huntington's disease	16	1.06E-03	1.47E-02	1.238
Oxidative phosphorylation	13	1.21E-03	1.53E-02	1.418
Parkinson's disease	13	1.48E-03	1.71E-02	1.727
Propanoate metabolism	6	2.80E-03	2.98E-02	3.260
Citrate cycle (TCA cycle)	6	3.26E-03	3.21E-02	3.774

Rank	A. Map Folders	pValue
1	Tissue remodeling and wound repair	1.7e-6
2	Energy metabolism and its regulation	1.7e-5
3	Inflammatory response	1.2e-4
4	Cell differentiation	2.6e-4
5	Protein degradation	3.8e-4
6	Calcium signaling	1.0e-3
7	Vasodilation	1.8e-3
8	Lipid Biosynthesis and regulation	4.2e-3
9	Vasoconstriction	4.7e-3
10	Blood clotting	9.8e-3
	B. Networks
1	Muscle contraction	9.4E-15
2	Cytoskeleton_Actin filaments	8.0E-10
3	Development_Skeletal muscle development	2.2E-09
4	Cytoskeleton_Regulation of cytoskeleton rearrangement	1.9E-07
5	Development_Cartilage development	6.5E-05
6	Cytoskeleton_Intermediate filaments	7.9E-05
7	Cell adhesion_Integrin-mediated cell-matrix adhesion	3.2E-04
8	Cytoskeleton_Cytoplasmic microtubules	4.8E-04
9	Cardiac development_BMP_TGF_beta_signaling	5.6E-04
10	Cell adhesion_Cell-matrix interactions	8.4E-04

To see interaction among genes differentially regulated by Clo, the network building tool in Metacore was used (with Analyze algorithm). This resulted in a list of ranked interaction networks (Additional file 8). Consistent with the functional ontology top list (Tables 4 and 5), the first five networks were enriched in genes involved in oxidative phosphorylation and muscle contraction (Additional file 1, Figure S6 A and B).

Comparison of BaP and Clo affected genes and pathways

The major pathways and processes specifically affected by BaP are related to its known genotoxic oxidative stress effects such as apoptosis, DNA-damage response, mitogenic signaling and response to un-folded proteins (Additional file 1, Figure S5 A and B). Pathways and processes specifically affected by Clo are energy metabolism and its regulation, muscle contraction, actin filaments, skeletal muscle development and regulation of cytoskeleton rearrangement (Additional file 1, Figure S7 A and B). Enriched pathways and processes significantly affected by both BaP and Clo are tissue remodelling and wound repair, inflammatory response, cell differentiation and cell adhesion (Additional file 1, Figure S7 A and B).

Comparison of GO Cellular Component (CC) in MetaCore generated different top lists of localizations for BaP and Clo-regulated gene products, reflecting differences in enriched pathways (Additional file 1, Figure S7C). The only localization that is significantly enriched by BaP alone is Golgi apparatus (Additional file 1, Figure S7C). Several localization GO terms were significantly enriched by Clo but not BaP. These include muscle components (e.g. contractile fibre and myofibril) and organelles parts (e.g. mitochondrial matrix), consistent with top enriched muscle contraction and energy metabolism pathways, respectively (Tables 4 and 5, Figure 5, Additional file 1, Figure S7 A-C). The two compounds also share top list of GO localizations, mainly extracellular region (Additional file 1, Figure S7C), because both affected many extracellular proteins, possibly including yet un-characterized structural proteins of the house which is composed of more than 20 glycoproteins proteins [6, 7]. Although not included in pathway and GO analyses, the extracellular house component proteins (oikosins) of Oikopleura are among the genes down-regulated by both compounds (Additional file 1, Table S2). GO analysis also allocated 35% of 125 genes down-regulated by both compounds (Additional file 1, Table S3A) to the GO-CC term extracellular region and 59% of them are annotated as glycoproteins. Only 24 genes were up-regulated by both chemicals (Additional file 1, Table S3B). Thirteen genes were oppositely regulated by BaP and Clo (Additional file 1, Table S3C).

Validation of tiling array results by qPCR

Fourteen genes differentially regulated in tiling arrays were selected for validation using qPCR (Additional file 1, Table S4). The results showed good agreement between the two expression assays (Pearson's correlation coefficient, r = 0.93) (Figure 6). The directions of changes in log transformed fold-changes in expression (treated/control) using the two methods were consistent for all genes with the exception of PA2GD in Clo-treated samples and CYP2R1 in BaP-treated samples (Figure 6). Good correlation between the two methods was achieved because in the tiling array design, each gene in the genome was represented by several overlapping probes along its entire length (see Additional file 1, Figure S3) resulting in a robust statistical quantification of expression levels as described in methods section.

Identification and sequence analysis of Oikopleura genes

Oikopleura genes were identified by both Hidden Markov Model (HMM) searches (HMMER 3.0, http://hmmer.org/) in predicted O. dioica protein database in NCBI using respective PFAM models of conserved domains, and by BLAST searches in the Oikopleura genome (sequenced to 14X coverage) and EST (about 300,000) databases. Mus musculus Swiss-Prot CYPs retrieved from NCBI database were used in BLAST searches and Oikopleura CYP hits were then re-BLASTed to identify all conserved CYPs in the genome. Then we compared with CYPs identified by profile Hidden Markov Model (HMM) searches and identified CYP gene models were manually curated. PAS domain proteins were identified by a combination of HMM searches and BLAST searches using Mus musculus, Danio reiro, Caenorhabditis elegans, Branchiostoma floridae, Ciona intestinalis, Strongylocentrotus purpuratus and Nematostella vectensis proteins retrieved from the NCBI database.

HMM searches and BLAST searches (using mouse Swiss-Prot sequences retrieved from the NCBI database) in O. dioica genome and EST databases were also used to identify transcription factors involved in Nrf2 pathway and other putative oxidative stress response signaling transcription factors (Additional file 1, Figure S2). Protein sequences for other organisms were retrieved from the NCBI database.

Phylogenetic tree analyses were performed on the Phylogeny.fr platform http://www.phylogeny.fr[53]. Multiple alignments were performed using MUSCLE. Alignments were manually edited and PhyML and TreeDyn programs were used to construct and visualize phylogenetic trees, respectively.

Chemical exposure

Preliminary experiments were performed to estimate concentration ranges and exposure times that induced gene expression changes (as determined using qPCR and RNA blot assays of few genes) without causing excess mortality using clofibrate (Clo) (Sigma-Aldrich, St. Louis, MO) and benzo[a]pyrene (BaP) (Sigma-Aldrich). In these experiments 0.2 μM BaP and 5 μM Clo were found to be optimal concentrations. The 24 hour LC50 value of Clo was found to be approximately 14 μM. BaP has very low solubility in seawater and we did not determine its LC50. About 130 four-days-old sexually immature animals (O. dioica) were placed in 1L glass beakers with filtered and UV-light treated seawater, fed with algal culture and kept in suspension by rotation of a paddle connected to an electric motor [8]. The compounds were dissolved in 1 ml Dimethyl sulfoxide (DMSO), diluted in 50 ml seawater, vigorously vortexed and added. DMSO alone (1 ml) was used as a control. Final concentrations of Clo in seawater were 1 μM (0.243 mg/L) and 5 μM (1215 mg/L). Because of very low solubility of BaP, the actual concentrations in the seawater were not known. BaP was added at 0.2 μM (50.4 μg/L) and 1 μM (252 μg/L) nominal concentrations. There were no large and visible precipitates in the 0.2 μM BaP concentration, although invisible precipitates should form. Precipitates formed in the 1 μM concentration of BaP that formed larger aggregates possibly due to enhanced nucleation. The aggregation appeared to have resulted in lower up-take of BaP by the animals in the 1 μM concentration than in the lower (0.2 μM) concentration. Thus, the animals in the higher BaP concentration were swimming normal and less stressed than animals exposed to the lower (0.2 μM) concentration. Well dispersed tiny precipitates in the 0.2 μM concentration appear to have been taken up by the animals which filter microscopic food particles and ingest. The animals were kept at room temperature and harvested after 10 hrs of exposure by picking using inverted plastic pipettes, rinsed in filtered seawater and centrifuged 3 min at 12,000 rpm. Then the pellet was frozen in liquid nitrogen and stored at -80°C until RNA extraction.

RNA extraction and dscDNA synthesis

Total RNA was isolated from the frozen animals using the RNeasy Mini Kit according to manufacturer's protocols (QIAGEN, Hilden, Germany). RNA concentration and quality was assessed using NanoDrop ND-1000 (NanoDrop Technologies, Wilmington, DE), and agarose gel electrophoresis. For each sample, 5 μg total RNA was converted to dscDNA using SuperScript Double-Stranded cDNA Synthesis Kit (Invitrogen, Carlsbad, CA), and dscDNA samples were submitted for microarray analysis.

Whole genome tiling arrays

Custom made O. dioica whole genome tiling arrays manufactured by Roche NimbleGen (NimbleGen, Madison, WI) for the Sars International Centre for Marine Molecular Biology were used for gene expression analysis using a competitive hybridization strategy. The oligonucleotide features (represented by over 2 million probe sequences, with a probe size range between 50-75 nucleotides) cover the entire 70 Mb genome with average overlap between adjacent probes of approximately 30 bases. Thus, each gene in the genome was represented by numerous overlapping probes along its entire length. Labeling (using Cy3-coupled random nonamers for treated samples and Cy5-coupled random nonamers for DMSO control samples), and two-color hybridization (using 15 μg of each labeled sample per array) and scanning using an Axon 4000B scanner (Axon Instruments Inc., Union City, CA) were performed at the University of Iowa using standard conditions described in the Gene Expression User Guide (Roche NimbleGen).

Tiling array data analysis

Array data generated from each group (including controls hybridized to the same array) were quantile normalized separately. Background subtraction was performed by subtracting the top 2% of random probe (negative control) intensity values from individual O. dioica probe intensity values. These random probes (which represent 13,936 of the 2,173,626 features on the array) were included on the array to help define background probe intensity levels. Only probes with positive intensity values after background subtraction were included in the analysis. For a gene to be considered as expressed, at least 50 percent of the probes covering exons of the gene had to have a positive intensity. A gene was assigned with a median intensity of the positive probes covering it. T-test was used to test whether there is any significant difference (alpha level = 0.05) between means of two groups (treated versus DMSO control). MIAME complaint array data has been deposited in NCBI's Gene Expression Omnibus (GEO) database (GEO accession GSE33818).

Genes with sum of the normalized median intensities (DMSO control and treated) lower than 3000 were further removed from the analysis because the signals were found to be unreliable after inspection of expression profiles using SignalMap browser (NimbleGen). To identify differentially regulated genes, the normalized median intensities were used to calculate expression ratios (treated/control) for each predicted gene. Genes with significantly different intensity levels (threshold p-value = 0.05, t-test) and with at least 1.5 fold expression changes were considered differentially regulated between treated and control groups. Applying the cut-off values, largely the 0.2 μM BaP and 5.0 μM Clo concentration samples were used for analysis of significant differential expression. The 1.0 μM BaP and 1.0 μM Clo were not optimal concentrations in inducing differential gene expressions. However, largely overlapping sets of differentially regulated genes were detected for the two concentrations of each compound with "dose-response" trends in the majority of cases. Less than 10% of the significantly differentially regulated genes were specifically regulated by 1.0 μM but not 0.2 μM BaP or 1.0 μM but not 5 μM Clo, and these were also included in the list of differentially regulated genes applying more stringent p-values (threshold p-value = 0.01, t-test). For genes for which the predicted protein was annotated with the best BLAST hit (see below), the analysis resulted in 762 (336 up-regulated and 426 down-regulated) and 630 (166 up-regulated and 464 down-regulated) genes differentially regulated by BaP and Clo, respectively. In addition, 425 and 446 genes that were not annotated (BLAST e-value ≥ 0.001) were differentially regulated by BaP and Clo, respectively.

Annotation and functional classification of differentially regulated genes

Predicted protein sequences of 17,112 genes in Oikopleura Genome Browser V3 https://www.genoscope.cns.fr/secure-nda/Oikopleura/cgi-bin/gbrowse/Oikopleura/ were BLAST searched against the mouse proteome in NCBI database (NCBI BLASTP) and each predicted protein was named according to the best mouse hit (E-value < 0.001). Thus, all subsequent pathway and functional analyses were performed on the mouse homologs (the best mouse BLASTP hits) of the differentially regulated O. dioica genes using mouse annotation databases. This was done because most tools available for functional and pathway analyses are developed mainly for mammalian model organisms and do not allow direct inputs of Oikopleura genes. The differentially regulated genes were analyzed using Gene Ontology (GO), KEGG and other tools in DAVID Bioinformatics Resources [34, 35]. The mouse protein_gi_accessions for 762 BaP-differentially regulated genes (Additional file 3) and 630 Clo-differentially regulated genes (Additional file 4) were submitted to DAVID database http://david.abcc.ncifcrf.gov/ and analyzed using mouse proteome background and default settings. In addition, the gene lists were subjected to more detailed pathway, process and network analyses using the MetaCore knowledge database and software suite (GeneGo, St. Joseph, MI) [36]. Predicted proteins with insufficient homologies to mouse proteins (BLASTP E-value ≥ 0.001) that accounted for 36% of BaP regulated genes (Additional file 9) and 41% of Clo regulated genes (Additional file 10) were not included in pathway analyses.

Quantitative real-time PCR (qPCR)

qPCR was performed for validation of differential expressions in the microarray experiments. A total of 14 genes found to be differentially regulated using microarray by one or both compounds (0.2 μM BaP and 5 μM Clo) were analyzed. One of each primer pair (Additional file 1, Table S4) was designed to span exon-exon junctions to avoid amplification of genomic DNA. The amplicon sizes were 80-150 bp. RNA polymerase B11a (RPB11a) (Additional file 1, Table S4) was used as a house-keeping gene for normalization (reference gene). In preliminary experiments using the same RNA samples, comparison of RPB11a with beta actin showed both are good housekeeping reference genes. The same RNA samples analyzed by microarrays were used. For each sample, total RNA (1.0 μg) was reverse-transcribed using SuperScript III First-Strand Synthesis System for RT-PCR in 20 μL reaction as described in the manufacturer's protocols (Invitrogen). Then the reaction was diluted 1:10, and 5 μL of the cDNA was used in 20 μL amplification reaction using FastStart SYBR Green Master Mix according manufacturer's protocols (Roche, Basel, Switzerland), and qPCR was performed using LightCycler 480 Real-Time PCR System (Roche). The thermal cycling conditions were as follows: an initial denaturation of 95°C for 10 min, followed by 45 cycles of denaturation at 95°C for 10 s and an annealing temperature of 55°C for 20 s and elongation at 72°C for 30 s. Negative controls with no reverse transcriptase enzyme were run for each primer pair and each sample was amplified in duplicate. Post-amplification melting curve analysis was performed to check specificity of products. The PCR products were also analyzed by agarose gel electrophoresis to verify the amplification of a single product of the right size. The relative quantification method [54] was used to calculate gene expression and expression level for each gene was expressed as fold-change relative to the DMSO control.