Resurrecting prehistoric parvalbumins to explore the evolution of thermal compensation in extant Antarctic fish parvalbumins

Most of our knowledge concerning the adaptation of proteins to varied thermal habitats comes from research on enzymes (Hochachka and Somero, 2002). According to the corresponding states theory, differences in the thermal sensitivity of Michaelis–Menten constants (K_m) of orthologous enzymes from poikilotherms adapted to different thermal regimes are driven by subtle changes in protein primary structure, which in turn lead to changes in conformational flexibility (Somero, 1978; Somero, 1983). These changes in conformational flexibility offset the function-altering effects of a change in temperature, thus conserving optimal function at physiological temperature. In thermodynamic terms, at lower temperatures less heat/energy is available to drive physiochemical reactions. Cold-adapted orthologs respond with a smaller net enthalpy change and a larger change in entropy than warm/temperate-adapted orthologs, which compensates for the reduced temperature/kinetic energy. This leads to a similar change in free energy and conserved binding ability at physiological temperatures (Somero, 1978; Somero, 1983; Somero, 1995; Jaenicke, 2000; Hochachka and Somero, 2002). Over evolutionary time scales, thermal compensation is accomplished through adjustments in protein primary structure. It appears that these compensating amino acid substitutions are generally excluded from the active site residues, which tend to be conserved to maintain substrate specificity (Wilks et al., 1988; Golding and Dean, 1998; Fields, 2001).

In prior work, the model derived from enzyme systems was extended to the non-catalytic Ca²⁺-binding protein parvalbumin (Erickson et al., 2005; Erickson and Moerland, 2006). Parvalbumins (PVs) are small (~10–12 kDa), acidic (pI~3–5) proteins that are unique to vertebrates and are present at high levels in the cytosol of fast-twitch skeletal muscle (Rall, 1996). The consensus view of PV function is that it acts as a soluble Ca²⁺ buffer in fast-twitch muscle cells. By allowing faster unloading of troponin C, the presence of PV leads to faster contraction/relaxation cycles (Rall, 1996). Isoforms of PV purified from white muscle of cold-adapted Antarctic fish of the Perciformes sub-order Notothenioidei and from temperate counterparts displayed thermal sensitivity patterns similar to those of enzymes of species adapted to different temperatures. Specifically, at a common measurement temperature Antarctic fish PVs showed a weaker binding affinity, as evidenced by a higher Ca²⁺ dissociation constant (K_d), than temperate counterparts, but at physiological temperatures function was highly conserved (Erickson et al., 2005; Erickson and Moerland, 2006). A specific structural mechanism leading to the thermal sensitivity pattern found in Antarctic fish PVs, however, has yet to be elucidated.

A member of the EF-hand family of cation-binding proteins, PV contains three EF-hand domains: two functional EF-hand ion-binding sites with high Ca²⁺ affinity and moderate Mg²⁺ affinity (the CD and EF sites; Fig. 1A) and a third (AB) that does not bind ions. Proteins of the EF-hand family consist of one or more pairs of the helix-loop-helix binding motif (Fig. 1). The Ca²⁺ coordinating residues are contained in a conserved 12 residue loop at relative positions 1, 3, 5, 7, 9 and 12 that correspond to a coordinate system: X, –X, Y, –Y, Z, –Z binding in a pentagonal bipyramidal arrangement (Fig. 1B,C) (Lewit-Bentley and Réty, 2000). The ion-binding domains, CD and EF, are highly conserved (Pauls et al., 1996) while the AB domain has a much more variable sequence and is considered to be the remnant of an ancestral binding site that has lost ion-binding function due to the loss of its paired EF-hand and two residues in its loop region (Cox et al., 1999).

The coordinating residues of PV are almost completely conserved from fish to humans. Analysis of an alignment of PVs from a variety of vertebrates reveals that a large proportion of substitutions are found in the AB domain away from the binding loops (A.C.W., unpublished observation). Evidence points to this region as being important for modulating PV affinity (Permyakov et al., 1991; Cox et al., 1999; Henzl et al., 2004). Recent work by Henzl and colleagues has shown that for rat PV, the ion-binding ability is quite sensitive to small sequence changes away from the active site. By altering the interactions and movements of the hydrophobic core, these substitutions affect the function of the highly conserved active sites (Henzl et al., 2008). The above evidence suggests that structural differences in the AB domain sequence, mediated through hydrophobic contacts in the protein core, can modulate affinity in the highly conserved functional binding sites. Thus, the overall ‘strategy’ observed in enzymatic systems in which adaptation is mediated by a small number of substitutions away from the ligand-binding site also appears to apply to non-catalytic proteins.

This study focuses on PV from fast-twitch skeletal muscle of Antarctic teleost fish of the suborder Notothenioidei. This group of Antarctic fish, which dominates the fish fauna in the Southern Ocean, has evolved from a benthic ancestor to fill a variety of niches in a unique environment (Eastman, 1993). The Antarctic clade of notothenioids first appeared around 35 million years ago with the subsequent adaptive radiation into the extant ‘species flock’ being mediated by the presence of antifreeze glycoproteins and a variety of other ecophysiological traits (Near et al., 2012). Water temperatures in this region have been below 5°C for the last 14 million years, and now the sub-zero temperature of the Southern Ocean is extremely stable (Eastman, 1993).

Typically, studies of thermal adaptation in enzymes have been conducted on proteins from closely related, often congeneric species that occupy different thermal niches (Holland et al., 1997; Fields and Somero, 1998). The polar and temperate species studied by Erickson et al. (Erickson et al., 2005) do not fall into the same study paradigm as the enzyme work. Instead of being congeneric or confamilial species, these polar (the notothenioids Gobionotothen gibberifron and Chaenocephalus aceratu) and temperate (carp Cyprinus carpio and bass Micropterus salmoides) species are highly divergent. Their last common teleost ancestor probably lived during the Triassic 200–250 million years ago (Patterson, 1993; Palmer, 1999; Nelson, 2006). Over this evolutionary distance neutral substitutions build up and can confound the typical simple sequence alignment analyses used to search for functionally adaptive substitutions (Bae and Phillips, 2004). In the present study, we sought to identify the amino acid substitutions that have occurred during the evolution of Antarctic fish that led to the characteristic polar thermal sensitivity pattern of notothenioid PV function.

We have used ancestral sequence reconstruction (ASR) as a tool to ‘walk’ backwards through sequence space to identify above the background of neutral substitutions the set of substitutions most likely to have caused a functionally significant shift in PV function during the evolution of Antarctic notothenioids in the frigid waters of the Southern Ocean. Originally described by Zuckerkandl and Pauling (Zuckerkandl and Pauling, 1965) and recently reviewed (Thornton, 2004), ASR has been used previously to successfully investigate environmental adaptation in an ecological context in a variety of protein systems (Thomson et al., 2005; Gaucher et al., 2008; Yokoyama et al., 2008).

Orthologs of PV corresponding to a series of ancestral fish between the last common ancestor of notothenioids and a temperate-adapted representative, carp, and the ancestral Antarctic notothenioid were reconstructed. Using three-dimensional structural modeling of these ancestral PVs we were able to track the changes that probably occurred during the evolution of PV in teleosts and discern which substitutions were most likely functionally neutral and which substitutions were most likely responsible for thermal adaptation. Using PV from the Antarctic fish G. gibberifrons as a template for site-directed mutagenesis, we resurrected the putative ancestral states and characterized the laboratory-generated ancestral PVs. Our findings support the parsimonious interpretation that these substitutions recapitulate the evolutionary trajectory of PV cold adaptation in Antarctic fish and are sufficient to explain the characteristic polar pattern of thermal sensitivity found for Antarctic notothenioid PVs (Erickson et al., 2005). Kinetic measurements of the rates of Ca²⁺ binding and dissociation for wild-type (WT) and mutant PVs, coupled with steady-state binding studies, suggest a structural mechanism for the cold adaptation of Antarctic fish PVs wherein ion binding is mediated by a subtle interplay of stabilization of the ion-bound and ion-free forms of PV.

Samples of Antarctic fish Chaenodraco wilsoni Regan 1914, Champsocephalus gunnari Lönnberg 1905, Chionodraco rastrospinosus DeWitt and Hureau 1979, Dissostichus mawsoni Norman 1937, Gobionotothen gibberifron Lönnberg 1905, Lepidonotothen nudifrons (Lönnberg, 1905), Notothenia coriiceps Richardson 1844, Notothenia rossii Richardson 1844, Parachaenicthys charcoti (Vaillant 1906), Patagonotothen ramsayi (Regan 1913), Pseudochaenichthys georgianus Norman 1937 and Trematomus hansoni Boulenger 1902 were collected by trawl or trap from the ARSV Laurence M. Gould in Dallman Bay, Antarctica in 2003. Dissostichus eleginoides Smitt 1898 samples were purchased from a local fish market. All fish were collected according to a protocol approved by the Animal Care and Use Committee of Florida State University (FSU) (Protocol no. 9304).

DNA sequencing of notothenioid fish full-length PV cDNA was performed using standard protocols (contact A.C.W. for detailed methods). DNA and protein sequence alignments were performed using Clustal W2 (Chenna et al., 2003) on the EBI webserver (http://www.ebi.ac.uk/Tools/msa/clustalw2/). An empirical Bayesian approach was used for ASR (Huelsenbeck and Bollback, 2001). A phylogenetic guide tree of 46 teleost fish PV sequences (supplementary material Table S1 gives accession numbers) was constructed with the topology conforming to known and well-supported evolutionary relationships (Briolay et al., 1998; Miya et al., 2003; Teletchea et al., 2006; Nelson, 2006; Li et al., 2007; Near and Cheng, 2008; Hertwig, 2008). Branch lengths were allowed to vary and were estimated from the data during the analysis. The ModelTest (Posada and Crandall, 1998; Posada, 2006) (http://darwin.uvigo.es/software/modeltest_server.html) and ProtTest (Abascal et al., 2005) (http://darwin.uvigo.es/software/prottest_server.html) servers were used to identify the appropriate models of nucleotide and amino acid sequence evolution, respectively. MrBayes v3.2 (Huelsenbeck and Ronquist, 2001; Ronquist and Huelsenbeck, 2003) was used for ASR of PV nucleotide sequences. Reconstruction of cDNA sequences was performed using the general time-reversible model with a gamma distribution of rates with a shape parameter, α, of 0.52311, and a proportion of invariable sites estimated by the program. Four chains were run for 1,000,000 generations with a sample frequency of 1000. The first 25% of the samples were discarded as after this period the four chains typically converged on a statistically equivalent set of parameters, which were used to estimate the consensus parameters. The ancestral sequences were then annotated from the stat files. Ambiguity in the reconstructions due to model choice and algorithm was investigated using the ASR program FASTML v.2.02 (Pupko et al., 2000; Pupko et al., 2002). Nucleotide sequences were reconstructed using the Jukes–Cantor model (Jukes and Cantor, 1969), the Yang codon model (Yang et al., 2000), Goldman–Yang codon model (Goldman and Yang, 1994) and an empirical codon model (Schneider et al., 2005). ProtTest identified the Jones–Taylor–Thornton amino acid substitution model (Jones et al., 1992) based on the PV alignment and this was also used in FASTML.

All PV homology models were constructed using the Swiss-Model webserver (Arnold et al., 2006) (http://swissmodel.expasy.org/). The crystal structure of C. carpio PV, pdb 4CPV, was used as the template for all models. Model validity was evaluated using the built-in methods of the Swiss-Model server. This includes the atomic empirical mean force potential ANOLEA which assesses packing quality of the model, and the GROMOS force field to assess the local quality of residues. Models were visualized with PyMol (The PyMOL Molecular Graphics System, Version 1.3, Schrödinger, LLC; www.pymol.org/) and VMD (Humphrey et al., 1996). Hydrogen atoms were added to the models in PyMol. Intermolecular distances were determined using the measurement wizard in PyMol. Protein-folding free energies were calculated in DeepView (Guex and Peitsch, 1997).

All cDNA cloning, site-directed mutagenesis, protein expression and purification were performed using established protocols (contact A.C.W. for detailed methods). Assay buffer and proteins were stripped of divalent cations using procedures modified from previously published protocols (Erickson et al., 2005; Erickson et al., 2006; Heffron and Moerland, 2008). To decalcify the assay buffer (20 mmol l⁻¹ Hepes, 150 mmol l⁻¹ KCl), a one-step procedure was used: 50 ml of assay buffer was mixed with 5% chelex-100 resin (Sigma Chemical, St Louis, MO, USA) in a 50 ml plastic conical tube covered with Parafilm (Pechiney Plastic Packaging Company, Chicago, IL, USA). Tubes were swirled on a rotary shaker overnight to chelate divalent cations. The tubes were then centrifuged in a table-top hanging-bucket centrifuge for 20 min at 1000 g to concentrate the chelex at the bottom of the tube. The buffer was then decanted into acid-washed 50 ml conical tubes, which were then covered with Parafilm. To adjust the pH of the assay buffer while avoiding contamination, small aliquots were removed and the pH measured. Then small aliquots of molecular biology grade, concentrated HCl (Sigma Chemical) were added to the stock buffer container to adjust the pH to 7.2.

Purified proteins were stripped of divalent cations using a four-step dialysis procedure. First, proteins were dialyzed against 1 l of assay buffer containing 5% chelex, 40 μmol l⁻¹ EDTA and 40 μmol l⁻¹ EGTA to scavenge cations, and 4 mol l⁻¹ urea to provide a denaturing environment. Second, the proteins were refolded by dilution in assay buffer with no urea containing 5% chelex, 40 μmol l⁻¹ EDTA and 40 μmol l⁻¹ EGTA. The third and fourth steps involved only assay buffer and 5% chelex. After dialysis, proteins were transferred to acid-washed plastic 15 ml conical tubes and stored at 4°C. Protein concentrations were then determined by the BCA assay (Pierce, Rockford, IL, USA). After determination of protein concentration, decalcified β-mercaptoethanol was added to the protein stocks to 14 mmol l⁻¹ to maintain a reducing environment.

A competition assay using the fluorescent Ca²⁺ indicator fluo-3 (Molecular Probes, Eugene, OR, USA) based on previously published methods (Eberhard and Erne, 1994; Erickson et al., 2005; Heffron and Moerland, 2008) was used to determine the Ca²⁺ K_d for all proteins. All titrations were performed using a Varian Cary Eclipse fluorescence spectrometer (Foster City, CA, USA) with internal temperature control. The excitation wavelength was 505 nm and the emission wavelength was 530 nm. Titrations were performed from 5 to 25°C in 5°C increments. All titrations were performed with β-mercaptoethanol added to 7 mmol l⁻¹ to maintain a reducing environment throughout the procedure.

Two rate constants define the interaction of Ca²⁺ with PVs – the off-rate or rate of dissociation (k_off) and the on-rate or rate of association (k_on). These rate constants provide mechanistic information about the nature of Ca²⁺ binding to PVs and the overall impact of temperature. Off-rates for Ca²⁺ were measured using terbium fluorescence as a reporter ligand based on established methods (Hou et al., 1991; Hou et al., 1992). PV has a higher affinity for terbium and the off-rate of Ca²⁺ can be measured as it is replaced by terbium and terbium fluorescence increases. Excitation of intrinsic phenylalanine residues in PV allows excitation of terbium through resonance energy transfer. Off-rates were measured using an Applied Photophysics (Surrey, UK) model SX.18 MV stopped-flow instrument housed in the laboratory of Dr Jonathan Davis at the Ohio State University Medical Center. The instrument has a mixing time of 1.4 ms at 15°C. PV samples were prepared as described above, including 1 mmol l⁻¹ DTT in the protein stocks. For off-rate measurements, 20 μmol l⁻¹ CaCl₂ was added to 10 μmol l⁻¹ PV. This solution was rapidly mixed with 500 μmol l⁻¹ TbCl₂ in the instrument and excited at 250 nm; emission was monitored through a 547 nm emission filter. Terbium fluorescence versus time curves were fitted with single and double exponential curves. A double exponential provided a better fit and these results are reported here. This analysis yields two k_off values – a fast k_off and a slow k_off. The fast k_off rates were used for data analyses and interpretation as these values correspond to previously published Ca²⁺ dissociation rate constants for other PV isoforms (Hou et al., 1992; Johnson et al., 1999; Lee et al., 2000). The second, slower rate constants are an order of magnitude slower than previously reported values. The Ca²⁺ association rate, k_on, was calculated using the relationship k_on=k_off/K_d.

A guide tree for ASR was constructed using the PV sequences of 46 teleost species (Fig. 3) including notothenioid PV sequences determined for this study (supplementary material Fig. S1). It should be noted that this tree only reflects topology or order of speciation. For visualization purposes the branch lengths in Fig. 3 were set to the arbitrary value of 1.

Sequences were estimated for five ancestors of modern, extant notothenioids (Fig. 4). This particular reconstruction scenario places the eel Anguilla japonica as the outgroup with Acanthomorpha containing the rest of the teleost orders represented in the tree. In order to discern the robustness of ASR, ambiguity in the reconstructions (defined as different possible reconstructions) was investigated by estimating sequences using different models of sequence evolution and a different program, FastML, which uses a different algorithm from MrBayes. A small amount of sequence ambiguity was found among the different models and programs compared with the GTR model used in MrBayes (supplementary material Fig. S2). The positions identified in this work as being functionally important (8 and 26) showed no ambiguity in the reconstruction.

Comparing the ancestral PV sequences reveals a set of substitutions between Perciformes ancestral PV (PAPV) and notothenioid ancestral PV (NAPV). This set of 16 substitutions correlates with paleontological and geological events highlighted in Fig. 3 including the cooling of the Southern Ocean to the current sub-zero temperatures. Of these 16 substitutions (Fig. 5), eight are conservative and eight are non-conservative.

In order to determine the subset of residues likely responsible for the observed thermal sensitivity profiles of extant notothenioid PVs and the potential structural and functional differences between the two ancestral PVs, homology models of both ancestral PVs were made using the major PV isoform from C. carpio, CCPV (PDB: 4CPV) (Kumar et al., 1990) as the template structure. This allows the substitutions to be viewed in a three-dimensional context. The sequence identities between the template and target structures are at least 78%. Models built with this level of sequence identity have been shown to be accurate and have resolutions comparable to that of crystal structures (Nayeem et al., 2006). Root mean square deviation (r.m.s.d.) of all the models compared with the template structure 4CPV, was within 0.05 Å, which is well within the resolution of the template structure of 1.5 Å. Additionally, evaluation by ANOLEA and GROMOS force fields within the Swiss-Model server showed a high level of quality of all the models.

In this study we focused on the eight non-conservative substitutions as these are most likely to affect function. Viewed in the tertiary structure, most of the substitutions reside outside of the active sites (supplementary material Fig. S3). One substitution at position 97 affects a Ca²⁺ coordinating residue. This residue, however, coordinates with its backbone oxygen and can accept a higher degree of variability than the other coordinating residues, which are almost completely conserved (Falke et al., 1994).

Virtual site-directed mutagenesis of G. gibberifrons parvalbumin (GBPV), an Antarctic fish PV which has been characterized functionally (Erickson et al., 2005), was used to determine the effect that each of the PAPV to NAPV substitutions had on the extant PV structure and potentially on function. In this process each substitution between PAPV and NAPV was substituted in reverse into the GBPV primary structure and then modeled to confirm that the virtual mutations would recapitulate the interactions present in PAPV (e.g. Fig. 6). Each of these models was evaluated as above.

Energies of unfolding of the homology models were calculated in DeepView using GROMOS force field calculations of free energy of folding as a representative measure of relative conformational flexibility/thermal stability (Table 1). Our results showed a correlation between this measure of relative stability and thermal habitat. The NAPV and GBPV had similar energies, which were both higher than those of the warm/temperate-adapted PAPV and the PV from the temperate representative (CCPV) and from the freshwater bass M. salmoides (MSPV). Most of the substitutions did not substantially affect the free energy of folding. The substitutions at position 8 and 26 individually reduced the energy of GBPV to approximately half-way towards that of the temperate/warm PVs. The double mutant (DM) PV had an energy value similar to that of the Perciformes PV structure, indicating that the PAPV structure, and potentially function, is recapitulated. Additionally, modeling of GBPV WT and DM show that the hydrogen bonds present in PAPV are introduced in the double mutant (Fig. 6). The complete conservation among all Antarctic notothenioid PVs (14 sequences) sampled here supports and underscores the importance of the evolutionary loss of these two hydrogen bonds in Antarctic fish PV as compensatory substitutions to the sub-zero temperatures of the Southern Ocean. Distance measurements in PyMol (supplementary material Fig. S4) show these substitutions to be located greater than 20 Å from the bound calcium in the ion-binding sites conforming to the paradigm seen in enzymes in which functionally modulating substitutions are located outside of ligand-binding sites (Fields, 2001; Hochachka and Somero, 2002).

The two ion-binding sites of PV have previously been shown to bind two Ca²⁺ molecules with no difference in affinity between sites (Pauls et al., 1993; Eberhard and Erne, 1994; Agah et al., 2003). This set of substitutions between the Periformes to notothenioid transition provides the search set for functionally adaptive substitutions. This has also been shown to be the case for Antarctic fish PV (Erickson et al., 2005). In the present effort, all PV titration curves were fitted with a hyperbolic function, which implicitly assumes single-site binding, i.e. independent binding of Ca²⁺ ion to each site.

Increased temperature produced a characteristic shift in the PV–Ca²⁺ titration curves (Fig. 7). This was apparent from the decreased steepness of the hyperbolic plots, which reflects weaker Ca²⁺ binding at higher temperatures. All four protein constructs showed this same general trend of decreased binding ability, as evidenced by a larger K_d value with increasing temperature, which is indicative of an exothermic process (Fig. 8). At each measurement temperature the WT protein gave the highest K_d. All three mutant PVs showed a right shift in thermal sensitivity pattern (Fig. 8), indicating increased Ca²⁺-binding ability. The WT PV construct showed a sharper increase in K_d above 20°C than the three mutants, indicating that it had a higher thermal sensitivity.

Two-way analysis of variance (ANOVA) was used to test for significant differences among the K_d values with respect to measurement temperature. For the ANOVA tests, K_d was the measurement variable and temperature and PV construct were the nominal variables. All K_d values were found to be significantly different for the four proteins (P<0.0001) at each measurement temperature except for the comparison of the K8N and K26N constructs, for which no significant difference was found (P=0.8647), indicating that K8N and K26N PVs show equivalent Ca²⁺-binding ability across measurement temperatures.

A representative raw data trace for the stopped-flow measurements is shown in supplementary material Fig. S5. While the two EF hand ion-binding sites of PV were assumed to be equivalent in analysis of K data, Ca²⁺ d dissociation rates using terbium fluorescence as a reporter displayed a better fit using a double exponential regression. As stated in Materials and methods, we report the faster component for k_off (Fig. 9). With the exception of the K26N construct at 25°C, all three mutants showed a right shift in k_off versus temperature curves compared with the WT construct. The three mutants showed a nearly identical shift in off-rates at each measurement temperature. Ca²⁺ on-rates, which were calculated as described above, are shown in Table 2. Because k_on values are derived from the ratio of two means of different size N values, there is no variance for the estimates of on-rates. Calculated k_on values for the WT, K8N and K26N constructs were virtually identical at corresponding temperatures. In contrast, k_on values for the DM PV were consistently higher than values for the other PV constructs. Note that the k_on value for the double mutant at 25°C was nearly identical to the corresponding value for the WT construct at 5°C.

Antarctica and its surrounding waters have been isolated for around 25 million years since the opening of the Drake Passage and the formation of the Antarctic Circumpolar Front (Eastman, 1993). Water temperatures in this region have been below 5°C for the last 14 million years, and remain very stable today (Eastman, 1993). On the Antarctic Peninsula, temperatures range from −1.8°C in winter to only 1.5°C in the summer, while on the Ross Ice Shelf the sea temperature is a constant −1.86°C year round (Eastman, 1993). Notothenioid fishes are extreme stenotherms that have developed a suite of adaptations for living in the frigid waters of their habitat, including antifreeze glycoproteins and metabolic enzymes tuned to function optimally in these extremely cold waters (Sidell, 2000; Coppes Petricorena and Somero, 2007). The diverse and endemic notothenioids provide an excellent study system to investigate thermal adaptation at the protein level.

It has been shown in several groups of enzyme orthologs that protein function correlates tightly with environmental temperature, which is thought to be due to adaptive amino acid substitutions that maintain protein function at an optimal level at physiological temperatures (Holland et al., 1997; Fields and Somero, 1998; Dong and Somero, 2009). The corresponding states theory (Somero, 1978; Somero, 1983; Somero, 1995; Jaenicke, 2000; Fields, 2001; Hochachka and Somero, 2002) predicts that differentially adapted orthologs (homologous protein isoforms separated by speciation) existing as populations of a series of related conformational states will sample a similar subset of states at their physiological temperature. This similarity in conformational states gives proteins the appropriate amount of flexibility to provide the necessary level of function in an organism's thermal habitat.

It is important to note that conformational flexibility can be considered a phenomenon distinct from thermal stability. Thermal stability, or macro-stability, is a protein's resistance to denaturation in the face of temperature or chemical denaturants. Thermal stability maintains the integrity of protein tertiary structure. Conformational flexibility, or micro-stability, is an inclusive term that describes a variety of protein movements including but not limited to global peptide movements or oscillations termed ‘breathing motions’, loop movements and transiently unfolded local regions – all of which affect the rigidity of the folded protein structure (Privalov and Tsalkova, 1979; Vihinen, 1987; Závodszky et al., 1998; Fields, 2001). While both correlate with an organism's physiological temperature, it has been shown experimentally through directed evolution studies that these two parameters can in fact be decoupled (Miyazaki et al., 2000; Wintrode et al., 2000). For proteins from cold-adapted organisms, the characteristic thermolability may be a result of a lack of selection for thermal stability as these proteins never experience temperatures that would cause denaturation (Fields, 2001). Changes in conformational flexibility, however, could be selected for maintenance of optimal protein function at environmental temperature (Fields and Somero, 1998; Fields et al., 2001).

With this framework in mind, we can make predictions about the opposing effects of temperature and adaptive substitutions on PV function in Antarctic notothenioid fish. Because of the exothermic nature of PV Ca²⁺ binding, decreased temperature would cause a decrease in Ca²⁺ K_d. In white muscle this increase in PV Ca²⁺ affinity would disrupt the correct timing of muscle contraction by diverting Ca²⁺ from binding to troponin C, which presumably could put fish at a disadvantage when trying to feed on swimming prey or when evading predators. Indeed, it has recently been shown that zebrafish swimming performance is correlated with PV content (Seebacher and Walter, 2012), indicating a role of PV in ecological fitness of teleost fish. For Antarctic fish PVs to function optimally we would expect to find evolutionary substitutions of amino acid residues compensating for the decreased ambient temperature. Such changes would have occurred sometime after the Antarctic notothenioids diverged from their temperate ancestor.

In addition to conservation of Ca²⁺-binding ability (Ca²⁺ K_d), Erickson and colleagues (Erickson et al., 2005) found that the cold- and temperate-adapted PV isoforms studied showed similar changes in the free energy of binding (ΔG). Overall, across all temperatures tested using isothermal titration calorimetry they found that the Antarctic isoforms displayed a less negative, i.e. less favorable, enthalpy change (ΔH), but had a more positive, i.e. more favorable, entropy change (ΔS) than the temperate isoforms, while ΔG remained relatively stable. In functional terms, as temperature decreases an exergonic reaction (ΔG<0) will be favored and in the case of PV will be manifested as increased binding ability or lower K_d. Indeed, the PV isoform from the temperate-adapted carp studied by Erickson and colleagues (Erickson et al., 2005), which has a more negative enthalpy change than the cold-adapted notothenioid PVs, has a lower K_d at all temperatures measured. The larger positive change in entropy (indicative of increased flexibility) may compensate for the less negative enthalpy change in notothenioid PVs, which experience an environment with less thermal energy to drive bond breaking and formation.

We used ASR, homology modeling and free energy of folding calculations to identify specific amino acid substitutions that could account for cold adaptation in notothenoid fish PVs. A composite phylogeny reflecting currently understood teleost systematics was constructed for ASR. This tree gives a broad sampling of teleost species. Nothenioidei are a sub-order of the large and diverse order of Perciformes. This area of the phylogenetic tree is of most interest for reconstruction of notothenioid ancestral sequences and this region of the tree is the most densely sampled. As a counter example, discerning the adaptive substitutions in PV from the polar cod, Borerogadus saida, would require a dense sampling of PV sequences from Gadiformes fish. By comparing the reconstructed PV sequences for the cold-adapted notothenioid ancestor and several of its teleost ancestors (Fig. 4) it was possible to track the amino acid substitutions that occurred during PV evolution along a trajectory towards adaptation to the cold Southern Ocean. Additionally, viewing these substitutions in a phylogenetic context allows the relative timing of the substitutions to be visualized.

With ASR we have attempted to ensure that we are comparing orthologs in our analysis. The first step was having a data set with only β-type PVs. A sequence alignment of 164 PV sequences from all classes of jawed vertebrates shows that α- and β-type PVs are distinguished by residues at position 19 and 67. Specifically, α-type PVs have a typical phenylalanine at position 19 and a small hydrophobic residue, valine or isoleucine, at position 67. β-Type PVs typically display a cysteine at position 19 and a phenylalanine at position 67 (A.C.W., unpublished observation). In addition to selecting only β-type PVs, we selected sequences from GenBank described to have been the major isoform isolated from the white muscle of teleost fish. While this does not guarantee the syntenic relationship of the analyzed PVs, it does give us a high degree of confidence that we are comparing orthologs. It should also be noted that many teleosts carry multiple genes for PV and express more than one β-type isoform in their white muscle. This opens up a possible alternative evolutionary scenario where an ancestral fish carried multiple isoforms of PV, each with a different thermal sensitivity pattern and differential expression based on habitat. This scenario assumes that being eurythermal is the ancestral condition for the Notothenioidei and the ancestral notothenioid already had the cold-adapted PV isoform prior to the extreme cooling of the Southern Ocean. Then, after cooling of the Southern Ocean the ancestral notothenioid lost expression of the other isoforms and now only expresses the cold-adapted isoform. A full analysis of the syntenic relationships among the somewhat complex PV phylogeny is beyond the scope of this work and alternative evolutionary scenarios do not refute our parsimonious interpretation: that substitutions identified by ASR and homology modeling are sufficient to explain the current thermal sensitivity pattern of Antarctic fish PVs.

Table 2 gives the predicted effect of each of these substitutions when viewed in a three-dimensional context. The substitutions at positions 8 and 26 stand out as being important as they introduce hydrogen bonds into the GBPV model. Hydrogen bonds have a negative enthalpy of formation and are stabilized by decreased temperature. Additionally, the presence of solvent-exposed, charged lysine residues in the cold-adapted PVs could cause a slight increase in flexibility compared with the solvent-exposed asparagines in the warm/temperate PVs. A shift from uncharged asparagines to positively charged lysines may destabilize protein structure through enhanced solvent interactions (Feller et al., 1997; Fields, 2001). Our phylogenetic and modeling results led to our hypothesis that reintroduction of the hydrogen bonds individually will provide an intermediate, but not necessarily additive decrease in conformational flexibility, and an intermediate shift in thermal sensitivity pattern. Specifically, the single mutants (K8N and K26N) would show tighter Ca²⁺ binding (measured as K_d) than the wild-type GBPV at each measurement temperature. Moreover, we proposed that the DM form would show a full conversion from the Antarctic WT thermal sensitivity pattern to that of the temperate counterparts. Essentially, this would represent the evolutionary steps that occurred, leading to the current thermal sensitivity pattern of Ca²⁺ binding displayed by Antarctic notothenioid PVs.

The steady-state Ca²⁺-binding measurements reported here support this hypothesis. At a common measurement temperature the single mutants showed a stronger binding ability than the WT. However, the Ca²⁺ K_d of the single mutants measured at 15°C was virtually identical to the value obtained for the WT protein at 5°C, suggesting conservation of function. The single mutants could correspond to a phenotypic state present in an ancestral, transitional fish inhabiting waters at an intermediate temperature. This cool/temperate-adapted ancestral fish could have had a PV with either of these single substitutions and maintained an optimal level of Ca²⁺ buffering at this period of intermediate cooling of the Southern Ocean.

Introduction of both hydrogen bonds at positions 8 and 26, creating the DM variant, caused a full conversion of binding ability to that of a warm/temperate-adapted PV. Again, there is conservation of function at proposed physiological temperatures. The DM displays a nearly identical K_d value at 25°C to that of the WT at 5°C. In effect, this reverse engineering of the warm-adapted state supports the hypotheses regarding the evolutionary steps toward thermal compensation in Antarctic fish PV. The evolutionary loss of two hydrogen bonds during the cooling of the Southern Ocean is sufficient to establish the current thermal sensitivity pattern of PV function in extant Antarctic fishes.

In addition to altered functional sensitivity, our K_d measurements suggest that the WT PV displays a lower thermal stability than the mutant forms. The WT K_d profile shows a sharp increase in K_d at the highest measurement temperature. Also, estimates of the free energy of folding correlate not only with habitat temperature (known for extant isoforms and proposed for recombinant proteins) but also with thermal sensitivity of function. Interestingly, these results suggest that in the case of Antarctic fish PVs, thermal stability and conformational flexibility (here inferred from changes in functional sensitivity) are directly linked. Combined with thermodynamic data (Erickson et al., 2005), the above functional and thermal sensitivity data suggest that altered flexibility is the driving force behind the changes in the thermal sensitivity profile of the mutant PVs. With the WT being equivalent to the muscle-purified G. gibberifrons PV, and the DM representing a warm/temperate-adapted PV isoform, it follows that the WT should have a smaller enthalpy change and a larger entropy change associated with Ca²⁺ binding than the DM. This is seen qualitatively from van't Hoff plots of the Ca²⁺ K_d data (supplementary material Fig. S6, Table S2). A more flexible structure for the WT is inferred from this larger entropy change. This inference is supported by the higher thermal sensitivity of the WT revealed by the functional data. Furthermore, the single mutants should have an intermediate shift in enthalpy and entropy change, and an intermediate level of flexibility, corresponding to the intermediate thermal sensitivity pattern displayed for Ca²⁺ K_d measurements. Again, this is supported qualitatively by calculations of ΔH and ΔS from the van't Hoff plots. More direct studies of conformational flexibility in WT and mutants through either simulation studies or empirical structural measurements would be needed to show that the observed changes in functional sensitivity to temperature in fact correlate with changes in conformational flexibility.

In order to obtain a more comprehensive model of cold adaptation of Antarctic fish PVs, we sought to further characterize mechanistic aspects of thermal compensation in Ca²⁺-binding ability. In addition to the steady-state Ca²⁺ K_d data, we determined the unidirectional rate constants of Ca²⁺ binding and dissociation. These two rate constants define the interaction of Ca²⁺ with PV and it has been shown that these parameters can directly affect muscle function (Hou et al., 1991; Hou et al., 1992). The equilibrium Ca²⁺-binding constant K_d is determined by the ratio of off-rate over on-rate (K_d=k_off/k_on). These rate constants provide mechanistic information about the nature of Ca²⁺ binding to PVs and the overall impact of temperature. The off-rates or rates of dissociation (given by k) describe the rate at which Ca²⁺ off leaves the ion-binding sites of PV. From this metric, we can infer the relative stability of the Ca²⁺-bound state. A slower k_off is indicative of a stabilized bound state. The other metric, k or Ca²⁺ on -binding rate, describes the rate of association of Ca²⁺ with PV. The rate at which Ca²⁺ loads onto the PV ion-binding site is indicative of the stability of the apo or unbound state of PV. A faster k_on suggests a loop structure that is closer to the bound state loop structure. Henzl and colleagues have demonstrated through NMR solution structures of rat PV isoforms that the structural similarity of the bound and apo-state binding loops (i.e. loop rigidity) correlates with Ca²⁺-binding affinity (Henzl and Tanner, 2007; Henzl and Tanner, 2008). It follows that a more stable apo state suggests a more rigid binding loop.

The combined steady-state and kinetic data show an interesting pattern. The K8N and K26N single mutant constructs had generally slower off-rates than the WT construct (Fig. 9), but very similar on-rates (Table 2). This suggests that the differences in K_d between the single mutants and the WT are mediated by a stabilization of the Ca²⁺-bound state by introduction of a single hydrogen bond in the AB domain without concomitant stabilization of the apo state. The DM construct, in contrast, did not show a further decrease in off-rate but had consistently higher k_on rates at all temperatures when compared with the WT and single mutant constructs, indicating that the further increase in binding affinity seen with both substitutions is mediated through a faster on-rate. This implies that the presence of two hydrogen bonds at a distance from the ion-binding sites potentially stabilizes the apo state of the DM construct. Another possible explanation for the changing on-rates, which are nearly diffusion limited, is a difference in surface charge in PV. The substitutions at position 8 and 26 change a positively charged residue, lysine, to a non-charged asparagines, which would make the overall charge of the PV molecule more negative, allowing a stronger charge interaction with the positively charged Ca²⁺ ion, thus increasing the on-rate.

The present study provides a possible structural and mechanistic basis for the thermal sensitivity pattern of PV isoforms found in Antarctic notothenioid fish. Changes in the structure of the AB domain of PV, distant from the functional ion-binding sites, modulate Ca²⁺-binding ability, presumably mediated through intramolecular contacts in the hydrophobic core. Decreased flexibility in the AB domain leads to stabilization of the binding loop either in the apo or bound state. Interestingly, the stabilization of the apo and bound state of PV appears to be decoupled. The reintroduction of a single hydrogen bond stabilizes the bound state while reintroduction of both hydrogen bonds leads to a stabilization of the bound and apo states.

In the context of the evolving ancestral notothenioid, as the Southern Ocean began to cool, one of two hydrogen bonds was lost from PV by an asparagine to lysine substitution at either position 8 or 26, leading to a destabilization of the PV apo state and a left shift in the thermal sensitivity pattern of Ca²⁺-binding ability to an intermediate stage. Function would have been conserved in this transitional ancestral PV (represented by K8N and K26N) found in an ancient notothenioid fish inhabiting cool/temperate waters. As the Southern Ocean continued cooling to the present frigid temperatures, the loss of a second hydrogen bond would have destabilized the bound state of PV, leading to the further left shift in the thermal sensitivity pattern of Ca²⁺ binding seen for the extant Antarctic fish PV, represented here by the WT construct.

We would like to thank Danielle Sandoz-Osmus for assistance during the sequencing phase of the project. We would also like to thank Dr Jonathan Davis and the Davis Lab at The Ohio State University Medical Center for training and assistance with stopped-flow spectrometry measurements. Dr P. Bryant Chase provided helpful comments on the manuscript. Dr Ross Ellington provided invaluable assistance at all stages of this project. The late Dr Bruce Sidell provided critical input in the initial stages of this project and made possible specimen collection through a National Science Foundation (NSF) Office of Polar Programs grant.

 
• ASR

  ancestral sequence reconstruction

 
• CCPV

  Cyprinus carpio parvalbumin

 
• DM

  double mutant

 
• K8N

  position 8 lysine to asparagine substituted parvalbumin

 
• K26N

  position 26 lysine to asparagine substituted parvalbumin

 
• MSPV

  Micropterus salmoides parvalbumin

 
• NAPV

  notothenioid ancestral parvalbumin

 
• PAPV

  Perciformes ancestral parvalbumin

 
• PV

  parvalbumin

 
• WT

  wild-type Gobionotothen gibberifrons parvalbumin