Reciprocal osmotic challenges reveal mechanisms of divergence in phenotypic plasticity in the killifish Fundulus heteroclitus

For aquatic organisms, salinity is one of the main factors limiting species distributions (Kolar and Lodge, 2002; Lee and Bell, 1999; Whitehead, 2010). Within fishes, large clades are generally restricted to either freshwater or marine habitats (Nelson, 2006; Vega and Wiens, 2012); however, some closely related species inhabit different osmotic niches and have divergent osmotic physiologies (Dunson and Travis, 1991; Plaut, 1998; Whitehead, 2010). It is thought that species from habitats characterized by variable salinities, such as estuaries, are those most able to invade novel osmotic niches (Lee and Bell, 1999). For fishes, this suggests that euryhalinity may enable the exploitation of diverse osmotic niches (Bamber and Henderson, 1988; Schultz and McCormick, 2013). We seek to identify the physiological and genetic mechanisms underlying euryhalinity, a plastic phenotype that enables exploitation of novel osmotic niches and serves to facilitate radiations across osmotic boundaries.

Species within the genus Fundulus occupy the continuum of osmotic niches and have undergone multiple invasions from marine to freshwater systems (Whitehead, 2010). Fundulus species are derived from an ancestral marine or estuarine species and harbor differences in physiological plasticity and tolerances to environmental salinity (Griffith, 1974; Whitehead, 2010). The estuarine fish Fundulus heteroclitus has further evolved extreme osmotic plasticity, allowing it to tolerate salinities ranging from fresh water to three times the salinity of seawater (Griffith, 1974), consistent with its broad osmotic niche spanning marine to freshwater environments (Hildebrand and Schroeder, 1928). A genetic discontinuity exists between geographically proximate populations along parallel salinity clines that transition between brackish and fresh waters of the eastern US despite no apparent barriers to gene flow. Freshwater populations show an increased capacity to acclimate to hypo-osmotic stress when compared with their brackish counterparts, consistent with adaptive divergence along this salinity gradient (Whitehead et al., 2011).

Phenotypic plasticity can evolve to permit population persistence in rapidly changing environments and can enable rapid convergence toward a new phenotypic optimum (Lande, 2009). Indeed, marine copepods that have recently invaded freshwater habitats have evolved expanded osmoregulatory plasticity (Lee et al., 2011). However, given enough time in a new stable environment, where plasticity is no longer favored, plasticity may be assimilated (Lande, 2009) because it can be costly to maintain (DeWitt et al., 1998). Osmoregulatory plasticity probably consists of the ability to effectively sense a change in the osmotic environment and transduce signals to alter physiological processes that lead to cell volume regulation and remodeling of ion-transporting epithelia. Since this osmotic compensatory response is energetically expensive (Kidder et al., 2006a,,b), deterioration of the sensing, signaling and effector mechanisms that support this response may be expected under relaxed selection (Lahti et al., 2009). While F. heteroclitus populations derived from freshwater habitats, but reared in the laboratory in sea water, have a higher capacity for hypo-osmoregulation relative to populations from brackish and marine sites (Whitehead et al., 2011), their tolerance to hyper-saline water is unknown. In killifish, migration into freshwater habitats may be associated with relaxation of stabilizing selection that maintains the phenotypic plasticity necessary to tolerate salinity fluctuations typical of estuarine habitats.

We employed a comparative experimental design to investigate physiological and transcriptomic responses over a time-course of acclimation for populations of F. heteroclitus native to fresh water (FW-native) and brackish water (BW-native) niches during both high- and low-salinity challenges. By combining the global approach of transcriptomics with physiology, we infer mechanisms underlying osmotic phenotypes through the identification of evolutionarily conserved responses. Because changes in selective pressures in different osmotic niches may lead to divergence in acclimation responses between populations, this study provides insight into the mechanisms responsible for adaptive divergence between osmotic niches (Whitehead, 2012).

As a result of the differential selective pressures in each environment, we hypothesized that each population would exhibit its most efficient acclimation in its native environment and that transcriptional responses could be used to illuminate the mechanisms underlying osmoregulation capabilities. Exposure to low salinity (hypo-osmotic) requires aquatic organisms to regulate osmolytes against an osmotic gradient, whereas in high salinity (hyper-osmotic), osmolytes must be actively excreted to counteract the passive entry of salts (Evans and Claiborne, 2006; Hwang et al., 2011). Therefore, one may predict that maintenance of homeostasis in the different osmotic environments relies on unique molecular pathways. However, one may also predict mechanisms that share a common response to osmotic stress regardless of the direction (hypo- or hyper-osmotic), including cell signaling, macromolecule damage control and other stress mechanisms (Fiol and Kültz, 2007). To understand the pathways contributing to euryhalinity in fishes, it is necessary to focus on both the common and divergent regulatory responses to hypo- and hyper-osmotic stress. Furthermore, transcriptomic responses that differ between populations may be mechanistically linked to divergence in physiological capabilities. We investigated three patterns of transcriptional response to salinity stress: (1) a response that is common in both salinity challenges and conserved between populations; (2) a response that differs between each challenge but is conserved between populations; and (3) a population-specific response that represents evolved divergence between the populations. Patterns 1 and 2 may facilitate identification of mechanisms that broadly contribute to euryhaline capabilities, whereas pattern 3 may identify mechanisms and pathways that underlie alternative osmotic phenotypes that have evolved in different osmotic niches.

Hypo and hyper-osmotic environments provide opposing challenges for aquatic organisms and require that the rates and directionality of active ion flux in transporting epithelia such as the gill be altered to adjust to each osmotic environment (Evans et al., 2005). Generally, failure to acclimate during a hypo-osmotic challenge manifests as dilution of body fluids such as blood, whereas failure to acclimate to hyper-osmotic challenge has a concentrating effect, each followed by varying degrees of recovery of osmotic homeostasis. In the hypo-osmotic challenge, the three-way ANOVA revealed a time×salinity interaction where both populations showed a parallel loss of plasma Na⁺ concentration followed by a rapid return to control values (P<0.001, Fig. 1A). There were differences between populations in the capacity to regulate extracellular fluid (ECF) Cl⁻ during hypo-osmotic challenge, where the BW-native population lost ECF Cl⁻ balance but the FW-native population did not (population×salinity interaction, P=0.04; Fig. 1B). For the BW-native population, post hoc analyses showed a salinity×time interaction (P=0.03), where plasma Cl⁻ concentration was significantly reduced by 6 h and did not return to pre-transfer control values (but was not significantly different from time-matched controls after 6 h post transfer), indicating that the BW-native fish struggled to recover Cl⁻ homeostasis even after a 2 week period of acclimation. These results indicate that the FW-native population maintains or recovers internal osmotic homeostasis better than the BW-native fish when exposed to a low salinity challenge, suggesting adaptation to its osmotic niche. Other killifish studies have contrasted Atlantic coastal populations from northern and southern regions, where the FW-native population shows greatest genetic similarity to fish from the northern clade in both nuclear (Duvernell et al., 2008) and mitochondrial (Whitehead, 2009) markers. Northern fish display improved ability to maintain internal ion homeostasis in response to low salinity challenge compared with southern fish (Whitehead et al., 2012b; Scott et al., 2004).

After hyper-osmotic exposure, no significant population-dependent interactions were detected for Cl⁻ in the three-way ANOVA (Fig. 1D). However, the subtle trend towards increased Cl⁻ for the BW-native population may warrant further investigation. In contrast, the regulation of plasma Na⁺ concentration differed significantly between populations, where the FW-native population lost ECF homeostasis but the BW-native population did not (time×population interaction, P=0.002; Fig. 1C). Post hoc tests showed that for the FW-native population alone, a salinity×time effect was detected (P=0.004), where plasma Na⁺ increased significantly at 6 h and 1 day but recovered by 3 days.

These physiological responses reveal differences in regulatory abilities between populations where each fails to maintain homeostasis of a different ion in different osmotic environments (Fig. 1C,D). The population divergence revealed by freshwater challenges is consistent with other studies, but the population variation revealed by brackish water challenges is not (Scott et al., 2008,, 2004; Scott and Schulte, 2005). Scott et al. (2004) reported F. heteroclitus population divergence in Cl⁻ regulation following freshwater challenges. However, previous high-salinity challenges showed no loss of osmolyte balance and no difference between northern and southern populations (Scott et al., 2008; Scott and Schulte, 2005). These studies acclimated animals to brackish water [3.5 parts per thousand (p.p.t.) or 10 p.p.t.] prior to a full strength seawater (∼35 p.p.t.) challenge, while animals in the experiments reported here were acclimated to fresh water (0.2 p.p.t.) before a half strength seawater (15 p.p.t.) challenge. Killifish acclimated above ∼1 p.p.t. maintain a seawater gill morphology, whereas those acclimated below ∼1 p.p.t. maintain a freshwater gill morphology (Whitehead et al., 2012b). Thus, we find that osmotic challenges spanning the salinity threshold across which gills are remodeled (∼1 p.p.t.) are more likely to reveal physiological differences between populations than challenges that do not.

This remodeling threshold also corresponds to a genetic break between the FW-native and BW-native populations in parallel salinity clines, consistent with salinity representing a barrier to gene flow in nature (Whitehead et al., 2011). It is unclear whether population differences in osmotic acclimation through gill remodeling are linked to fitness differences in alternative osmotic habitats. While the high energetic costs associated with osmoregulation and remodeling may have long-term effects (Kidder et al., 2006a,,b), it is possible that other mechanisms are contributing to the maintenance of the genetic break, such as compromised fertilization or embryonic development, which are currently under investigation.

Our results are consistent with extended physiological plasticity at hypo-osmotic extremes in populations that have invaded freshwater habitats, accompanied by loss of salt tolerance capabilities towards the hyper-osmotic end of the salinity continuum when compared with the BW-native population. Although gain of hypo-osmotic plasticity is probably adaptive as estuaries transition to fresh water, loss of hyper-osmotic plasticity may be a result of relaxation of selection for osmotic plasticity in an environment that no longer experiences high salinity (Lahti et al., 2009). Conversely, divergence may be adaptive in fresh water because of the costs associated with maintaining a plastic phenotype (DeWitt et al., 1998). At the macro-evolutionary time-scale, repeated speciation from saltwater to freshwater environments within Fundulus is accompanied by loss of hyper-osmotic plasticity (Whitehead, 2010).

That BW-native fish showed a reduced ability to regulate ECF Cl⁻ in fresh water, while FW-native fish showed a reduced ability to regulate ECF Na⁺ in brackish water, suggests that cell–cell junctions are involved in the divergence of physiological abilities. In F. heteroclitus, gills consist of a heterogeneous mix of cells, including impermeable pavement cells, mitochondrion-rich cells (MRCs) and accessory cells (ACs) (Laurent et al., 2006). Upon salinity challenge, fish remodel their gill morphology by changing the type of MRC and altering the abundance of ACs to regulate ion concentrations (Hwang and Lee, 2007). In killifish, the seawater gill type contains recessed MRCs and ion regulation in high salinities relies on the active excretion of Cl⁻ in MRCs and the passive diffusion of Na⁺ through leaky cell–cell junctions between MRCs and ACs (Wilson and Laurent, 2002). In contrast, the killifish freshwater gill comprises flat-surfaced MRCs with large external surface area (Katoh et al., 2001; Whitehead et al., 2012b). In fresh water, Na⁺ is actively imported but Cl⁻ is not; rather, plasma Cl⁻ concentration is maintained by preventing loss through tight cell–cell junctions (Laurent et al., 2006; Wood and Laurent, 2003) or by stimulating dietary Cl⁻ uptake (Wood et al., 2010). Scott et al. (2004) showed that the ability to transition between gill types is an integral aspect of acclimating to different osmotic conditions. Our results suggest that during this remodeling, the ability to transition between tight and leaky junctions may limit osmotic acclimation (Fig. 1). These data provide an intriguing area for future research and may suggest a trade-off in the regulation of tight junction permeability, where the ability to quickly transition to an impermeable tight junction in fresh water comes at the cost of efficient transition to leaky junctions in brackish water and vice versa.

Of the 4566 transcripts included in our analysis, we focused on three patterns of gene expression response during acclimation to osmotic challenge. The first set included genes showing a transcriptional response during acclimation that was common and parallel between both hypo- and hyper-osmotic challenges and between FW-native and BW-native populations (Fig. 2). These genes had a significant (P<0.01) effect of time only, no significant time-course interaction with challenge direction (hypo or hyper-osmotic) or population (268 transcripts fit this criterion and are detailed below). The second set of genes included 478 transcripts that showed expression responses that were different between hypo- and hyper-osmotic challenges, but conserved between populations (Fig. 3). These genes had significant time-course×challenge interaction, but no significant time-course×population or population×challenge interaction and are described below. Third, we identified 360 transcripts that are likely to be associated with differential acclimation abilities between populations (see below, Fig. 4). This set included genes that were significant for the three-way interaction (time-course×population×challenge) as well as genes with significant two-way interactions including the population term (time×population, population×challenge).

While exposure to hyper- versus hypo-osmotic environments requires different physiological responses to maintain homeostasis, a core set of compensatory responses is common to general osmotic stress (Bonga, 1997). Accordingly, we screened for genes in F. heteroclitus that showed a change in expression over time that was consistent during acclimation in both hypo- and hyper-osmotic challenges, and consistent between populations (referred to as ‘time only’ transcripts; Fig. 2). A caveat to consider with this set is how the experimental design may impact our ability to detect a common response between challenges. The challenges differ in their extent; the hypo-osmotic challenge spans a greater osmotic gradient than the hyper-osmotic challenge. This lack of congruence could result in different expression because of the extent rather than direction of exposure. Therefore, this analysis is conservative and identifies fundamental responses because it captures genes whose expression is consistent across both opposite and varying degrees of salinity challenge.

When this set of genes was analyzed together there was enrichment for ‘positive induction of I-kappaB kinase/NF-kappaB cascade’ [Gene Ontology database (GO) ID:0043123; P=0.028], ‘anatomical structure homeostasis’ (GO:0060249; P=0.067), ‘chromosome’ (GO:0005694; P=0.001) and ‘cell cycle’ (GO:0007049; P=0.027). These terms indicate a general stress response, particularly ‘positive induction of I-kappaB kinase/NF-kappaB cascade’. NF-κB is a family of transcription factors that coordinates immune response to infection, stress and injury (Ghosh et al., 1998; Xiao and Ghosh, 2005) and whose function is conserved across vertebrates (Kim et al., 2011). A sign test allowed us to segregate these genes into those that were up-regulated (117 transcripts), down-regulated (105 transcripts) and irregular through time (46 transcripts) (Fig. 2). The sets that were up-regulated and that were irregular through time were not significantly enriched for any functional categories (Fig. 2A,C). However, genes that were down-regulated were enriched for the terms ‘chromosome’ (GO:0005694, P<0.001) and ‘cell cycle’ (GO:0007049, P<0.001) (Fig. 2B). The differential expression of these terms and the associated genes indicates the regulation of cell division and growth. Altering the cell cycle is a common mechanism in response to stress because it allows sufficient time to repair DNA damage prior to replication (Kaufmann and Paules, 1996). Evidence for regulation of the cell cycle in response to osmotic stress has been identified in mammalian kidney cells (Kültz and Chakravarty, 2001), tilapia (Kammerer and Kültz, 2009) and leopard shark (Dowd et al., 2010). The data presented here suggest a similar change in regulation of the cell cycle, but additional studies are necessary to confirm this.

When transitioning between fresh and salt water, it is necessary to withstand and respond to the immediate challenges of osmotic stress (e.g. cell swelling or shrinkage) while simultaneously remodeling both the apical morphology of the gill as well as the basolateral or apical locations of transporters and ion channels (Fiol and Kültz, 2007; Marshall et al., 2002). Here, we highlight genes, organized by functional categories, which showed a different transcriptional response during acclimation between brackish and freshwater challenges, but were conserved in their response between populations.

The acclimation of F. heteroclitus between fresh and salt water is energetically costly, requiring 6–10% of the total energy budget (Kidder et al., 2006a,,b). For genes with salinity-specific responses (Fig. 3) there was significant (P<0.001) enrichment for GO terms associated with energy production, including ‘mitochondrial inner membrane’ (GO:0005743), ‘mitochondrion’ (GO:0005739), ‘generation of precursor metabolites and energy’ (GO:0006091) and ‘oxidative phosphorylation’ (GO:0006119). Two of these enriched terms – ‘generation of precursor metabolites and energy’ and ‘mitochondrion’ – were included in the set of genes showing up-regulation late in the hypo-osmotic treatment and late down-regulation in the hyper-osmotic treatment (Fig. 3A). This set of genes includes creatine kinase (CKM), which plays a major role in cellular energy regulation and has been implicated in acclimation to hypo-osmotic environments (Wallimann et al., 1992; Whitehead et al., 2012b). Also included are NADH dehydrogenase subunits, ATP synthase subunits and cytochrome B and C subunits, all of which are components of the mitochondrial oxidative phosphorylation energy production pathway (Newmeyer and Ferguson-Miller, 2003).

The difference in regulation of mitochondrion-related genes between treatments suggests that hypo-osmotic exposure requires elevated energy production relative to hyper-osmotic exposure (Fig. 3A). Kidder et al. (2006b) found that osmoregulation following seawater-to-freshwater transfer required less of the total energy budget of F. heteroclitus than freshwater-to-seawater transfers and that transfer to brackish water from either sea water or fresh water should be the least energetically costly of the three osmoregulation scenarios. Our results support this prediction; the up-regulation of mitochondrion-related genes following hypo-osmotic challenge and the concurrent down-regulation in hyper-osmotic challenge suggest that F. heteroclitus requires higher energy demands at 0.2 p.p.t. than at 15 p.p.t. (Fig. 3A). Alternatively, osmotic challenge requires switching between sea and freshwater MRC types, which requires mitochondrial biogenesis and therefore synthesis of proteins associated with energy production (Marshall et al., 2002; Scott, 2004). Thus, MRC production may result in the patterns observed here, though additional research would be needed to verify this finding.

Genes that were up-regulated in hypo-osmotic treatments but down-regulated in hyper-osmotic treatments were enriched for the GO terms ‘tight junction’ (GO:0005923, P=0.008), ‘integral to membrane’ (GO:0016021; P=0.009) and ‘ion binding’ (GO:0043167; P=0.03) (Fig. 3B). Permeability of the paracellular pathway in transporting epithelia is an important aspect of controlling ion homeostasis. The freshwater gill of F. heteroclitus actively accumulates Na⁺, but not Cl⁻, and relies on the regulation of epithelial tight junction proteins, such as claudins and occludins, to regulate the efflux of ions across the gill (Ernst et al., 1980; Patrick et al., 1997; Sardet et al., 1979). Claudins and occludins are involved in tight junction formation and play a role in paracellular ion transport (Furuse et al., 1993; Van Itallie and Anderson, 2006). Claudin 3, claudin 4 and occludin showed significant up-regulation in hypo-osmotic conditions coupled with a slight down-regulation under hyper-osmotic challenge (Fig. 3B, supplementary material Fig. S1). These genes were consistently up-regulated in response to hypo-osmotic challenge in multiple populations of F. heteroclitus (Whitehead et al., 2011,, 2012), as well as in other euryhaline Fundulus species (Whitehead et al., 2013), highlighting their importance for osmotic acclimation.

Genes involved in the active transport of ions showed regulation in response to salinity. Na-K-Cl co-transporter (NKCC1) and its co-localizing protein Focal adhesion kinase 1 (FAK1) trended towards down- and up-regulation in response to hypo- and hyper-osmotic challenges, respectively (Fig. 3C, supplementary material Fig. S1). NKCC1 is responsible for inward Cl⁻ transport across the basolateral membrane, after which Cl⁻ moves passively out of MRCs across the apical membrane via the cystic fibrosis transmembrane regulator (Hoffmann et al., 2007; Marshall et al., 2005). FAK co-localizes with NKCC1 and its phosphorylation activates NKCC1 (Hoffmann et al., 2007). We also observed significant up- and down-regulation of a Na⁺/K⁺-transporting ATPase (ATP1A1) in response to hypo- and hyper-osmotic challenge, respectively (Fig. 3A, supplementary material Fig. S1). Two V-type proton ATPase subunits (ATP6V1E1 and ATP6V0E1) showed a significant change in expression, but the directionality of this change was not consistent through time (Fig. 3A, supplementary material Fig. S1). V-type ATPases and ATP1A1 are located on the basolateral membrane of MRCs and are responsible for creating the negative electrochemical gradient necessary for uptake of Na⁺ in fresh water and thus up-regulation of these genes is expected under hypo-osmotic exposure (Marshall, 2002).

Genes putatively involved in calcium regulation are prevalent in this gene set and show up- and down-regulation in hypo- and hyper-osmotic exposure, respectively. These include otopetrin 1 (OTOP1), solute carrier family 25 (Slc25a25), calpain-1 (CAPN1), advillin (AVIL), gelsolin (GSN) and annexin A5 (ANXA5) (Fig. 3B, supplementary material Fig. S1). Intracellular levels of calcium play an important role in response to osmotic stress, because sensing of the external osmotic environment is achieved primarily though calcium binding and influx of intracellular calcium contributes to regulatory volume decrease (RVD) following swelling (Fiol and Kültz, 2007; Pasantes-Morales et al., 2006). Calcium binding by this set of proteins may respond to and regulate changes in intracellular ion concentrations. For example, OTOP1 regulates intracellular calcium concentrations (Hughes et al., 2007; Kim et al., 2010) and because an increase in intracellular calcium is essential in signaling RVD, it is possible that OTOP1 influences RVD following hypo-tonic swelling. This protein has been implicated in the divergence between F. heteroclitus and its congener, F. majalis; these species differ in their hypo-osmotic capabilities and OTOP1 may play a role in this physiological divergence (Whitehead et al., 2013). Although OTOP1 expression is associated with osmotic divergence at a macro-evolutionary scale (see Whitehead et al., 2013), microevolutionary divergence in OTOP1 expression is not yet evident within F. heteroclitus (supplementary material Fig. S1).

GSN, AVIL and ANXA5 mediate cellular responses to changes in calcium concentrations by binding and responding to free intracellular calcium (Fiol and Kültz, 2007; Gerke et al., 2005; Kothakota et al., 1997). These proteins are involved in actin modification and up-regulation of GSN and annexins related to ANXA5 has been shown to increase actin turnover rates and play a role in regulation of NKCC and gill remodeling under osmotic stress (Cunningham et al., 1991; Dos Remedios et al., 2003; Fiol and Kültz, 2007). AVIL has been implicated in the regulation of intracellular cytoskeletal organization and fibroblast formation, with its up-regulation indicative of cellular differentiation (Shibata et al., 2004). Taken together, these data suggest that regulation of the actin-based cytoskeleton is involved in the modification of MRC morphology when transitioning from high to low salinities and may be controlled by the interaction of calcium and the aforementioned calcium binding proteins (Daborn et al., 2001; Fiol and Kültz, 2007).

The physiological data are consistent with adaptive divergence in acclimation capabilities between the FW-native and BW-native populations. To explore the mechanisms involved in this divergence, we focused on 360 genes showing population-dependent divergence in expression during osmotic acclimation: that is, those genes that show any statistically significant population dependency (interaction) in their transcriptional responses during acclimation (Fig. 4). Despite acclimation to common lab conditions, some of the genes that differ between populations could be a product of developmental plasticity. However, the population-variable genes that we identify here probably underlie evolved differences between populations and provide insight into the evolutionary mechanisms behind adaptation to osmotic environments. This is because populations are genetically distinct and exhibit adaptive divergence in ionoregulatory physiologies (Whitehead et al., 2011) and almost no genes differ in expression between field-caught killifish, laboratory-acclimated animals and animals bred in the lab for 10 generations (Scott et al., 2009), indicating that most differences in gene expression between populations are heritable.

Transitioning between osmotic environments is energetically costly and the ability to efficiently manage these costs is likely to be important for successful acclimation. Within this gene set, we found enrichment for GO terms ‘mitochondrial component’ (GO:0044429; P<0.001), ‘cellular macromolecular complex subunit organization’ (GO:0034621; P=0.009) and ‘response to DNA damage stimulus’ (GO:0006974; P=0.03). Genes underlying the mitochondrial response, which consisted mainly of those involved in oxidative phosphorylation, increased in expression in the BW-native population in response to the low-salinity challenge. However, for these same mitochondrial genes, the FW-native population showed relatively little transcriptional response to low salinity (supplementary material Fig. S2). This result is consistent with the hypothesis that acclimation to low salinity environments is more energetically costly for the BW-native than the FW-native population. Additionally, genes contributing to enriched GO terms ‘macromolecular complex subunit organization’ and ‘response to DNA damage stimulus’ show a pattern of expression where differential transcription is observed in FW-native fish in response to hyper-osmotic challenge, but no response is seen in BW-native fish (supplementary material Fig. S2). Taken together, these patterns of expression divergence are consistent with the physiological data, indicating that each population is adapted to osmotic conditions within its native environment.

Cell volume regulation in the gills is one of the most immediate challenges of osmotic stress. Hypo-osmotic and hyper-osmotic environments can lead to cell volume increases and decreases, respectively. Cell swelling and shrinking can inhibit cellular functioning and growth and lead to cell death; RVD and regulatory volume increase (RVI) serve to maintain cell volume homeostasis (Hoffmann et al., 2009). We found population divergence in the regulation of genes that may be involved in RVD and RVI. Versican (VCAN), a protein that attracts cations and helps to retain cellular water (Faggian et al., 2007), showed significant decrease and increase in response to hypo- and hyper-osmotic challenge, respectively, in both populations. This gene trended towards a three-way time-course×population×challenge interaction because of a return to control values under hyper-osmotic challenge by the BW-native population. In contrast, no recovery to control levels of expression was observed in the FW-native population (P=0.019, supplementary material Fig. S1). Similarly, the water channel aquaporin 3 (AQP3) was up- and down-regulated in hypo- and hyper-osmotic exposures, respectively (supplementary material Fig. S1). BW-native fish again trend towards a return to control expression levels in the hyper-osmotic challenge while expression in the FW-native fish did not recover to control levels. Similar responses of AQP3 have been observed in previous studies of F. heteroclitus (Jung et al., 2012; Whitehead et al., 2012b), as well as other species of killifish (Kozak et al., 2014; Whitehead et al., 2013). Together with regulation of VCAN, these data suggest divergence in RVI and RVD responses between the BW-native and FW-native populations.

There is significant population divergence in the regulation of genes involved in polyamine synthesis including ornithine decarboxylase 1 (ODC1), spermine-binding protein (SBP), arginase 2 (ARG2) and transglutaminase 1 (TGM1) (supplementary material Fig. S1). ODC1 is the rate-limiting enzyme in polyamine synthesis that converts ornithine to putrescine (a polyamine) and has been implicated as playing a role in hypo-osmotic acclimation in mussels and brine shrimp (Lockwood and Somero, 2011; Rhee et al., 2007; Watts et al., 1996). Production of ornithine is controlled by ARG2 and thus, levels of ARG2 and ODC1 directly correlate to polyamine levels (Rhee et al., 2007). Conversely, SBP has an inverse relationship with polyamine levels as it binds free cellular polyamines, limiting their cellular availability (Tabet et al., 1993). TGM1 can also reduce free cellular polyamines because it facilitates protein–polyamine conjugation (Greenberg et al., 1991). The regulation of these genes indicates that synthesis of polyamines increases under hypo-osmotic but decreases under hyper-osmotic conditions. Additionally, these genes tend to exhibit transient transcriptional responses for BW-native fish during hyper-osmotic challenge. Polyamines may function to stabilize protein–nucleic acid interactions (PNAIs) during changes in osmolyte concentrations following fluctuations in cell volume (Capp et al., 1996). While osmolytes assist in maintaining PNAIs, under hypo-osmotic conditions osmolyte concentrations are decreased and polyamines are increased. Conversely, under hyper-osmotic conditions osmolytes are abundant and it may be energetically efficient to reduce polyamine synthesis. As osmolyte concentrations return to normal levels, we expect polyamine expression to return to baseline values. Alternatively, polyamines may play a role in regulatory volume decrease under hypo-osmotic stress (Watts et al., 1996). Increased expression under low salinity stress coupled with decreased expression under high salinity may work to maintain appropriate cell volume during acclimation. Transcriptional data from multiple experiments on Fundulus species (Whitehead et al., 2011,, 2012,, 2013) and on two species of mussel (Lockwood and Somero, 2011), indicate that regulation of the polyamine synthesis pathway is involved in osmotic acclimation, though the mechanistic role remains unknown.

An intriguing observation emerges in the set of genes that differ among the populations (Fig. 4). Fig. 4A shows the heat map of all 360 genes through the experiment, while Fig. 4B,C summarize the patterns of expression for this set of genes using principal components (PC) analysis. The first principal component (PC1, which accounts for 39.7% of expression variation) indicates a pattern of persistent change in expression for populations challenged with their non-native salinity, but transient changes in expression (expression returns to pre-transfer values) when populations are challenged with their native salinity. A similar pattern is observed in PC2 (20.1% of the expression variation), though it is much more pronounced for the BW-native population (Fig. 4C). These expression patterns are consistent with the hypothesis of ‘transcriptomic resilience’ as proposed by Franssen et al. (2011). The idea of transcriptomic resilience posits that the return of gene expression patterns to control levels (transient expression or transcriptional elasticity) is predictive of mechanisms that underpin adaptive physiological resilience to environmental stress (Franssen et al., 2011). That is, during exposure to native conditions to which organisms are adapted, osmotic perturbation is attenuated compared with that in non-native conditions and the transcriptional response follows suit. Here, for transcriptional patterns represented by PC1, we observe evidence for transcriptomic resilience, suggesting adaptation for both populations in their native environment (Fig. 4B). The loss of resilience for FW-native fish in PC2 is contrasted by the continued resilience of BW fish in hyper-osmotic conditions (Fig. 4C). This pattern is manifested in several physiologically relevant genes, including AQP3 and those involved in polyamine synthesis (supplementary material Fig. S1). The emergence of polyamine-synthesizing genes as a resilient set further strengthens the evidence that they serve an integral function in osmotic acclimation.

The principal component analysis of these population-dependent genes also suggests that salinity specific genes are more likely to be the targets of adaptive variation than those genes that show a common response to osmotic stress. In both PC1 and PC2, the response trajectories are opposite between hyper- and hypo-osmotic challenge (Fig. 4). This suggests that general responses to osmotic stress are not diverging between the populations, but rather that specific responses to either high or low salinity are driving population differences.

The divergence between populations in the magnitude and duration of the perturbation in plasma ion concentration suggest that proteins controlling paracellular permeability, particularly tight junction proteins (e.g. claudins, occludins), would be differentially regulated between populations. However, these genes do not show a pattern of population-dependent expression (supplementary material Fig. S1). A number of factors could explain this discrepancy. First, claudins are a large gene family (54 genes identified in Danio rerio, Baltzegar et al., 2013), but only five were included on our microarray, making it possible that claudins important for population-dependent acclimation were missing from the current analysis. Additionally, a recent study on two other Fundulus species shows that coding sequence and gene expression may evolve independently in response to osmotic environment (Kozak et al., 2014). Population divergence in tight junction physiology may be underpinned by protein variation rather than transcriptional variation. Last, claudin function can vary based on post-translational modifications such as phosphorylation, interactions with other claudins and interactions with scaffold proteins (Findley and Koval, 2009; Umeda et al., 2006). Ongoing RNA-seq experiments may offer more nuanced insight into the role of tight junction protein sequence and regulation in osmotic physiological divergence.

Since the FW-native and BW-native populations share an ancestor with a marine physiology (Duvernell et al., 2008), we conclude that the FW-native population has evolved expanded physiological abilities to accommodate low salinity stressors, consistent with previous findings (Whitehead et al., 2011). Concurrently, the FW-native population has lost some of the plasticity that enables acclimation to high salinities. Lande (2009) modeled how natural selection, following shifts to novel habitats, favors expansion of physiological plasticity in the short term, followed by genetic assimilation of that plasticity over the long term (∼3000 generations). The FW-native population likely invaded its current freshwater niche ∼15,000 generations ago (Duvernell et al., 2008). Thus this invasion may represent the assimilation phase of Lande's model and suggests that the FW-native population may have experienced an initial expansion of plasticity upon invading the freshwater habitat, where resilience to the saline end of the osmotic continuum has since been assimilated.

Studies on freshwater invasions by other groups agree with Lande's predictions for changes in plasticity because some recent invasions have shown expansions of plasticity. For example, Lee et al. (2011) found an expansion of plasticity after a marine to freshwater invasion by the copepod Eurytemora affinis. This habitat shift of E. affinis occurred only 150–300 generations ago, such that insufficient time may have accumulated for assimilation of plasticity. In contrast, older invasions may demonstrate the later phase where plasticity has been assimilated. For example, the marine or anadromous threespine stickleback, Gasterosteus aculeatus, has repeatedly and independently radiated into freshwater habitats throughout the northern hemisphere (Bell, 2001). Studies focusing on growth, reproduction and osmoregulation, among other measures, have found reduced survival, growth and overall health of the freshwater populations in sea water compared with ancestral marine or anadromous populations; a pattern that is consistent with assimilation of plasticity in the freshwater populations (DeFaveri and Merilä, 2014; Marchinko and Schluter, 2007; McCairns and Bernatchez, 2010). Research on salmonids (Morgan and Iwama, 1991), Alewives (Velotta et al., 2014) and Arctic char (Staurnes et al., 1992), among others, also show a reduction of plasticity after the invasion of a static freshwater environment.

We present evidence for mechanisms underlying population divergence within F. heteroclitus and recent studies on other species of killifish demonstrate that these mechanisms also underlie divergence at the macro-evolutionary level. For example, Kozak et al. (2014) found salinity-specific divergence between genes involved in cell junctions, V-type H⁺ ATPases and aquaporins between two sister species of killifish that occupy different osmotic niches. Whitehead et al. (2013) found divergence in expression between F. heteroclitus and the closely related, but predominantly marine, species F. majalis for genes involved in polyamine synthesis, cell volume regulation and aquaporins. This inter-species divergence matched the patterns described here between populations of F. heteroclitus, supporting the idea that micro-evolutionary processes are fine-tuning the same physiological pathways that underpin macro-evolutionary divergence in fundamental osmotic niches within the genus. Moreover, this intra-species divergence may be indicative of the first steps of a speciation event and representative of how osmotic adaptive divergence has proceeded in other species of fish.

The data presented here extends previous studies showing that FW-native and BW-native populations of F. heteroclitus have diverged according to the osmotic environment (Whitehead et al., 2011b). After expansion into freshwater niches, populations derived new osmotic compensatory phenotypes, where hypo-osmotic resilience was expanded, but hyper-osmotic resilience was reduced. Physiological data suggested that a trade-off in the ability to regulate freshwater versus brackish water tight junctions may play a role in limiting osmotic tolerance in F. heteroclitus; however, transcriptomic data did not confirm this pattern. Transcriptomic data provided insight into the mechanisms underlying general euryhaline phenotypes, particularly those genes involved in cell volume regulation, ion transport and energetic processes. Gene expression data also revealed genes that may be underlying evolved differences between the populations and suggest that the regulation of polyamines may play an essential role in this divergence. Additional studies focusing on sequence divergence combined with full genome expression profiles would help to elucidate the molecular mechanisms underlying the shifts in plasticity between these two populations. Moreover, our data suggest that the BW-native and FW-native populations represent the early stages in adaptive divergence across osmotic boundaries. Expanding this comparative study to other species within the genus Fundulus, including more distantly related species, could help identify the mechanisms repeatedly underlying diversification across osmotic boundaries.

Fish were sampled from two previously studied Fundulus heteroclitus (Linnaeus 1766) populations native to fresh water and brackish water along the Potomac River in the Chesapeake Bay region (Whitehead et al., 2011). BW-native adult fish were collected from Point Lookout State Park, near Scotland, MD (38°3′10.90″N, 76°19′34.38″W), which has an average 25 year salinity of 11.8 p.p.t., but can vary annually between ∼5.9 and 17.6 p.p.t. (www.chesapeakebay.net). FW-native adult fish were collected from Piscataway Park, near Accokeek, MD (38°41′42.18″N, 77°3′10.38″W) where salinity has not exceeded 0.12 p.p.t. since monitoring began in 1986 (www.chesapeakebay.net). All fish were housed at Louisiana State University in a recirculating system with biological, mechanical and UV filtration and acclimated for a minimum of 2 months at either 32 p.p.t. or 0.2 p.p.t. salinity (made from reverse-osmosis purified water mixed with Instant Ocean Artificial Sea Salt). Nitrogenous wastes (ammonia, nitrate, nitrite) were monitored using commercial test kits (API, Mars Fishcare Inc., Chalfont, PA) and remained non-detectable during experiments. Light cycles were 12 h:12 h light:dark and water temperature was maintained at 22°C.

For the hypo-osmotic challenge, fish from both populations were acclimated to 32 p.p.t. and subsequently transferred to 0.1 p.p.t. water. Five individuals were sampled pre-transfer and at each of the following time points after transfer: 6 h, 1 day, 3 days, 7 days, 14 days. For the hyper-osmotic challenge, fish acclimated to 0.2 p.p.t. and were transferred to brackish water (15 p.p.t.). Though there was no mortality after 32 p.p.t. to 0.1 p.p.t. transfer, range-finding experiments indicated that the reciprocal transfer from 0.1 p.p.t. to 32 p.p.t. caused limited mortality. As such, we chose 15 p.p.t. as the challenge for our hyper-osmotic transfer because that salinity crosses both the gill remodeling boundary and the isosmotic point, caused no mortality such that transcriptomic analysis would capture acclimation rather than morbidity responses and approximates the highest salinities found at Point Lookout State Park. Sampling (N=5) occurred pre-transfer and 6 h, 1 day, 3 days, 10 days, 17 days post-transfer. Control fish maintained at 32 p.p.t. and 0.2 p.p.t. were sampled at the same time points for both challenges. Blood was collected by caudal puncture using heparinized hematocrit capillary tubes and immediately spun at ∼1500 g for 1 min in order to isolate plasma. Plasma was stored at −20°C prior to sodium and chloride ion measurements. Fish were sacrificed by decapitation before gills were removed, preserved in RNAlater (Ambion) and stored at −20°C prior to RNA extraction for gene expression analysis.

Plasma sodium concentration was measured using flame atomic absorption spectroscopy and plasma chloride was measured using a modified mercuric thiocyanate method (Zall et al., 1956). Plasma data were analyzed using analysis of variance (ANOVA) with ‘time’, ‘salinity’ and ‘population’ as main effects with all interaction terms specified. Significance levels were set to 0.05 and post hoc comparisons were performed using Tukey's HSD method. All statistical analyses were performed in R (R Core Team, 2013). For the hypo-osmotic exposure, some of the plasma data have been reported previously (Whitehead et al., 2011).

Total RNA was extracted from gill tissue using TRIzol reagent (Invitrogen, Life Technologies Corp.) and antisense-RNA (aRNA) was prepared with the MessageAmp II aRNA Amplification Kit (Ambion, Life Technologies Corp.). The quality of total RNA and aRNA was assessed by microfluidic electrophoresis (Experion, Bio-Rad Laboratories). aRNA was coupled with Alexa Fluor dyes (Alexa Fluor 555 and 647, Molecular Probes, Inc.) and competitively hybridized to custom F. heteroclitus microarrays (Agilent Technologies; Whitehead et al., 2012a,b,, 2013) according to the manufacturers' protocols. Samples were hybridized in a randomized loop design with each challenge (i.e. hyper- and hypo-osmotic) processed in an independent loop. Each sample was hybridized twice including a dye swap. Slides were hybridized for 18 h at 60°C then washed and scanned using a Packard BioScience Scanarray Express (PerkinElmer, Waltham, MA, USA). Images were processed using ImaGene software (BioDiscovery, El Segundo, CA, USA).

Data processing, including normalization and statistical analyses, was performed in JMP Genomics (SAS Inc., Cary, NC, USA). Spots that were too bright (saturated fluorescence intensities) or too faint (<2 standard deviations above background fluorescence) were excluded from normalization and analysis. Remaining data were normalized using lowess followed by mixed model and quantile methods. Statistical analysis was performed using mixed models, where main effects included population (BW-native or FW-native), challenge (hypo-osmotic or hyper-osmotic) and time (pre-transfer and 0.25, 1, 3, 7, or 10 days post-transfer). Dye (Alexa 555 or 647) was specified as a fixed effect and array was specified as a random effect. Five biological replicates were included within each treatment and these replicates were specified as random effects. The statistical model was fully specified, including all interaction terms between main effects. Treatment effects were considered significant at P<0.01. Hierarchical clustering and principal component analysis (PCA) were performed in MEV (Saeed et al., 2006) and gene ontology enrichment analysis was performed using DAVID Bioinformatics Resources (Huang et al., 2008).

We thank Shujun Zhang, Charlotte Bodinier and Yanling Meng for assistance with osmotic challenge experiments and Jennifer Roach for assistance with transcriptomics data collection.