The transcriptome-wide effects of exposure to a pyrethroid pesticide on the Critically Endangered delta smelt Hypomesus transpacificus

Table S1. Gene descriptions, gene symbols and relative log2 fold changes for the 69 microarray features that were significantly different between the lowest exposure concentration (0.69 μg l) and the control group (ANOVA; q < 0.05) after 96 h of exposure to permethrin. Features that were significantly different from the control group are in bold. The order of the genes in the table is the same as in the heat map.


INTRODUCTION
Assessments of contaminant impacts generally rely on the use of well-studied model test species that have been approved by regulatory agencies for use in toxicity testing (e.g. fathead minnow Pimephales promelas and rainbow trout Oncorhynchus mykiss in fresh water, sheepshead minnow Cyprinodon variegatus in salt water). However, the responses in model test species may be different than those species of conservation concern depending on the contaminant or the species (Dwyer et al. 2005). While the use of surrogate species rather than endangered species for toxicological testing is technically preferable and sometimes necessary, the ability to assess contaminant impacts on endangered relative to surrogate species requires a comparative evaluation of sensitivities to a particular contaminant or classes of contaminants. Comparative sensitivity studies conducted with endangered species and model test species have been used to determine sensitivity factors that could be applied for the protection of species listed under the United States Endangered Species Act (ESA; e.g. Dwyer et al. 2005). These sensitivity factors are useful for determining potential relative effects of exposure for certain species that are not amenable to or available for laboratory study; however, when possible, a direct approach of conducting toxicological studies on cultured species of conservation concern is potentially most informative (Dwyer et al. 2005). Comparative sensitivity assessments have mainly focused on acute lethal exposures, examining lowest observed effect concentrations (LOEC) or no observed effect concentrations (NOEC), which are generally higher than the concentrations aquatic organisms experience in the environment and may not be as sensitive to sub-lethal exposure as changes in physiology, behavior and gene expression. Many biological processes such as tissue differentiation, cell proliferation and growth may be impacted by sub-lethal contaminant exposure (Mazurais et al. 2011), and therefore potentially provide a more sensitive response to low exposure concentrations than standard test endpoints and may ultimately be more predictive of longer-term adverse effects.
The delta smelt Hypomesus transpacificus is a pela gic fish species endemic to the Sacramento−San Joaquin Rivers and the San Francisco Estuary, CA, USA. This species has experienced a decline in abundance since the 1980s (Moyle et al. 1992, Bennett 2005, with a significant step decline in 2000 (Feyrer et al. 2007). The delta smelt was classified as threatened under the Federal and State ESA in 1993, listed as endangered under the California ESA in 2010 (CDFW 2014), and classified as Critically Endangered by the International Union for Conservation of Nature in 2014 (IUCN; www.iucnredlist.org/details/10722/0). Factors contributing to the species' decline include habitat de gradation and loss, competition with non-native species, decreased food availability due to changes in community structure, and alterations in the physiochemical properties of the system (reviewed in Bennett 2005. Harmful effects are also likely caused by contaminants entering the system from wastewater treatment plant effluent and surface runoff (Kuivila & Foe 1995, Thompson et al. 2000, Werner et al. 2000, Kuivila & Moon 2004, Brooks et al. 2012. Among the numerous contaminants of concern in surface waters throughout North America are pyrethroid insecticides. Permethrin, a type of pyrethroid pesticide, was originally registered for use by the United States Environmental Protection Agency (US EPA) in 1979, and re-registered in 2006(US EPA 2009. It is broadly used in mosquito abatement programs and for pest control on feed crops, livestock, pets, and on urban wooden structures and buildings. It is also used as a pharmaceutical for the treatment of head lice and scabies on humans. Consequently, there are multiple potential sources of permethrin entering aquatic systems via surface runoff or wastewater discharge. The half-life of permethrin in the water column ranges between 19 and 27 h, however when adsorbed to sediments, permethrin can persist for longer than 1 yr (Imgrund 2003). Permethrin can be highly toxic to marine, estuarine and freshwater aquatic species. Pyrethroid pesticides disrupt the nervous system in vertebrates through prolonged opening or closing of ion channels, which can result in convulsions, paralysis and eventual mortality (Burr & Ray 2004, Werner & Moran 2008. At lower relative exposure concentrations, sub-lethal neurotoxic effects in fish include impaired swimming ability (Beggel et al. 2010, Goulding et al. 2013) and a reduced ability to avoid predators (Floyd et al. 2008).
Detectable levels of permethrin have been re por ted throughout the Sacramento−San Joaquin watershed (Holmes et al. 2008), with reported concentrations of up to 0.094 µg l −1 in tributaries of the Sacra mento River (Bacey et al. 2004), and in excess of 1 µg l −1 in storm drains that discharge into the Sacramento− San Joaquin Rivers Delta (Zhang 2010). As a pelagic fish species, delta smelt would primarily be ex posed to permethrin in the water column, and because of the relatively short half-life, the permethrin exposures that delta smelt experience in the environment may be acute. Detections of permethrin in California waterways generally occur during the rainy season (~November to April; Weston & Lydy 2010), which also coincides with the timing of delta smelt returning to their freshwater habitat to spawn (~December to May; Bennett 2005). Therefore, it is important to determine the risks posed to delta smelt in their early developmental stages.
With advances in the development of molecular tools and bioinformatics resources, alterations in cellular processes resulting from exposure to contam -inants can be relatively rapidly determined. Highthroughput gene expression analyses, such as RNA se quencing and microarray approaches, allow transcriptome-wide responses to be evaluated. These effect-based responses provide a greater understanding of the impacts of exposure to contaminants, which can potentially be used towards predicting adverse outcomes at higher levels of biological organi zation (Ankley et al. 2010, Connon et al. 2012a. Prior studies on the delta smelt have reported the development and application of cDNA microarrays for the assessment of contaminant exposures (Connon et al. 2009, 2011a,b, Hasenbein et al. 2014); however, we recently sequenced the delta smelt trans criptome and developed a new oligonucleotide microarray. Here, we present the use of this tool towards determining the sublethal, transcriptome-wide effects of acute exposure to permethrin on an early life-stage of the Critically Endangered delta smelt.

Permethrin exposures
Delta smelt (~43 d post hatch), were obtained from the University of California, Davis Fish Conservation and Culture Laboratory facility in Byron, CA. Fish were given 24 h to recover from transport and handling stress in aerated, filtered (1 µm) culture facility water. A total of ten ~44 d old fish were randomly transferred into one of 32 glass vessels (four 3.5 l replicate vessels containing aerated, filtered culture facility water for each exposure treatment, a control and a solvent control). A stock solution of permethrin (CAS No. 52645-53-1; Chem Service) was prepared in methanol and spiked into treatment solutions (culture facility water) to achieve nominal concentrations of 1.25, 2.5, 5.0, 10.0, 20.0 and 40.0 µg l −1 . Treatment solutions were spiked at test initiation, and again at each 24 h renewal when 80% of the water was exchanged. Fish were exposed for 96 h. Mortality was recorded daily when dead organisms and debris were removed from the test vessels. At the end of the exposure, all surviving fish (from the solvent control and the 1.25, 2.5, 5.0 and 10.0 µg l −1 [nominal] treatments) were sacrificed with a lethal dose of tricaine methanesulfonate (MS-222), frozen in liquid nitrogen and stored at −80°C for subsequent assessments. Throughout the 24 h recovery and the 96 h exposures, fish were fed 3 times daily with live Artemia spp. nauplii.
Comprehensive Environmental Toxicity Information System (CETIS) by Tidepool Scientific Software was used to calculate lethal effect concentrations (i.e. LC 50, the concentration that was lethal for 50% of the individuals in the experiment), LOEC and NOEC. Mortality patterns were also analyzed using a Kruskal-Wallis (KW) test. Water quality parameters were measured at test initiation, before and after each water renewal and at test termination for mean (± SD) electric conductivity (743.8 ± 4.7 µS cm −1 ), temperature (15.4 ± 0.1°C), dissolved oxygen concentration (9.35 ± 0.13 mg l −1 ), pH (8.03 ± 0.03) and total ammonium concentration (0.017 ± 0.041 mg l −1 ). Water hardness (164 mg l −1 ) and alkalinity (92 mg l −1 ) were measured at test initiation.

Analytical chemistry
Sub-samples of permethrin treatments and controls (1 l) were preserved with 15 ml of 99.9% di chloromethane (Fisher Scientific) at test initiation and shipped on ice overnight to the California Department of Fish and Wildlife Water Pollution Laboratory in Rancho Cordova, CA. Extractions were performed within 24 h of arrival and analyzed using gas chromato graphy with mass spectrometry and ion-trap de tection. Permethrin recovery was approximately 50% of nominal concentrations, which is within the normal range of recovery for permethrin (see Wheelock et al. 2005). Nominal to corresponding measured isomer concentrations were 2.5 = 1.37 (0.52 cis, 0.85 trans) µg l −1 ; 5.0 = 2.56 (0.99 cis, 1.57 trans) µg l −1 ; 10.0 = 4.84 (1.89 cis, 2.95 trans) µg l −1 ; 20.0 = 12.88 (4.87 cis, 8.01 trans) µg l −1 ; and 40.0 = 24.94 (9.54 cis, 15.40 trans) µg l −1 . The lowest nominal concentration (1.25 µg l −1 ) was not measured and was thus estimated to be 0.69 (0.27 cis, 0.43 trans) µg l −1 based on the higher concentration measurements. The ratio of trans to cis isomers was consistent among treatments, ranging from 1.56 to 1.63, which is important because the trans isomer of permethrin has been shown to be more toxic than the cis isomer in rainbow trout (Glickman et al. 1981).

Microarray development
We developed an oligonucleotide microarray (Agilent Technologies) using transcriptome sequencing data from delta smelt that were exposed to a number of stressors (i.e. temperature, salinity, hypoxia, turbidity, pH, exercise, olfactory, ammonia, copper, permethrin, esfenvalerate) that was designed to capture a diversity of mRNA transcripts related to the stress response. Total RNA was extracted from tissue samples (i.e. brain, heart, gill, kidney, liver, spleen, skin, muscle) or frozen whole organisms using a combination of TRIzol Reagent (Life Technologies) and RNeasy kits (Qiagen), including on-column DNAse digestion, as per manufacturer's guidelines. RNA sample integrity was verified using a 2100 Bioanalyzer (Agilent Technologies) that indicated high quality RNA (RNA integrity number, RIN values > 9). Samples were pooled to represent an equal proportion of RNA in tissue (adults) and whole organisms (larval and juvenile fish) to normalize and maximize the identification of rare transcripts. Transcriptome sequencing was conducted at the Center for Geno mics and Bioinformatics (CGB) at Indiana University, Bloomington, IN. Two sequencing libraries were constructed from 5.0 µg of pooled total RNA normalizing transcript representation and optimizing transcriptome sequence read lengths. Libraries were se quenced using a Roche 454 GS-FLX Titanium system. Clean reads were assembled into contigs using Newbler software (Roche) and annotated by the bioinformatics core at CGB. Custom 4 × 44K oligo nucleo tide microarrays were designed utilizing Agilent specific software (https:// earray. chem. agilent. com/ earray). The array consisted of a total of 39 652 microarray features (i.e. spots on the microarray), with a total of 18 409 duplicate delta smelt features representing 12 595 unique delta smelt genes, plus non-specific control duplicate probes. There were a higher num ber of delta smelt features than unique genes because some of the available sequence fragments represented different gene regions, therefore there were different microarray features that represent different segments of the same gene on the microarray.

Tissue processing and microarray analysis
Gene expression assessments were conducted on individual fish from the solvent control and the 0.69, 1.37, 2.56 and 4.84 µg l −1 treatments. Mortality was > 90% at the 12.88 and 24.94 µg l −1 test concentrations, and therefore these treatments were not in cluded in the gene expression analyses. Due to the small size of the larval fish, individual whole fish were homogenized using a pestle in 1 ml of Trizol Reagent following manufacturer's protocols, followed by RNA Cleanup using RNeasy columns (Qiagen) with on-column DNAse digestion. Total RNA con centration and purity (260/280, 260/230) were as ses sed using a Nanodrop spectrophotometer (Thermo Fisher Scientific), and RNA integrity was determined visually by running samples on an agar ose gel. Samples were stored at −80°C until used in the gene expression analyses.
In total, 200 ng of total RNA was amplified and labeled with Cy3 using the Low Input Quick Amp Labeling kit (Agilent Technologies). A total of 20 single color array hybridizations (n = 4 ind. treatment −1 ) were performed with custom designed 4 × 44K oligonucleotide arrays using Gene Expression Hybridization Kits (Agilent Technologies). To minimize technical artefacts, labeling reactions were performed simultaneously, and individuals from each treatment were randomized between slides. Prior to hybridization, 1.65 µg of dye-labeled cRNA sample was combined with 2.2 µl of 25× fragmentation buffer, 11 µl of 10× Gene Expression Blocking Agent and nucleasefree water to bring the final volume to 55 µl. The fragmentation mix was incubated at 60°C for 30 min, cooled on ice for 1 min, and the reaction was stopped by adding 55 µl of 2× Hi-RPM Hybri di zation Buffer. Samples were centrifuged for 1 min, placed on ice, 100 µl of the mix was loaded onto gasket slides, and the microarray slides were placed on top of the gasket slides. Each slide and gasket slide combination was placed in a single hybridization chamber and hybridized to the arrays for 17 h at 65°C followed by a wash with Gene Ex pression Wash Buffer 1 at room temperature and a wash with Gene Expression Wash Buffer 2 heated to 37°C (Agilent Technologies). Fluorescent images were scanned using an Axon GenePix 4000B Scanner and the analysis software GenePix Pro (Molecular Devices). The images were quantified using Feature Extraction software v.11. 5.1.1 (Agilent Tech no logies). The microarray data have been de posited in the Gene Expression Omnibus (www. ncbi. nlm. nih. gov/geo/) with the ac ces sion number GSE67521.
Normalization of the microarray data and statistical analysis were performed in GeneSpring v.12.6 (Agilent Technologies). Features that were less than 2.5 times the intensity of the background for that array were filtered out prior to normalization. Data was quantile normalized and log 2 transformed prior to statistical analysis. Differences in expression of features between treatment groups were assessed by ANOVA at q < 0.05 (a false discovery rate corrected p-value) followed by a Student Newman-Keuls (SNK) post-hoc test. It is important to note that there are 2 different overlapping microarray features for the same gene sequence that were treated as separate genes in the statistical analyses. Therefore, multiple copies of each gene fragment that were significant in our analyses showing the same directional change in expression were used to confirm the re sponse pat-terns and distinguish true patterns from random noise in the data. Functional analysis was performed using Blast2GO (Conesa et al. 2005). Analysis was performed considering all gene ontology (GO) categories; however, to reduce the large number of significantly enriched GO terms, only biological process terms are presented here as they sufficiently represented the overall functional analysis results. Functional groups were considered significantly enriched in gene lists using Fisher's exact tests at a false discovery rate (FDR) < 0.05. Redundancies in the GO term lists were reduced using ReVigo (Supek et al. 2011).

cDNA synthesis and quantitative PCR
One µg of total RNA was used to synthesize cDNA in 96-well plates. For each sample, nuclease-free water was added to bring the total volume up to 12 µl, followed by the addition of 2 µl of gDNA wipeout buffer and an incubation at 42°C for 2 min. For each cDNA synthesis reaction, 4 µl 5× Quantiscript RT Buffer, 1 µl RT Primer Mix, and 1 µl Quantiscript RT (Qiagen) were added to each sample for a total reaction volume of 20 µl. Reactions were carried out at 42°C for 30 min, followed by 95°C for 3 min and a hold at 4°C. The cDNA was stored at −20°C until quantitative PCR (qPCR) analysis.
We designed qPCR assays for 22 genes of interest (Table 1) to validate some of the microarray expression patterns and to verify the non-monotonic response patterns we observed in the microarray results due to exposure to permethrin (n = 10 to 13 fish treatment −1 , which included the same fish used in the microarray analysis). A 1:5 dilution of cDNA was used as a template in the qPCR reactions. Assays were performed on an automated fluorometer (ABI HT 7900 A FAST Sequence Detection System; Applied Biosystems) in 384-well plates using 12 µl reaction volumes composed of 6 µl of Maxima Probe/Rox qPCR Master Mix (Thermo Fisher Scientific), 0.6 µl primer-probe mix (forward primer, reverse primer, the appropriate Roche fluorescent probe and nuclease-free water), 0.4 µl of nuclease-free water and 5 µl of diluted cDNA. Cycling conditions were 2 min at 50°C, 10 min at 95°C, 40 cycles of 15 s at 95°C and 60 s at 60°C, and amplifications were analyzed using SDS v.2.4 software (Applied Biosystems). Relative expression of target genes was determined using the 2 −∆∆Ct method and are presented as relative fold change (Livak & Schmittgen 2001). Target gene expres sion was normalized to 3 reference genes (ACTB, GAPDH, and RPL7). Of the reference genes, only RPL7 was newly designed for these qPCR assessments; ACTB and GAPDH had been de signed from previous studies (Connon et al. 2011a, Hasenbein et al. 2014). Stability of the reference genes across experi mental treatments was initially as ses sed using the program GeNorm (Vandesompele et al. 2002). Differences in expression levels were as sessed using a 1-factor ANOVA followed by a SNK post-hoc test (p < 0.05). Where the data did not fit the assumptions of an ANOVA, a KW test was used followed by a Dunn's post-hoc test (p < 0.05).

Acute toxicity
The 96 h LC 50 for measured permethrin concentrations was 4.07 µg l −1 (95% CI = 3.43−5.50 µg l −1 ), with a LOEC of 2.56 µg l −1 (87.5% survival) and a NOEC of 1.37 µg l −1 . The LC 10 was 2.38 µg l −1 (95% CI = 1.78− 2.83 µg l −1 ). Survival was significantly reduced at higher exposure concentrations, with 35.0, 3.0 and 0% survival at 4.84, 12.88 and 24.94 µg l −1 , respectively (p < 0.001; Fig. 1). The LC 50 for delta smelt was compared with those of different marine and freshwater model test species and other species of conservation concern (Table 2). The LC 50 of the freshwater model test species (fathead minnow and the bluegill sunfish, Lepomis macrochirus), were 2.3 to 3.8 and 1.9 times higher, respectively, whereas the common marine test species (sheepshead minnow) had an LC 50 that was 4.2 times higher than delta smelt. Estuarine model test species (inland silverside Menidia beryllina and Japanese medaka Oryzias latipes; a 48 h exposure) had LC 50 levels that were 6.8 and 2.7 times higher, respectively, than delta smelt. LC 50 estimates for rainbow trout varied due to the life-stage tested and the exposure temperature; however, the LC 50 estimate on the youngest life-stage was 0.81 times that of delta smelt. Juvenile rainbow trout at comparable exposure temperatures as in the present study had LC 50 estimates that were 1.58 to 1.72 times higher than delta smelt. Similar to the larval rainbow trout, larval zebrafish Danio rerio had an LC 50 that was 0.61 times that of delta smelt; however, the rainbow trout and zebrafish tests were both conducted in fresh water as opposed to the low salinity conditions in the present study. When comparing the LC 50 estimates for other species of conservation concern that are available in the literature, the LC 50 of delta smelt was greater than 5 species (ranging from 1.22 to 2.58 times greater), less than 3 species (ranging from 1.46 to 6 times greater than the delta smelt), and effectively the same as the Cape Fear shiner Notropis mekistocholas (LC 50 = 4.16 µg l −1 ).

Microarray assessments
There were 3342 differentially expressed microarray features at q < 0.05. Because of the large number of significant features, we focused on those that were significantly different from the control group after the post-hoc SNK test, which reduced the number of significant features to 2278 (Fig. 2). There were 69 microarray features (55 separate genes) that differed between the lowest exposure concentration (i.e. 0.69 µg l −1 ) and the control group (23 features unique to that comparison; Fig. 3, see Table S1 in the Supplement at www.int-res. com/articles/ suppl/ n028 p043 _ supp.pdf for gene descriptions and names). This is in sharp contrast to the 1993 features that were significantly different between the fish exposed to 4.84 µg l −1 permethrin and the control group. A complete summary of the number of microarray features that were significantly different in all pairwise comparisons between each treatment group is available in Table S2 in the Supplement. Interestingly, for the 1.37, 2.56 and 4.84 µg l −1 exposure groups, there were more microarray features that were sig ni ficantly different from the 0.69 group than were significant different from the control group.
Many of the 69 significant features in the 0.69 µg l −1 treatment group had opposite expression patterns (i.e. directional change) from those found at the higher exposure concentrations (Table S1). For example, genes involved in proteosomal breakdown (PSMA1; CAND1; PSMD6; OCIAD1; NDFIP2), protein synthesis (ELOF1; MTERFD3; RNF10), tyrosine catabolism (HPD), and a gene with an unknown function were downregulated at 0.69 µg l −1 ; however, they were significantly upregulated in the 4.84 µg l −1 group relative to the controls. Conversely, the genes IGSF10 (involved in osteogenesis) and BDH2 (iron homeostasis and potentially associated with a stress response) were significantly upregulated in the 0.69 µg l −1 group, but were significantly downregulated in the 4.84 µg l −1 group relative to controls.
Functional analyses were performed on the significant gene lists between the control group and each exposure concentration to identify the main trends in the data (Table 3). There were genes involved in several biological processes associated with proteosomal protein breakdown and immune responses that were significantly differentially expressed in the 0.69 and 4.84 µg l −1 exposure groups. Interestingly, many of these genes demonstrated opposite response patterns between the 0.69 and 4.84 µg l −1 groups. Other GO categories, which were enriched in the significant gene lists, were associated with cell signaling pathways, metabolic processes, cell proliferation and an unfolded protein response.
Fourteen microarray features (12 different genes) were statistically significant and had the same directional change in expression at each exposure concentration relative to the control. Multiple copies of ATG101 (involved in autophagy and protein binding) and 1 copy of SPAST (involved in cytoskeleton and cell division processes) were significantly upregulated at every exposure concentration. Conversely, genes involved in protein synthesis (MYCT1; TRA2A; MRPS16), binding (PCDH2AB6, cell adhesion and calcium binding; BAP18, DNA binding; METTL21A, protein binding and general stress), metabolism (NDUFS8), cell signaling (SEZ6), and in an immune response (DB-1) were all significantly downregulated in all the exposure treatments relative to the control. In addition to these 12 common genes, the gene ZNF106A, involved in metal ion binding, was downregulated at all exposure concentrations except at 2.56 µg l −1 . Lastly, SLC 22A2, which is involved in ion transport and uptake of organic compounds from blood circulation, was upregulated at every exposure concentration except at 2.56 µg l −1 .  Table 2. Comparisons of the permethrin LC 50 (concentration that was lethal for 50% of the individuals in an exposure treatment) estimates on freshwater and marine model test species and species of conservation concern in the early life-stages. FW: fresh water; HS: high saline ponds; SW: salt water; YOY: young of year; IUCN: International Union for Conservation of Nature (IUCN listing as per www.iucnredlist.org/); na: not available qPCR assessments Several qPCR assays were designed to examine the effects of permethrin exposure on specific genes associated with the functional categories that were enriched in the microarray data set, or that are particu-larly relevant to developing fish (Fig. 4, see Table 1 for gene names). All responses were significant at p < 0.001 unless indicated otherwise. In all cases, the qualitative patterns matched the responses observed in the microarray data; however, most of the statisti-cally significant responses in the qPCR assays were only detected in the 2.56 and 4.84 µg l −1 treatment groups. For example, the following were all signifi-cantly differentially expressed between the 2.56 and 4.84 µg l -1 groups and the control and the 0.69 µg l −1 group: genes involved in general metabolism (APOA-I-2, cholesterol transport; FABP1, fatty acid binding and transport; IDH1, citric acid cycle), protein break-down (CTSD, lysosomal protein breakdown; PSMA1 and PSMC3, proteasomal protein breakdown; ZFAND2B, associated with proteosomal degradation), protein synthesis (BYSL, cell proliferation and tran-scription regulation; CEBPD, an inducible transcrip-tion factor that regulates many genes involved in an immune response and may be associated with apop-tosis processes in stressed fish, Sleadd & Buckley 2013; TP53INP2, a transcription factor and regulator of autophagy), structural collagens in the extracellular matrix (COL1A2; COL1A1B, p < 0.01) and in a general stress response (HSC71, however is responsive to ex-posure to contaminants such as heavy oil and organochlorines; Deane & Woo 2006, Nakayama et al. 2008) (Fig. 4). These patterns were consistent with those observed in the microarray analysis.  Control 0.69 µg l -1 1.37 µg l -1 2.56 µg l -1 4.84 µg l -1 1.5 -1.5 Fig. 3. Heat map showing the gene response patterns of larval delta smelt Hypomesus transpacificus exposed to different concentrations of permethrin that were significant at q < 0.05 (ANOVA) in the microarray analysis. The 69 genes that were significantly different between the lowest exposure concentration (0.69 µg l −1 ) and the control group are enlarged and labeled (see Table S1 in the Supplement at www. int-res. com/ articles/ suppl/ n028 p043_ supp. pdf for gene descriptions and full names of the 69 genes). Data presented are log 2 normalized expression values (scaled from -1.5 to 1.5); blue represents a downregulation and yellow represents an upregulation  Table 3. Gene ontology (GO) terms associated with the biological process category that were significantly enriched at a false discovery rate (FDR) corrected p < 5E −2 in the significant gene list between each treatment and control group (ANOVA; q < 0.05) of larval delta smelt Hypomesus transpacificus. Because there were many more GO terms that were significantly enriched and unique to the 4.84 µg l −1 vs. control group comparison than in the other comparisons, only those unique to that comparison that were significant at p < 1E −5 are presented here Fig. 4. Changes in gene expression in larval delta smelt Hypomesus transpacificus after 96 h of exposure to permethrin as determined by quantitative PCR. There is a reference line at 1.0 because a negative mean relative fold change is impossible using the 2 −∆∆Ct method for qPCR analysis; therefore, values >1 indicate an upregulation and values <1 indicate a downregulation of a gene. Statistical differences at p < 0.05 from the control and the 0.69 µg l −1 group are indicated by ** and statistical differences from only the 0.69 µg l −1 group are indicated by * Interestingly, many of the genes assayed in the 2.56 and 4.84 µg l −1 groups were significantly different from the lowest exposure concentration (i.e. 0.69 µg l −1 ), but not significantly different from the control group, which was in contrast to the patterns observed in the microarray. These patterns may potentially suggest a non-monotonic response in these genes (Fig. 4). This response pattern was observed in genes associated with apoptosis (CASP3, an executioner caspase involved in the final stages of apoptosis, p < 0.01; CASP9, which activates the executioner caspases, p < 0.01), immune response (BNIP3L, which interacts with viral responsive genes and can induce apoptosis), and protein breakdown (CTRB2, protein catabolism; HINT1, regulates proteasomal degradation in addition to regulating transcription and TP53-mediated apoptosis). Only HSP30, an inducible molecular chaperone that may be upregulated in response to an abundance of ubiquitinated proteins (Young & Heikkila 2010), did not follow a dose-response pattern and was only significantly upregulated in the 4.84 µg l −1 group (p < 0.05) compared with the 0.69 µg l −1 group, which suggests a severe stress response in those fish.

DISCUSSION
We used larval delta smelt to assess the acute lethal and sub-lethal responses to a range of permethrin exposure concentrations in a Critically Endangered fish species. The LC 50 for delta smelt was 4.07 µg l −1 , which was comparable to other species of conservation concern for which a 96 h LC 50 estimate for permethrin was available, and was generally lower than those of model test species. It is important to note that direct comparisons between different species' relative sensitivity to permethrin is confounded by the fact that many experimental factors, such as the temperature or salinity of the exposures, how the exposures were administered (i.e. a pulse exposure or flow-through exposure), how the estimates of permethrin concentrations were determined (i.e. measured or nominal) and the life-stage of the fish tested, can all influence the LC 50 results. This highlights the need to explicitly report these details in the methods to improve the ability to qualitatively compare responses across studies. To truly test relative sensitivity to contaminants, endangered species and the relevant model test species should be tested at the same facility, similar to previous work conducted by Dwyer et al. (2005); however, exposures should also be conducted under similar experimental conditions and at comparable life-stages.

Effects of permethrin on the delta smelt transcriptome
We sequenced the transcriptome of the Critically Endangered delta smelt and developed an oligonucleotide microarray that was used for an effect-based analysis of exposure to permethrin. At permethrin concentrations that resulted in 12.5 and 65.0% mortality (2.56 and 4.84 µg l −1 , respectively), cellular responses were associated with proteosomal degradation and apoptosis that collectively suggest a severe stress condition for those fish. At the lowest tested concentration of 0.69 µg l −1 , there was in creased expres sion of genes associated with metabolic and develop mental processes. There were 12 genes in the microarray analysis that showed the same directional change and relative fold change in expression at all concentrations, suggesting that these genes could be developed as biomarkers of a general response to exposure to permethrin (see Table S1 in the Supplement). There were multiple copies of ATG101, which is involved in stabilizing the expression of ATG13 (autophagy-related protein 13) and protecting ATG13 from proteosomal breakdown to contribute to autophagy in the cell (Mercer et al. 2009), that were upregulated at every exposure concentration. This potentially identifies the role of auto phagy in responding to exposure to permethrin in delta smelt and suggests that ATG101 may be a useful biomarker of this response. In addition to these 12 genes, SLC22A2 (also known as organic cation transporter 2; OCT2), a polyspecific cation transporter found in the kidney and liver that is involved in up take and excretion of various drugs and toxins (Koepsell et al. 2007), was upregulated due to exposure to permethrin. This likely indicates that SLC22A2 is involved in a key biological process for physiologically responding to permethrin in delta smelt. We also detected significant responses in lily-type lectin-2 (upregulated relative to the higher concentrations) and cathep s i nD ,c a s p a s e9a n dc a spase 3 (downregulated relative to the higher concentrations) in the qPCR analysis in the present study. Previous work on the effects of the pyrethroid insecticide esfenvalerate on larval delta smelt found an upre gulation of fish-egg lectin and downregulation of cathep sin S-like, procathepsin B, caspase 3 and caspase 1 at an exposure concentration of 0.0625 µg l −1 (the lowest exposure concentration in that study; Connon et al. 2009). Ad ditionally, Connon et al. (2009) reported that genes involved in iron homeostasis, immune re sponse, and growth and development were differentially ex pres sed in delta smelt ex posed to esfenvalerate, similar to the functional analysis results in the present study. Collectively, the consistent gene expression patterns observed in the 2 studies suggest that these re sponses may be general indicators of exposure to pyrethroid insecticides.
Consistent with the neurotoxic properties of permethrin in fishes, we observed changes in the expression of genes associated with neurological function. We found an upregulation of EDP1, which is often found in the extracellular fluid in the brain of fish and may be associated with learning and neuronal regeneration (Shashoua 1991), similar to the effects of esfenvalerate on ependymin in larval delta smelt (Connon et al. 2009). There was also an upregulation of CHL1, a cell adhesion protein in the extracellular matrix that has a role in nervous system development (Kamiguchi & Lemmon 1997), with exposure to permethrin. Lastly, we detected a downregulation of SEZ6, a gene associated with brain development and cognitive function that may have cell signaling properties (Herbst & Nicklin 1997, Osaki et al. 2011. We used whole-body homogenates to conduct the transcriptome-wide assessments; therefore, there may be additional neurological genes affected by permethrin that would be measured specifically in the brain that we were unable to detect. Collectively, these results suggest that exposure to permethrin may affect brain function in developing delta smelt. We detected differential expression of genes involved in growth and development in the present study. There was a downregulation of genes associated with the structure of the extracellular matrix and bone growth (e.g. COL1A2 and COL1A1B) at 4.84 µg l −1 in both the microarray and qPCR assessments that may indicate an increased likelihood of developmental issues in delta smelt exposed to concentrations of permethrin >1.37 µg l −1 . Permethrin has been shown to cause developmental deformities (e.g. crooked bodies) in embryonic and larval zebrafish at concentrations ≥300 µg l −1 (Yang et al. 2014). We also detected an upregulation of several genes involved in maintaining the cytoskeleton in the microarray results at a much lower exposure concentration (i.e. 0.69 µg l −1 ) than in Yang et al. (2014). Additionally, we detected an upregulation of IGSF10 (which is involved in bone formation) at 0.69 µg l −1 , yet IGSF10 was significantly downregulated at 4.84 µg l −1 . Collectively, these responses may indicate that exposure to permethrin during critical stages of development may lead to longer-term impacts on growth and development for this species.
There was increased expression of genes associated with immune responses at the lowest exposure concentration, yet decreased expression of those genes at the higher concentrations. Immune function has been suggested to be an important indicator of non-lethal adverse impacts on organisms in ecotoxicological research (reviewed in Dunier & Siwicki 1993, Bols et al. 2001, and changes in the expression of genes involved in immune function have been demonstrated in delta smelt exposed to different contaminant stressors (Connon et al. 2009, 2011b, Hasenbein et al. 2014. Permethrin exposure elicited significant immune system and inflammatory re sponses in the present study, which have previously been reported in fish in association with pyrethroid insecticide exposures (Eder et al. 2004, Connon et al. 2009, Jin et al. 2010, Guardiola et al. 2014). Exposure to sub-lethal concentrations of pyrethroids has also been demonstrated to directly affect susceptibility to pathogens in fish (Clifford et al. 2005, Eder et al. 2007, which may be linked to the permethrinin duced changes in the expression of genes associated with an immune response. Decreased expression of immune genes may also be linked with processes associated with mortality, as we detected decreased expression of immune genes at exposure concentrations that resulted in significant mortality, which is consistent with patterns shown in dying Pacific salmon (Jeffries et al. 2012). The fish exposed to 4.84 µg l −1 of permethrin were potentially nearing death (only 35% of the fish survived the 96 h exposure) and this may have contributed to the decreased expression of immune function genes. It is important to note that genes involved in immune function may also be affected by osmoregulatory stress (e.g. Norman et al. 2014). Because the primary mode of action of permethrin on fish is to impair ion channel functioning (Burr & Ray 2004, Werner & Moran 2008, it is possible that the fish were also ex periencing osmoregulatory stress that could have in fluenced the expression of genes in volved in immune function. The gene that showed the largest magnitude of response in the qPCR analysis was HSP30; however this only occurred at the 4.84 µg l −1 exposure concentration. Because HSP30 can be upregulated in response to ubiquitinated proteins (Young & Heik kila 2010), it suggests that the severity of the stress response crossed a critical threshold that activated a HSP30 response at the highest exposure concentration. Previous proteomics work showed that genes involved in ubiquination and proteosomal de gra dation were upregulated in fathead minnow brains exposed to 7.5 µg l −1 nominal concentrations of permethrin for 72 h (Biales et al. 2011). We detected an increase in the expression of genes involved in ubi -quination and proteosomal degradation with ex posure concentrations ranging from 1.37 to 4.84 µg l −1 , which suggests an overall increase in ubi quinated proteins in these fish.

Non-monotonic responses to permethrin
Delta smelt gene expression as determined by the microarray analysis responded in a non-monotonic pattern. In several instances, the directional response in the gene expression at 0.69 µg l −1 was the opposite of those seen at higher exposure concentrations. For example, 10 genes that were significantly downregulated in the 0.69 µg l −1 group were upregulated at 4.84 µg l −1 relative to the control group (see Table S1). Conversely, 3 genes were significantly upregulated in the 0.69 µg l −1 group but were downregulated at 4.84 µg l −1 relative to the control group. This non-monotonic response pattern is supported by the fact that there were more microarray features differentially expressed between all exposure concentrations (i.e. 1.37, 2.56 and 4.84 µg l −1 ) and the 0.69 µg l −1 group than with the control group -a pattern that was verified for some genes in the qPCR assays. Nonmonotonic responses have been previously reported in delta smelt following exposure to esfenvalerate (Connon et al. 2009), another pyrethroid insecticide, where functional genomics analyses indicated altered expression of genes involved in growth and development, apoptosis, and neuromuscular function corresponded with behavioral responses.
Sub-lethal exposure concentrations do not always follow classical dose-response curves, a phenomenon that is important to the field of ecotoxicology (Vandenberg et al. 2012). Non-monotonic, or bi-modal, response patterns are common (Conolly & Lutz 2004) and must be taken into consideration when assessing the sensitivity of a species as they re present the biochemical processes that respond to a specific stimulus. The response peaks in non-monotonic responses, U-shaped or inverted U-shaped, are often near or below the NOECs as seen in the present study, which may demonstrate a mismatch between acute lethal assessments (i.e. LC 50 ) and more sensitive sub-lethal responses that are potentially predictive of long-term impacts. These response types are commonly re ported in ecotoxicology studies that examine endo crinological endpoints (e.g. Ankley et al. 2009, Brander 2013, Viñas & Watson 2013 and may indicate deleterious effect thresholds that im pact longer-term developmental and reproductive outcomes. Such thresholds are also observed in num erous biochemical pathways following exposure to a variety of contaminant classes in fishes (e.g. Dorts et al. 2011a,b, Wang et al. 2013. Non-monotonic re sponses likely arise as a result of a variety of cellular mechanisms depending on the biological pathway im pacted (Conolly & Lutz 2004), for example, through different feedback mechanisms (Ankley et al. 2009) or through interactions with cellular receptors (Vandenberg et al. 2012).

Future research and applications
One of the aims of this study was to enhance the transcriptomic data available for studying the endangered delta smelt, facilitating the development of tools to examine the physiological processes that may be affected by environmental and anthropogenic stressors. The long-term goal of this effort is the application of these tools in the assessment of multiple stressors (such as contaminant exposure and disease) as well as the interaction of contaminants with abiotic stressors (such as temperature and salinity) that impact habitat suitability for this species. Effect-based approaches, such as transcriptomic analyses, that evaluate impacts on biological processes provide complementary tools to toxicity identification evaluation, and when combined with chemical analyses can be used to identify contaminant classes present in the environment. Effect-based tools, such as gene expression profiling, have successfully been used to determine pollutant sources (e.g. Hasenbein et al. 2014) and to assess the health status of wild fish (e.g. Connon et al. 2012b, Jeffries et al. 2014. Permethrin is one of many contaminants that have the potential to impinge on the recovery of delta smelt in the Sacramento−San Joaquin Rivers Delta. Additionally, the interactions of contaminants in mixtures and with environmental parameters are complex. In order to directly and comprehensively evaluate the impact of multiple stressors on this sensitive species, it would be necessary to conduct a multitude of exposures to environmental water samples, either ex situ or in situ, utilizing a suite of molecular, physio logical and reproductive endpoints. To successfully conduct an experiment of this nature, it is necessary to assess the impacts of relevant abiotic factors a priori to differentiate between stressor types or contaminant classes. Within the context of determining the effects of contaminants on species of conservation concern in the wild, it is not necessarily the identification of the contaminant stressor(s) that is important but rather the overall effect that water quality has on the species. However, identification of contaminant sources and classes of contaminants found in a system is important from a regulatory viewpoint. Thus, effect-based assessments conducted over a broad range of sites within the delta smelt's habitat could potentially determine contaminant inputs to the system that may be detrimental to the recovery of the species. The present study builds on the suite of molecular tools that are now available to study this Critically Endangered species, and to aid in conserving and managing this important key indicator species of ecosystem health in the Sacramento−San Joaquin Rivers Delta and San Francisco Estuary.