Salinity- and population-dependent genome regulatory response during osmotic acclimation in the killifish (Fundulus heteroclitus) gill

Estuaries are dynamic environments, where physical parameters such as temperature, hypoxia and salinity can change rapidly, stochastically and periodically, and can vary across small spatial scales. Successful estuarine residents must adjust their physiology to compensate for this environmental variation. The killifish Fundulus heteroclitus is able to adjust its physiology to compensate for extremes in all of these physical parameters (Burnett et al., 2007), but is perhaps most notable for its ability to acclimate to salinities from nearly four times that of seawater down to freshwater (Griffith, 1974). Within Atlantic coast estuarine habitats, F. heteroclitus is the numerically dominant vertebrate (Yozzo and Smith, 1997; Sweeney et al., 1998), and their extreme osmotic plasticity enables them to exploit osmotic niches that span the continuum from marine to freshwater. We apply a comparative approach, together with a physiological genomics experimental design, to explore the functional genomic mechanisms that enable this extreme physiological flexibility and account for subtle differences in hypo-osmotic tolerance.

Though much is known about the molecular effectors of osmotic acclimation, comparatively little is known of the signaling and regulatory networks in fish gills that integrate and transduce environmental cues to initiate the physiological acclimation response (Fiol and Kültz, 2007). Though transcription-independent signal transduction is an important component enabling the osmotic compensatory response (Kültz and Avila, 2001), transcriptomics tools offer a holistic top-down discovery-based approach for exploring systems-level mechanisms involved in this complex physiological response.

Our experimental design includes comparison of two populations that vary in hypo-osmotic tolerance and exposure to two levels of hypo-osmotic stress, and follows physiological and transcriptomic responses throughout an acclimation time course. Though individuals are physiologically flexible, F. heteroclitus populations vary in their ability to acclimate to extreme osmotic challenges. For example, freshwater populations compensate for hypo-osmotic challenge more efficiently than brackish or coastal populations despite long-term rearing in seawater (Whitehead et al., 2011a), and northern coastal populations are more tolerant of hypo-osmotic challenge than southern coastal populations (Scott et al., 2004), indicating evolved intraspecific divergence in osmotic plasticity. We compared a northern coastal and a southern coastal population in a ‘common-garden’ experimental design to distinguish physiological and genomic responses to osmotic challenge that are under strong evolutionary constraint from those that are more prone to evolutionary divergence. We challenged populations with two levels of hypo-osmotic stress. At the less extreme stress level, populations did not differ in physiological response, whereas at the more extreme stress level, population differences emerged. We were particularly interested in genes that differed in their transcriptional response between the two stress levels, as traits that show a ‘dose-response’ may be more functionally relevant to osmotic challenge than traits that show a more general response. Finally, we sampled fish at multiple time points during an acclimation period to capture the temporal dynamics of this multi-phased physiological adjustment. Measurement of plasma ions and electron microscopy of gill morphology was used to characterize the nature and time course of common as well as population- and salinity-dependent adjustments at the organ and organ-system levels. Transcriptomics was used to characterize the dynamics of genome response to osmotic challenge, and to generate hypotheses about the important genes, pathways and higher-order networks that enable the extreme physiological flexibility that we observed.

Adult Fundulus heteroclitus (Linnaeus 1766) were collected from coastal New Hampshire (NH population) and coastal Georgia (GA population) and acclimated to full-strength seawater (32 p.p.t.) in a common re-circulating system for at least 6 months prior to experiments. Re-circulating systems were equipped with mechanical, biological and ultraviolet filtration. For osmotic challenge experiments, fish were manually transferred to 32 p.p.t. (control), 0.4 p.p.t. or 0.1 p.p.t. re-circulating systems and maintained at a density of one fish (4 g average) per gallon. Water chemistry for the 0.1 p.p.t. treatment was as follows: [Na⁺]=293±8.7 μmol l^–1, [Cl^–]=1250±32.8 μmol l^–1, [Ca²⁺]=337±11.7 μmol l^–1, [Mg²⁺]= 117±3.8 μmol l^–1 and [K⁺]=58±2.6 μmol l^–1 (N=8). The 0.4 and 32 p.p.t. treatments were approximately 1.1 and 91.4% of full-strength seawater, respectively. The concentrations of nitrogenous waste products, which were measured approximately once every 2 days (N=8), were always non-detectable. The photoperiod was 12 h:12 h light:dark, and temperature was kept at ambient (approximately 21°C).

Tissues were sampled pre-transfer (32 p.p.t.; 0 h), and at 6 h, 24 h, 72 h (3 days), 168 h (7 days) and 336 h (14 days) post-transfer to hypo-osmotic treatments for gene expression and physiology measures. Tissues were also collected from the 32 p.p.t. transfer time course for physiology measures. We chose these salinities because several studies have suggested that the gills of F. heteroclitus transition between a seawater and freshwater morphology at approximately 0.4 p.p.t. (Copeland, 1950; Philpott and Copeland, 1963). Exposure to 0.4 p.p.t. is the most extreme challenge where both populations were capable of maintaining homeostasis, whereas challenge down to 0.1 p.p.t. revealed compromised homeostatic maintenance and population variation in the pace of recovery of osmotic balance. Exposure to less severe hypo-osmotic challenge may reveal different homeostatic mechanisms, and may be the focus of future studies. At each sampling time point, plasma was collected by caudal puncture and was immediately isolated for ion measurements, one gill basket was excised and preserved for electron microscopy, and the other gill basket was preserved in RNAlater (Ambion, Austin, TX, USA) for transcriptomics. This experimental procedure is directly comparable to that of previous studies (Whitehead et al., 2011a; Whitehead et al., 2011b). All fish were treated in accordance with approved institutional animal care and use protocols.

Fish were killed by cephalic dislocation and sampled immediately for blood by cardiac puncture using ammonium-heparinized hematocrit capillary tubes. Blood was spun at approximately 2000 g for 2 min, and plasma was collected and frozen in liquid nitrogen prior to processing. Plasma was analyzed for osmolality using a freeze-point osmometer (Micro-osmometer 5004, Precision Systems, Natick, MA, USA), sodium by flame atomic absorption spectrophotometry (Varian AA240FS, Agilent Technologies, Santa Clara, CA, USA), and chloride using a modification of the mercuric-thiocyanate method (Zall et al., 1956). All plasma chemistry data are expressed as means ± s.e.m. For plasma chemistry, we used a one-way ANOVA with a Student–Newman–Keuls test, with an alpha level of 0.05 (SPSS version 17, IBM, Armonk, NY, USA). Fish gills were washed briefly in deionized water and then each gill arch was isolated. The first two left gill arches were immersed in fixative buffer (2% glutaraldehyde, 1% formaldehyde in 0.1 mol l^–1 cacodylate) for 4 h, rinsed in 0.1 mol l^–1 cacodylate buffer containing 0.02 mol l^–1 glycine, post-fixed in 2% osmium tetroxide for 1 h, and rinsed in water. Samples were dehydrated with an ethanol series, dehydrated with liquid CO₂ in a Denton critical point dryer (Cherry Hill, NJ, USA), mounted on aluminum scanning electron microscope (SEM) stubs, and coated with a gold:palladium solution (60:40 v:v) in an Edwards S150 sputter coater (Ashbord, Kent, UK). Random images of the afferent edge of gill filaments were captured using a JSM-6610 high vacuum mode SEM (Tokyo, Japan).

Total RNA was isolated from gill tissues following homogenization in Trizol reagent (Invitrogen, Carlsbad, CA, USA) and purified using spin columns. Antisense RNA (aRNA) was prepared using the MessageAmp II aRNA amplification kit (Ambion, Austin, TX, USA), and high-quality total RNA and aRNA was verified using microfluidic electrophoresis (Experion instrument and reagents, Bio-Rad Laboratories, Hercules, CA, USA). Purified aRNA was then coupled with Alexa Fluor dyes (Alexa Fluor 555 and 647, Molecular Probes, Eugene, OR, USA), coupled aRNA was cleaned up using spin columns, and pairs of samples were competitively hybridized to custom oligonucleotide microarrays (Agilent Technologies, Santa Clara, CA, USA) following the manufacturer’s protocols. Five biological replicates were hybridized per treatment, and each replicate sample was hybridized twice (including a dye swap) in a loop design where the treatment source of paired samples was balanced across treatments. Samples were enclosed in a SureHyb microarray hybridization chamber (Agilent Technologies) and incubated at 60°C for 18 h, following which slides were washed according to the manufacturer’s protocols and scanned using a Packard BioScience Scanarray Express (PerkinElmer, Waltham, MA, USA), and images were processed using ImaGene (BioDiscovery, Inc., El Segundo, CA, USA). The microarray platform was the same as that used in Whitehead et al. (Whitehead et al., 2011a; Whitehead et al., 2011b). Probes (N=6800) were designed from F. heteroclitus expressed sequence tag (EST) sequences and included all ESTs that could be annotated. The 6800 probes (representing 6800 UniGene target sequences) were printed in duplicate, including control spots, within the Agilent 15K-element 8-plex platform (Agilent eArray design ID 021434). Spots that were too bright (signal saturated) or too faint (signal less than two standard deviations above background intensity) were excluded from normalization and analysis.

Raw data were LOWESS-normalized, followed by mixed model normalization where ‘dye’ was treated as a fixed effect and ‘array’ was treated as a random effect, and normalized data were log₂-transformed prior to statistical analysis. Normalized data were analyzed using mixed model ANOVA in JMP Genomics (SAS Institute, Cary, NC, USA), where main effects included ‘population’ (NH or GA), ‘salinity’ (0.4 or 0.1 p.p.t. transfer) and ‘time’ (0, 6, 24, 72 and 168 h). An effect of salinity was detected by comparing expression post-transfer with that pre-transfer, rather than with time-matched no-transfer controls. The model was fully specified, including all three two-way interaction terms and the single three-way interaction term. ‘Dye’ was treated as a fixed effect, and ‘array’ and the five replicate individuals within each treatment were treated as random effects. Main effects and interactions were considered statistically significant at P<0.01, and false discovery rates were estimated using the q-value (Storey and Tibshirani, 2003) in R (R Foundation for Statistical Computing, Vienna, Austria). Network analysis was performed using Ingenuity Pathway Analysis (Ingenuity Systems Inc., Redwood City, CA, USA), and gene ontology enrichment analysis was performed using the Functional Annotation Clustering tool within the DAVID Bioinformatics Database (Huang et al., 2009).

Upon hypo-osmotic shock, blood plasma is expected to become diluted if fish are not capable of regulating net flux by either reducing diffusive ion loss or stimulating active ion uptake. Remarkably, F. heteroclitus did not lose osmotic balance when transferred from seawater (32 p.p.t.) to 0.4 p.p.t. (Fig. 1) at any time point post-transfer, suggesting an excellent propensity for ionoregulatory compensation. Only when severely challenged by exposure to 0.1 p.p.t. did fish temporarily lose osmotic homeostasis, and differences between populations emerged. Plasma osmolality and Na⁺ concentrations dropped immediately (Fig. 1A,B). Recovery to levels of 32 p.p.t. controls was within 3 days in the NH population but required up to 2 weeks in the GA population. The NH population was especially proficient in maintaining plasma Cl^– concentrations, as concentrations did not differ from controls at any time post-transfer even at 0.1 p.p.t. (Fig. 1C), whereas the GA population lost Cl^– balance immediately at 0.1 p.p.t., and again required up to 2 weeks for plasma Cl^– to return to control concentrations (Fig. 1C). This population variation is consistent with that found in a previous study, indicating that northern populations of F. heteroclitus tolerate hypo-osmotic shock better than southern populations (Scott et al., 2004), and that enhanced freshwater tolerance is associated with improved plasma Cl^– regulation. Furthermore, the northern clade of F. heteroclitus shares a similar genetic background with freshwater populations (Smith et al., 1998; Duvernell et al., 2008; Whitehead, 2009), and these freshwater populations are also more tolerant of hypo-osmotic shock, and refractory to impairment of plasma Cl^– concentrations, than populations with a southern genetic background (Whitehead et al., 2011a).

Although the mechanistic basis of regulation of Cl^– loss is not well understood, tight junction proteins, which regulate the conductance and the relative ionic permeability of the paracellular pathway, are likely influential in conferring a low paracellular loss of Cl^– in the face of osmotic stress (Wood and Grosell, 2008; Wood and Grosell, 2009). In addition to reducing ion efflux following hypo-osmotic exposure, restoration of ion homeostasis necessitates activation of active ion transport at the gill, a process most typically performed by mitochondrion-rich cells in the gill epithelium. Remodeling of the gill epithelium is a well-known core component of the osmotic acclimation response in euryhaline fishes. The gills in the NH population made a quick and full transition from a seawater to a freshwater phenotype only 1 day after hypo-osmotic exposure, whereas the GA population maintained a mixed seawater and freshwater gill mitochondrion-rich cell phenotype even after 14 days (Fig. 2). These results are consistent with those of Scott et al. (Scott et al., 2004) and Whitehead et al. (Whitehead et al., 2011a), who also demonstrated subtle population differences in the expression of gill cells on the epithelial surface. In our study, many genes associated with cell cycle regulation, cell growth and cell morphology were upregulated early and transiently, suggesting a putative role in epithelial plasticity during osmotic acclimation.

Acclimation to hypo-osmotic environments was associated with a substantial and complex genome expression response in gills of F. heteroclitus. One-third of transcripts included in the analysis (1177 transcripts) showed differential expression with time (P<0.01, q<0.009). A more nuanced analysis of hypo-osmotic transcriptional response explored acclimation-associated transcription that was dependent on population (NH vs GA) or salinity (0.1 vs 0.4 p.p.t.). A core set of 979 genes showed an acclimation time-course expression response that was common across both populations and both salinities (no significant interactions; Fig. 3A). A population-dependent acclimation time course response was detected for 131 transcripts (interaction between main effects ‘time’ and ‘population’, P<0.01, q<0.175; Fig. 3B). A salinity-dependent acclimation time course response was detected for 218 transcripts (interaction between main effects ‘time’ and ‘salinity’, P<0.01, q<0.143; Fig. 3C). The general pattern of temporal expression that we observed for genes showing a salinity-dependent time course response was a transient change in expression level at 0.4 p.p.t., but a sustained change in expression at 0.1 p.p.t. (Fig. 3C).

Aquaporin proteins are integral membrane proteins that facilitate the movement of water and other small solutes across biological membranes (Cerda and Finn, 2010). We found that aquaporin 3 (AQP3) was quickly upregulated upon hypo-osmotic challenge (Fig. 4A), with mRNA more strongly upregulated in the NH population than in the GA population immediately (6 h) following transfer, and upregulation tending to be more sustained at 0.1 p.p.t. but transient at 0.4 p.p.t. (P<0.01). NH fish retained osmotic homeostasis at 0.4 p.p.t. and AQP3 was only transiently upregulated, whereas GA fish at 0.1 p.p.t. took the longest to recover osmotic homeostasis and upregulation of AQP3 was sustained. Previous studies have shown transcriptional upregulation of AQP3 paralogs in the gills of freshwater-acclimated teleosts, though differences in cellular localization in branchial epithelia are evident among species. Some studies have described both basal and apical localization, suggesting a putative role in transepithelial water flux to minimize swelling or dehydration (Cutler and Cramb, 2002; Lignot et al., 2002). AQP3 expression may also facilitate other physiological functions, including ammonia or urea permeability, osmoreception or regulatory volume decrease. Challenge with hypo-osmotic conditions poses the immediate crisis of cell swelling, requiring the regulation of cell volume. Recently, AQP3 has been shown to facilitate preferential cellular uptake of water under hypotonic conditions, which facilitates an early increase in cell volume and a concomitant compensatory response aimed at restoring extracellular osmolality (Watanabe et al., 2005; Watanabe et al., 2009). Alternatively, AQP3 may be critical for the regulatory volume decrease (RVD) required after hypotonically induced cellular swelling, rather than in stimulating the initial cellular swelling as described in mammalian sperm (Chen et al., 2011). In general, RVD is achieved by activating transport systems that mediate the efflux of osmolytes out of the cell so that water will passively follow (Hochachka and Somero, 2002). Interestingly, we found that inositol monophosphatase (IMPA1) and versican (VCAN) were both downregulated in response to hypo-osmotic conditions (Fig. 4B). IMPA1 regulates synthesis of the intracellular osmolyte myo-inositol (Michell, 2008) and VCAN attracts cations to retain cellular water and mediate regulatory volume increase (Faggian et al., 2007). This downregulation showed a population-dependent trend, although it was not statistically significant (P=0.08 for IMPA1, P=0.11 for VCAN). Reduction of intracellular osmolytes, together with the expression of basolateral AQP3, may contribute to RVD following the initial cellular swelling that is likely associated with significant reduction in plasma osmolytes.

Genes involved in the synthesis of polyamines (putrescine, spermine and spermidine) were also regulated in response to hypo-osmotic challenge. Ornithine decarboxylase 1 (ODC1) is the rate-limiting enzyme in polyamine synthesis, and showed a pattern of population- and salinity-dependent time-course expression highly correlated with AQP3 expression (Fig. 4C). ODC1 also shows species-specific regulation among native and invasive mussels upon hypo-osmotic challenge (Lockwood and Somero, 2011). Expression of arginase 2 (ARG2), the enzyme that synthesizes ornithine, showed a similar pattern of salinity-dependent expression and peak expression was delayed relative to ODC1 (Fig. 4C). Regulation of two additional genes was consistent with a role of increased polyamines in the hypo-osmotic compensatory response. Expression of a spermine-binding protein (SBP) quickly dropped to nearly non-detectable levels, and transglutaminase 1 (TGM1), which facilitates polyamine-to-protein conjugation, quickly increased in expression (Fig. 4D). Expression of these genes is consistent with increased cellular free polyamine levels. Polyamines serve many cellular functions in a wide range of organisms (Rhee et al., 2007). One of these functions is to stabilize protein–protein and protein–nucleic acid interactions against changes in intracellular osmolyte levels (Capp et al., 1996). As proposed previously (Whitehead et al., 2011a), increased cellular polyamines may be important for stabilizing genome function in the face of rapid changes in osmolyte concentrations during osmotic challenge in fish, and may serve a similar role in bivalves (Lockwood and Somero, 2011).

Previous studies have indicated that regulation of histone gene expression is associated with osmotic acclimation in killifish (Whitehead et al., 2011b; Whitehead et al., 2011a), and we detected several histone transcripts with transient upregulation in response to hypo-osmotic challenge (Fig. 4E,F), two of which show salinity-dependent upregulation (Fig. 4E). Hypo-osmotic conditions can cause swelling of the nucleus and alter the structure and arrangement of chromatin, thereby potentially disrupting genome function (Finan and Guilak, 2010). Histone proteins contribute to maintaining the three-dimensional structure of chromatin, and their transient upregulation during the early crisis of cell swelling may be associated with conformational maintenance of chromatin.

The group of genes showing a salinity-dependent time-course response (Fig. 3C) was enriched with genes associated with the gene ontology term ‘transmembrane transport’ (P=0.02) and included ion-transporting genes. Two subunits of the sodium/potassium-transporting ATPase (ATP1A1 and ATP1A3) showed salinity-dependent upregulation at 1 day post-transfer or later, where they were more highly upregulated in both populations at 0.1 p.p.t. than at 0.4 p.p.t. (Fig. 4G). These results are consistent with an increased requirement for, and stimulation of, active ion uptake in the gill upon hypo-osmotic transfer. In addition, the V-type proton ATPase (ATP6V1E1) was moderately upregulated (Fig. 4G). This protein has been putatively localized to the basolateral epithelium of gill cells, where it may aid in transepithelial Na⁺ uptake in freshwater or perhaps other diverse functions associated with freshwater acclimation, including acid–base regulation (Patrick et al., 1997; Patrick and Wood, 1999; Scott et al., 2005).

In contrast to the sodium/potassium and proton ATPases, which are upregulated later in the time course, claudin 4 (CLDN4) was rapidly upregulated (Fig. 4H). CLDN4 followed a salinity-dependent response (P=1.0E–4) where upregulation was sustained at 0.1 p.p.t. but was transient at 0.4 p.p.t. Fundulus heteroclitus is unusual among euryhaline fish insofar as they do not actively accumulate chloride from freshwater. Rather, they rely on limiting chloride efflux to compensate for hypo-osmotic environments (Patrick et al., 1997; Patrick and Wood, 1999; Scott et al., 2004), or acquiring the ion from dietary sources. Claudins are an expansive family of cellular tight junction proteins integral to the regulation of paracellular permeability and ion selectivity (Van Itallie and Anderson, 2006). Recent studies have described the differential expression of several claudin genes in the ion-transporting epithelia of teleost fish subjected to osmotic challenges (Bagherie-Lachidan et al., 2008; Bui et al., 2010; Clelland et al., 2010). Although poorly described, claudins may function to limit chloride efflux through the paracellular tight junctions of Fundulus gills during hypo-osmotic transfer, providing a potential mechanism for altering the relative ionic permeability of gills during physiological transitions. Indeed, the group of genes showing a consistent time-course response (Fig. 3A) is enriched for genes associated with the gene ontology term ‘cell junction’, and includes occludin (OCLN) and two additional claudin proteins (CLDN4 and CLDN7) that were quickly upregulated upon hypo-osmotic challenge (Fig. 4I).

One of the most dramatically upregulated genes in our analysis was otopetrin 1 (OTOP1; Fig. 4J). OTOP1 showed a salinity-dependent transcriptional response (P=6.0E–12), where expression peaked by 24 h at 0.4 p.p.t. but at 0.1 p.p.t. continued to increase out to 7 days. OTOP1 is a protein that regulates intracellular calcium levels (Hughes et al., 2007), and expression knockdown interferes with otolith development in zebrafish (Hughes et al., 2004). Hypo-osmotic challenge is accompanied by an increase in intracellular calcium during cell swelling in various animal models, implying a role in regulatory volume decrease (Pasantes-Morales et al., 2006), and calcium is important for sensing and signal transduction during osmotic stress in diverse taxa (Fiol and Kültz, 2007). Though the dramatic salinity-dependent regulation of OTOP1 is intriguing, especially given the important and diverse roles of calcium during osmotic stress, the cellular localization of the protein in fish gills remains unknown, as does its functional relevance for osmotic acclimation.

Focal adhesion kinase (PTK2) downregulation indicated a salinity-dependent trend (P=0.08; Fig. 4K) upon hypo-osmotic challenge. The physiological relevance of this is unclear, but is intriguing because FAK phosphorylation is important for regulation of chloride channels during osmotic compensation in killifish MR cells (Marshall et al., 2009), and cystic fibrosis transmembrane conductance regulator chloride channels are deactivated and downregulated in freshwater (Singer et al., 1998).

The physiological transition from a seawater physiology to a freshwater physiology is energetically expensive, estimated to consume approximately 10% of the total energy budget in F. heteroclitus (Kidder et al., 2006a; Kidder et al., 2006b). We observed salinity-dependent upregulation of three probes for creatine kinase (two for CKM and one for CKB; Fig. 4L), as well as of glycogen phosphorylase (PYGM) and glucose-6-phosphate isomerase (GPI), which are involved in carbohydrate metabolism (Fig. 4M). For each of these genes, upregulation was more dramatic or more sustained at 0.1 p.p.t. than at 0.4 p.p.t., indicating that the more severe osmotic stress imposed a greater energetic cost to the animals. Mitochondrial electron transport chain genes were also upregulated, implying an increase in oxidative phosphorylation (see below). Gill PYGM was transiently upregulated in tilapia gills during seawater transfer (Chang et al., 2007), and this response correlated with enhanced gill ATPase activities; genes associated with energy production were also regulated in response to hypo-osmotic challenge in a euryhaline goby (Evans and Somero, 2008). These results attest to the high energetic demands required for the gills to transition between different osmotically induced phenotypes.

The core acclimation-associated genes (Fig. 3A) and the set of genes showing a salinity-dependent time-course response (Fig. 3C) were both enriched for genes associated with the gene ontology term ‘mitochondrion’, including various mitochondrial transport proteins (e.g. TIM8B, TOM40, SLC25A20 and SLC25A13; supplementary material Table S1). Intriguing but difficult to interpret are three observations about regulation of the mitochondrial genome and components of the mitochondrial electron transport chain (mtETC). First, we observed significant downregulation (P=3.0E–11) of a transcript with perfect sequence match to the F. heteroclitus mitochondrial D-loop (mt-Dloop; Fig. 5, line graph), showing a trend towards salinity dependence (P=0.08) and population dependence (P=0.05). This mt-Dloop transcript was also downregulated in response to hypoxia, confirmed by quantitative PCR (Flight et al., 2011). Second, we observed a coordinated gradual upregulation of most of the subunits encoding mtETC proteins, consistent with the high demand for ATP during osmotic acclimation (Fig. 5, heat map). However, this trend only applied for mtETC components that are encoded by the nuclear genome; the subset of mtETC transcripts encoded by the mitochondrial genome showed transient downregulation. Mitochondrial genome replication and mitochondrial transcription are both initiated by transcriptional machinery primed in the D-loop (Bonawitz et al., 2006), and D-loop transcripts are adenylated (Sbisà et al., 1990; Vijayasarathy et al., 1995; Camasamudram et al., 2003), so a decrease in D-loop transcript abundance could reflect decreased mtDNA synthesis and/or decreased transcription of the mitochondrial genome. However, this appears inconsistent with our observed increase in expression of components of the mtETC. That is, if mitochondrion abundance decreases, one would expect a parallel decrease in mtETC components. Also curious is the apparent conflict between coordinated upregulation of nuclear-encoded mtETC transcripts but transient downregulation of mitochondrion-encoded mtETC transcripts, as nuclear and mitochondrial regulation of mtETC gene transcription is tightly coordinated in mammals (Scarpulla, 2006) (though this regulation has not been as well studied in fish). Perhaps mitochondrion-encoded transcripts were in fact upregulated in coordination with nuclear-encoded transcripts, but the density of our temporal sampling was insufficient for capturing this coordinated regulation because of the shorter half-life of transcripts that are expressed in the mitochondrion (Murphy and Attardi, 1973; Gelfand and Attardi, 1981). Regulation of mitochondrion turnover and electron transport during osmotic acclimation merits more focused study because mitochondrion densities differ between seawater and freshwater chloride cells (Laurent et al., 2006), and regulation of ion flux is tightly coupled with electron transport-mediated synthesis of ATP (Bradley and Satir, 1981).

Genes that are differentially regulated during acclimation were enriched for the gene ontology terms ‘transcription regulator activity’ (P=0.009) and ‘transcription factor activity’ (P=0.02), consistent with the transcriptional response to hypo-osmotic challenge in a euryhaline goby (Evans and Somero, 2008). OSTF2 is a transcription factor that is quickly upregulated in response to hyper-osmotic challenge in tilapia and shares greatest sequence identity with the TSC-22 family of transcription factors (Fiol and Kültz, 2005). Transcripts for the killifish homolog to OSTF2 and killifish TSC22D3 showed coordinated early and transient upregulation upon hypo-osmotic challenge (Fig. 4N), correlated with expression of transcription factor JUNB. Though these transcription factors are quickly upregulated, their precise roles in mediating hypo-osmotic acclimation are unclear. TSC22D3 is induced by glucocorticoids and glucocorticoid receptor is involved in mediating osmotic acclimation (Shaw et al., 2007), and TSC22D3 plays a role in cellular protection from apoptosis (Asselin-Labat et al., 2004). Because OSTF1 contains putative binding sites for cyclin-dependent kinase 1 (CDK1) and MAP kinases, it may be involved in regulation of the cell cycle and apoptosis. JUNB can be activated by mechanical changes in membrane shape during hypo-osmotic exposure (Rauch et al., 2002), and expression is induced by hypo-osmotic exposure in keratin-defective epithelial cells along with a group of transcription factors that contribute to cell cycle regulation (Liovic et al., 2008). Perhaps regulation of these transcription factors is protective, where temporary growth arrest permits repair and prevents proliferation of damaged cells (Liovic et al., 2008). Iodothyronine deiodinase is not a transcription factor, but is the major activator of thyroid hormone, which is a well-known activator of transcription. Though the role of thyroid hormone in osmoregulation is controversial (McCormick, 2001), we found significant population-dependent upregulation of iodothyronine deiodinase 1 (DIO1; Fig. 4O), suggesting a role in salinity acclimation that may contribute to explaining observed population variation in hypo-osmotic tolerance (Fig. 1), consistent with adaptive divergence in expression reported in Whitehead et al. (Whitehead et al., 2011a). DIO proteins are important for the activation and deactivation of thyroid hormones (Orozco and Valverde-R, 2005). Members of the DIO family have been cloned from fish gills (Sanders et al., 1999), contain osmotic response elements in promoter regions in killifish (Lopez-Bojorquez et al., 2007), and deiodination activity is detectable in killifish gills (Orozco et al., 2000). Intriguingly, thyroid hormone signaling pathways have adaptively diverged between marine and freshwater populations of stickleback (Kitano et al., 2010).

The transcription factors HES1 and SPI1 were downregulated in response to hypo-osmotic challenge (Fig. 4P). HES1 is a transcriptional repressor and an effector of the notch-signaling pathway. Notch signaling is important for regulating proliferation and differentiation decisions and cell fate, where signaling can either inhibit or induce differentiation. HES1 regulation may reflect the involvement of notch signaling in the gill tissue remodeling that we observe. SPI1 is a transcription factor that regulates haematopoiesis, but also upregulates expression of cyclin D2 (Vicari et al., 2006), which also is important for cell-cycle regulation and cell proliferation.

Many genes involved in apoptosis were regulated in response to hypo-osmotic challenge, but gene ontology terms associated with ‘apoptosis’ were not enriched in our set of time-course responsive genes. The role of apoptosis during osmotic acclimation in fish gills is not clear. Change in intracellular potassium concentration, which is a key component of regulatory volume decrease following cell swelling (Chara et al., 2011), can trigger apoptosis (Bortner et al., 1997). However, upon hypo-osmotic challenge, killifish may suppress apoptosis to allow regulatory volume decrease to proceed followed by cellular remodeling, or may activate apoptosis if regulatory volume decrease is ineffective and cells start to accumulate damage. Laurent et al. (Laurent et al., 2006) suggested that during freshwater acclimation, seawater-type mitochondrion-rich cells may be eliminated by apoptosis, necrosis or exfoliation, or that they may be simply covered over by neighboring pavement cells. Many biochemical pathways can mediate apoptosis and signaling can be complicated and cell-type dependent (Strasser et al., 2000), and in the present study we found both upregulation and downregulation of several key pro-apoptotic genes, making broad interpretations difficult. However, considering the importance of tissue remodeling necessary for acclimation, and that osmotic stress-responsive signaling pathways are associated with apoptosis (Kültz and Avila, 2001), including those implicated in our network analysis (see below), we propose that regulation of apoptosis may be an important component of hypo-osmotic acclimation in killifish, a proposition that merits more focused study.

We observed significant regulation of several key mediators of apoptosis. Cytochrome c (CYC) transcription was rapidly upregulated by 6 h (P=1E–27; Fig. 4Q), was population (P=0.0005) and salinity dependent (P=0.0001), and time-course expression across populations and salinities was highly correlated with other genes associated with cell volume regulatory phenomena (e.g. AQP3 and ODC1; Fig. 4A,C). Release of mature CYC protein from the mitochondrion is a key trigger for apoptosis through enhanced activation of caspase, and upregulation of CYC transcription precedes induction of apoptosis (Chandra et al., 2002). In killifish, caspase-8 (CASP8) expression showed population-dependent upregulation (Fig. 4R) upon hypo-osmotic challenge. CASP8 is a key regulator of the apoptotic cascade, is centrally located within the network of apoptotic machinery components (Di Pietro et al., 2009), and increased expression is associated with enhanced sensitivity to apoptosis (Jiang et al., 2008). Cyclin-dependent kinase 1 (CDK1) upregulation was correlated with CASP8 expression (Fig. 4R). Activity of CDK1 is regulated by transcription and phosphorylation, and CDK1 activation is linked to both cell-cycle regulation and regulation of apoptosis (Castedo et al., 2002). Transcripts for the protein FADD and the apoptosis regulatory protein Siva (SIVA1) were both downregulated during the acclimation time course (Fig. 4S). FADD is pro-apoptotic through CASP8 activation, and SIVA1 is also pro-apoptotic through the caspase-dependent pathway (Py et al., 2004).

We exploited gene and protein interaction network analyses as an exploratory tool for revealing regulatory pathways that may activate and coordinate the complex genomic responses that facilitate physiological acclimation. Network analysis seeks connections between genes of interest often from experimental evidence, such as protein–protein interactions or protein–nucleic acid interactions, to provide greater scope for biological systems-level functional inference. We focused on connections among three sets of genes: (1) genes that showed a common (‘core’) transcriptional response during acclimation (‘time’ main effect with no population or salinity interaction), (2) genes that showed a salinity-dependent time-course response (interaction between ‘time’ and ‘salinity’ main effects), and (3) genes that showed a population-dependent time-course response (interaction between ‘time’ and ‘population’ main effects). These are the three sets of genes represented by different heat maps in Fig. 3, but for network analysis we included only the subset of genes with functional annotation. We identified hubs that united more than 15 of these genes into an interaction network with one degree of separation. Onto each network we mapped whether each gene showed a ‘core’ transcriptional time-course response (pink symbols; Fig. 6), a salinity-dependent time-course response (red symbols) or a population-dependent time-course (green symbols). We expected this analysis to reveal well-known osmotic stress response pathways, but in particular we were interested in identifying the subset of hubs that were significantly enriched for genes with salinity- and population-dependent expression; regulation from these hubs may be particularly important for imparting the extreme physiological plasticity observed in F. heteroclitus.

One major cluster of highly connected hubs included well-known components of the cytokine and kinase cellular signaling pathways (Fig. 6). This, the largest cluster that emerged from our analysis, united 12% of the time-course responsive genes, but was not significantly enriched for genes with salinity- or population-dependent expression (P=0.12). This network links 10 hubs, including cytokines TGF-β1 and TNFα, the receptor tyrosine kinase ERBB2, phosphatidylinositol 3-kinase (PI3K), the protein kinases PKC, JNK, ERK1/2 and P38 MAPK, and the transcription factors NFκB and AP-1, and each of these hubs is highly connected with one another (blue dashed lines; Fig. 6). Many of these hubs are well-known components of osmotic stress response pathways, confirming the validity of our networking approach. For example, cell swelling upon hypo-osmotic challenge activates stress-activated protein kinases (Niisato et al., 1999), and MAPK-mediated signaling regulates the cell cycle as an important component of the osmotic stress response (Kültz and Burg, 1998). Indeed, osmotic challenge quickly regulates ERK1, JNK and P38 MAPK protein abundance in killifish gill and opercular epithelium (Kültz and Avila, 2001; Marshall et al., 2005), hypo-osmotic challenge has been shown to induce MAPK signaling in an estuarine goby (Evans and Somero, 2008) and MAPK signaling has been shown to affect regulatory volume decrease in various other fish models (reviewed in Chara et al., 2011).

In addition to MAPK signaling, other inter-connected signaling pathways including TGF-β and TNFα/NFκB are also involved in cell proliferation, cell differentiation and cell death (Massague, 1998; Liu, 2005). Intriguingly, only one base pair distinguishes the osmotic response element (ORE) and the NFκB binding sites, and Iwata et al. (Iwata et al., 1999) found that TNFα induces aldose reductase by stimulating NFκB binding to the ORE of aldose reductase. Perhaps TNFα/NFκB activation can stimulate expression of other genes with ORE binding sites. TNFα affects tight-junction permeability (Johnson, 2006), which in turn regulates ion efflux from gill epithelia in hypo-osmotic conditions. TNFα is involved in cellular localization of occludin (Sappington et al., 2003), decreases expression of rat occludin in sertoli cells (Siu et al., 2003) and downregulates occludin in astrocytes (Wachtel et al., 2001). JUNB expression was quickly upregulated upon hypo-osmotic challenge (Fig. 4N) and increased JUNB transcription is a sensitive and early genomic response in cells stimulated by TGF-β1 (Pertovaara et al., 1989). JUNB induction accompanies TGF-β1-mediated induction of ODC in H-ras transformed cell lines (Hurta et al., 1993), and this pattern of co-expression between ODC and JUNB was observed in our data as well (Fig. 4C,N).

This network (Fig. 6) is also enriched for genes involved in regulation of apoptosis (P=0.006). MAP kinases are well-known mediators of cell death through activation of apoptotic cascades. TNFα can also induce apoptosis by triggering activation of CASP8, and osmotic stress can induce TNFα-mediated apoptosis in human cell lines (Franco et al., 2002). In contrast, NFκB can play anti-apoptotic roles through activation of apoptosis inhibitors (Liu, 2005). Again, the fact that molecular signaling governing apoptosis is enormously complicated makes it difficult to conclude activation or repression from transcript profiling alone. However, given the necessary damage control and tissue remodeling necessary for hypo-osmotic acclimation, observed changes in transcript abundance for key mediators of apoptosis (Fig. 4Q–S), and implication of signaling pathways known to be key mediators of cell death (Fig. 6), regulation of apoptosis may be an important component of hypo-osmotic acclimation in killifish and merits more focused study.

Our analysis revealed two sets of hubs that were enriched for genes showing salinity- and/or population-dependent patterns of expression during acclimation to hypo-osmotic challenge. Of all the hubs identified in our analysis, HNF-4α united the most genes (11% of the genes included in this analysis), and this network was the most significantly enriched for genes showing salinity- (P=0.0004) and population-dependent expression (P=0.0009; Fig. 7A). HNF-4α is a transcription factor crucial for liver development, though it is expressed in some mammalian osmoregulatory tissues, such as kidney and intestine, for which its function remains to be determined (Grigo et al., 2008). The gene is also expressed in fish gill, where zinc deficiency induces transcription of HNF-4α and of genes enriched with HNF-4α binding elements (Zheng et al., 2010). Though a role in fish osmoregulation is novel, a homolog of HNF-4α appears crucial for coordinating osmoregulatory functions of excretory cells in nematode worms. In Caenorhabditis elegans, HNF-4α regulates the expression of many subunits of the vacuolar ATPase (vATPase) (Hahn-Windgassen and Van Gilst, 2009). The vATPases are proton pumps that are found in all eukaryotes (Finbow and Harrison, 1997), and in fish gills they are important in energizing transepithelial Na⁺ uptake in freshwater-acclimated fish and in acid–base balance. In C. elegans, RNA interference of HNF-4α compromises the animal’s ability to compensate for acute osmotic stress (Hahn-Windgassen and Van Gilst, 2009). In F. heteroclitus, we detected salinity-dependent regulation of ATP6V1E1 (Fig. 4G) and population-dependent regulation of ATP6V1D (supplementary material Table S1).

A second set of mutually connected hubs, including the hormone β-estradiol, transcription factors MYC and TP53, and cytokine interferon-gamma (IFNG), was significantly enriched for genes with salinity- (P=0.0001) and population-dependent (P=0.013) expression (Fig. 7B). Estradiol injection typically has inhibitory effects of osmoregulatory capacity in fish (McCormick, 1995), including killifish (Mancera et al., 2004). Consistent with this, transcript for estradiol 17β-dehydrogenase, which catalyzes the synthesis of estradiol, was significantly decreased upon hypo-osmotic challenge in a euryhaline goby (Evans and Somero, 2008). IFNG is a cytokine primarily involved in regulation of immune system genes in macrophages, but it also exhibits anti-proliferative and pro-apoptotic functions (Schroder et al., 2004). Similarly, transcription factor TP53 affects cell fate and proliferation by inducing growth arrest or apoptosis. MYC also regulates cell growth, cell proliferation and apoptosis, and has been shown to be downregulated upon hypo-osmotic challenge in two species of mussel (Lockwood and Somero, 2011).

It is important to note that all of the sets of hubs highlighted in (Figs 6, 7) are interconnected (Fig. 8A). Interestingly, the least connected hub is HNF-4α, which is also the hub most significantly enriched for genes with salinity- and population-dependent expression, and is the only hub to unite more of the genes included in our analysis than expected by chance (Fig. 8B). Perhaps this relative lack of interconnectedness with other hubs that are responsive to osmotic challenge is associated with salinity-dependent flexibility in response of HNF-4α signaling, and with the likelihood of divergence of that signaling response among populations. Indeed, less connected components of molecular networks appear subject to fewer selective constraints than highly connected components, facilitating ‘evolvability’ (Yamada and Bork, 2009).

In addition to non-transcriptional responses (such as protein phosphorylation), it is clear that regulation of genome expression at the transcriptional level is a key mediator of acclimation to changing osmotic environments. In the experiments described here, two populations of killifish differ in their abilities to compensate for osmotic challenge (Fig. 1), in the capacity of their gill epithelia to transition into a freshwater state, and in their differential expression of genes. Most of the genes highlighted above are expressed differently between the NH and GA populations, and this divergence could have been governed by neutral or adaptive evolutionary forces (Whitehead and Crawford, 2006a; Whitehead and Crawford, 2006b). That alternate osmotic physiologies have evolved repeatedly and quickly within the genus (Whitehead, 2010) implies that the genomic mechanisms that govern osmoregulation can also evolve quickly. One might predict that for genes whose expression regulation is particularly important for osmotic acclimation, adaptive forces should govern their divergence among populations that have adapted to different osmotic environments. It is unclear whether increased freshwater tolerance is currently or historically adaptive in the coastal NH population. However, expression patterns for many of the genes highlighted above (AQP3, ATP1A1, CYC, DIO1, HES1, IMPA1, JUNB, ODC1, OSTF1, SBP, TGM1, TSC22D3 and VCAN) appear to have diverged by adaptive evolutionary mechanisms among populations of killifish distributed along a natural salinity gradient (Whitehead et al., 2011a). That is, many of these genes show patterns of adaptive divergence in gene expression in populations that occupy exclusively freshwater habitats and that have derived adaptive physiological acclimation to hypo-osmotic challenge. This emphasizes the physiological importance of expression regulation during acclimation for these genes.

ACKNOWLEDGEMENTS

Thanks to Benjamin Dubansky, Michael Truax, Lee McChesney, John Waldron and Alex Whelan for assistance with animal care.