Epistasis constrains mutational pathways of hemoglobin adaptation in high-altitude pikas.

A fundamental question in evolutionary genetics concerns the roles of mutational pleiotropy and epistasis in shaping trajectories of protein evolution. This question can be addressed most directly by using site-directed mutagenesis to explore the mutational landscape of protein function in experimentally defined regions of sequence space. Here, we evaluate how pleiotropic trade-offs and epistatic interactions influence the accessibility of alternative mutational pathways during the adaptive evolution of hemoglobin (Hb) function in high-altitude pikas (Mammalia: Lagomorpha). By combining ancestral protein resurrection with a combinatorial protein-engineering approach, we examined the functional effects of sequential mutational steps in all possible pathways that produced an increased Hb-O2 affinity. These experiments revealed that the effects of mutations on Hb-O2 affinity are highly dependent on the temporal order in which they occur: Each of three β-chain substitutions produced a significant increase in Hb-O2 affinity on the ancestral genetic background, but two of these substitutions produced opposite effects when they occurred as later steps in the pathway. The experiments revealed pervasive epistasis for Hb-O2 affinity, but affinity-altering mutations produced no significant pleiotropic trade-offs. These results provide insights into the properties of adaptive substitutions in naturally evolved proteins and suggest that the accessibility of alternative mutational pathways may be more strongly constrained by sign epistasis for positively selected biochemical phenotypes than by antagonistic pleiotropy.

Introduction

The adaptive evolution of protein function may typically involve the sequential fixation of individual amino acid mutations with each substitution producing an incremental improvement in the selected biochemical phenotype. However, there are two interrelated factors that can potentially affect fixation probabilities and constrain the number of selectively accessible pathways to high-fitness genotypes: Mutational pleiotropy (where a given mutation perturbs multiple phenotypes) and epistasis (nonadditive interactions between mutations). As amino acid mutations often affect multiple aspects of protein biochemistry (Taverna and Goldstein 2002; Wang et al. 2002; DePristo et al. 2005; Bloom et al. 2006; Weinreich et al. 2006; Tokuriki et al. 2008; Tokuriki and Tawfik 2009), mutations that improve one aspect of protein function may simultaneously compromise other functions. Pleiotropic trade-offs can also give rise to a context-dependence of mutational effects, as the phenotypic effect of a given mutation may be contingent on the existence of compensatory mutations at other sites in the same protein.

During the adaptive evolution of a particular protein function, epistasis can affect the fixation probabilities of function-altering mutations and influence trajectories of protein evolution in two ways: 1) The sign or magnitude of the mutation’s effect on the selected phenotype may be conditional on genetic background, and 2) the sign or magnitude of the mutation’s pleiotropic effects may be conditional on genetic background. Both forms of epistatic interaction can influence whether evolution is more likely to follow some pathways rather than others (Weinreich et al. 2006; DePristo et al. 2007; Ortlund et al. 2007; Poelwijk et al. 2007; Bridgham et al. 2009; Lozovsky et al. 2009; Carneiro and Hartl 2010; da Silva et al. 2010; Kryazhimskiy et al. 2011, 2014; Kvitek and Sherlock 2011; Rokyta et al. 2011; Salverda et al. 2011; Ostman et al. 2012; Weinreich and Knies 2013; de Visser and Krug 2014; Harms and Thornton 2014).

A compelling system for investigating the pleiotropic and epistatic effects of function-altering mutations is provided by the evolution of increased hemoglobin (Hb)-O2 affinities in vertebrate species that have adapted to environmental hypoxia. Among vertebrates, elevated Hb–O2 affinities have evolved in multiple species that have adapted to chronic O2 deprivation at high altitude (Storz 2007; Weber 2007; Storz and Moriyama 2008; Storz et al. 2009; Storz, Runck et al. 2010; Natarajan et al. 2013; Projecto-Garcia et al. 2013; Revsbech et al. 2013). Under severe hypoxia, an increased Hb–O2 affinity can help preserve an adequate level of tissue oxygenation by enhancing pulmonary O2 uptake while simultaneously maintaining the pressure gradient that drives O2 diffusion from capillary blood to the tissue mitochondria (reviewed by Storz, Scott et al. 2010; Mairbäurl and Weber 2012). Even slight increases in Hb–O2 affinity can efficiently safeguard arterial O2 saturation under hypoxia, thereby complementing numerous other physiological adjustments in the cardiopulmonary system (Turek and Kreuzer 1976; Turek, Kreuzer, Ringnalda 1978; Turek, Kreuzer, Turek-Maischeider, et al. 1978; Mairbäurl 1994).

Mammalian Hb is a heterotetramer, consisting of two α-type subunits and two β-type subunits that each contain a heme group—a porphyrin ring with a ferrous (Fe2+) iron atom capable of reversibly binding a single dioxygen molecule. The Hb tetramer undergoes an oxygenation-linked transition in quaternary structure, whereby the two semirigid α1β1 and α2β2 dimers rotate around each other by 15° during the transition between the deoxy (T-state) and oxy (R-state) conformations (Perutz 1972, 1989; Baldwin and Chothia 1979). This oxygenation-linked shift in quaternary structure between the T- and R-states is central to the allosteric function of Hb as an O2-transport molecule, and underpins the subunit cooperativity of O2-binding as well as the sensitivity of O2-affinity to allosteric effectors such as protons, Cl- ions, and organic phosphates. Allosteric effectors have the effect of reducing Hb–O2 affinity because they preferentially bind to deoxyHb, thereby stabilizing the low-affinity T structure through the formation of additional hydrogen-bonds and salt bridges within and between subunits (Bettati et al. 1983; Perutz 1989).

There are numerous possible mutational changes that can alter Hb–O2 affinity, but many or most such changes have deleterious pleiotropic effects on other structural or functional properties (Bellelli et al. 2006; Bonaventura et al. 2013; Varnado et al. 2013). For example, active site mutations that alter the polarity or hydrophobicity of the distal heme pocket can produce significant changes in ligand affinity, but such mutations can reduce structural stability and/or increase rates of spontaneous heme oxidation (autoxidation), which renders the oxidized Fe3+ “metHb” incapable of O2-binding (Varnado et al. 2013). Similarly, mutations that increase Hb–O2 affinity by suppressing sensitivity to allosteric effectors sacrifice an important mechanism of phenotypic plasticity because changes in intraerythrocytic concentrations of organic phosphates and other effectors are rendered less effective as a means of modulating blood–O2 affinity in response to transient changes in metabolic demand or environmental O2 availability. Thus, mutations with identical main effects on Hb–O2 affinity may still be unequal in the eyes of natural selection if they vary in the magnitude of antagonistic pleiotropy (Otto 2004; Streisfeld and Rausher 2011).

To assess the roles of pleiotropy and epistasis in the evolution of protein function, we used a combinatorial protein engineering approach to dissect the mechanistic basis of evolved differences in Hb–O2 affinity between two sister species of North American pikas (Mammalia: Lagomorpha) that have contrasting altitudinal range limits. Our comparison included the American pika (Ochotona princeps), an alpine specialist, and the closely related collared pika (O. collaris), which has a predominantly lowland distribution in Alaska and northwestern Canada. Ochotona princeps inhabits the highest summits in the contiguous United States at elevations of approximately 4,300 m where the partial pressure of O2 (PO2) is 60% of the sea level value (∼96 vs. ∼160 torr). Pikas originated and diversified in the Palearctic (Erbajeva 1994), and O. princeps and O. collaris represent the only two extant pika species in North America. The two species diverged approximately 4–8 Ma (Lanier and Olson 2009). Although O. princeps is mainly found above timberline in high-altitude talus habitats (Smith and Weston 1990; Hafner 1994), fossils and biogeographic evidence indicate that the species also inhabited lower elevations during Quaternary glacial cycles (Mead 1987; Hafner 1993; Galbreath et al. 2009, 2010; Lanier and Olson 2013). Because of the dramatic difference in elevational range limits between O. princeps and O. collaris, comparisons between these two closely related species provide an opportunity to identify the molecular basis of evolved differences in Hb function that may have adaptive significance.

Using ancestral protein resurrection and site-directed mutagenesis, we measured the structural and functional effects of sequential mutations in all possible pathways connecting the resurrected ancestral (low-affinity) Hb and the derived (high-affinity) Hb of the high-altitude species. In addition to measuring oxygenation properties of each engineered recombinant Hb (rHb) mutant, we also measured pleiotropic effects on functional properties that can potentially trade-off with Hb–O2 affinity.

The experiments were designed to answer the following questions: 1) What fraction of the evolved change in Hb–O2 affinity is attributable to substitutions at each step in the pathway? 2) Are the functional effects of mutations conditional on genetic background? For example, do the functional effects of individual mutations depend on the sequential order in which they occur? And 3) Do affinity-enhancing mutations have deleterious pleiotropic effects on other structural or functional properties? The experimental results revealed that the sign and magnitude of mutational effects on Hb–O2 affinity were highly dependent on the genetic background in which they occurred, but affinity-altering mutations produced no significant pleiotropic effects on other aspects of Hb function.

Results and Discussion

For the experimental analysis of Hb function, we collected specimens of O. princeps from a high-altitude locality in the Southern Rockies and specimens of O. collaris from lowland localities in interior Alaska (see Materials and Methods). We experimentally confirmed that both species express a single tetrameric α2β2 Hb isoform during adult life, consistent with genomic evidence that O. princeps possesses single copies of postnatally expressed α- and β-globin genes (HBA and HBB, respectively; fig. 1A). We sequenced both HBA alleles and both HBB alleles from multiple individuals per species. The sequence data revealed that the Hbs of the two species are distinguished by five amino acid substitutions in the β-chain subunits (sites β5, β58, β62, β123, and β126; fig. 1B). Using blood samples from individuals with known HBA and HBB genotypes, we tested for differences in the oxygenation properties of purified Hbs from highland O. princeps (n = 4) and lowland O. collaris (n = 4) that were distinguished by the five β-chain substitutions. To gain insights into functional mechanisms that may be responsible for evolved differences in O2-binding properties, we measured O2-affinities of purified Hbs under four treatments: In the absence of allosteric effectors (stripped), in the presence of Cl− ions (added as 0.1 M KCl), in the presence of 2,3-diphosphoglycerate (DPG, at 2-fold molar excess over tetrameric Hb), and in the simultaneous presence of KCl and DPG. The (KCl+DPG) treatment most closely approximates in vivo conditions in mammalian red blood cells (Mairbäurl and Weber 2012), so from this point forward we focus on comparisons of Hb–O2 affinity in the presence of anionic effectors.

F

Pika globin genes. (A) Genomic structure of the α- and β-globin gene clusters of Ochotona princeps (GenBank accession numbers NT_167265 and AC234505, respectively). Isoelectric focusing experiments revealed that adult O. princeps and O. collaris both express a single, tetrameric α2β2 Hb isoform that incorporates the α-chain products of the HBA gene, and β-chain products of the HBB gene. Products of the remaining α- and β-like globin genes are not incorporated into functional Hb tetramers in definitive erythrocytes. Based on data from other eutherian mammals, products of three of the remaining globin genes are expected to encode α-type (HBZ [ζ-globin]) or β-type (HBE [ε-globin] and HBG [γ-globin]) subunits that are incorporated into structurally distinct Hb isoforms that are expressed at different stages of prenatal development. (B) Alignment of HBA and HBB amino acid sequences from seven pika species, including the two North American sister species, O. princeps and O. collaris. Orthologous sequences from the European rabbit (Oryctolagus cuniculus) are provided for comparison. In addition to the five β-chain substitutions at sites 5, 58, 62, 123, and 126 that distinguish the Hbs of O. princeps and O. collaris (denoted by boxes), there is one nearly fixed difference at site 78 in the α-chain.

The experiments revealed that the Hb of O. princeps exhibits a uniformly higher O2-affinity than that of O. collaris, as indicated by the lower values of P50 (the partial pressure of O2 at which Hb is 50% saturated; table 1 and fig. 2A). The Hb of O. princeps exhibited a significantly higher intrinsic O2 affinity than that of O. collaris, as indicated by the lower P50 for stripped Hb (P50(stripped)), and this difference persisted in the presence of Cl− ions (P50(KCl)), in the presence of DPG (P50(DPG)), and in the simultaneous presence of both allosteric effectors (P50(KCl+DPG); table 1 and fig. 2A). Measures of allosteric regulatory capacity (fig. 2B) revealed that the species differences in P50(KCl+DPG) were exclusively attributable to differences in intrinsic O2 affinity; there were no appreciable differences in the responsiveness to allosteric effectors. This was confirmed by dose–response curves for DPG, which demonstrated that Hb–O2 affinities of both species had identical association constants for DPG binding (fig. 3). Hbs of both species were also characterized by highly similar cooperativity coefficients (table 1).

Table 1.

O2-Affinities (P50, torr) and Cooperativity Coefficients (n50) (mean ± SEM) of Purified Hbs from High- and Low-Altitude Pika Species (Ochotona princeps and O. collaris, respectively).

	O. princeps (high altitude)	O. collaris (low altitude)
	(n = 4)	(n = 4)
P50 (torr)
Stripped	5.20 ± 0.29	6.23 ± 0.16
+KCl	10.07 ± 0.65	11.66 ± 1.32
+DPG	10.51 ± 0.63	12.91 ± 0.31
KCl + DPG	11.24 ± 0.69	13.24 ± 0.33
n50
Stripped	1.99 ± 0.21	2.35 ± 0.05
+KCl	2.53 ± 0.20	2.54 ± 0.06
+DPG	2.41 ± 0.09	2.56 ± 0.22
KCl + DPG	2.51 ± 0.16	2.58 ± 0.05

Note.—O2-equilibrium curves were measured at 0.3 mM heme in 0.1 M HEPES buffer at pH 7.4, 37 °C, in the absence (stripped) and presence of added KCl (0.1 M) and/or DPG (DPG/Hb tetramer ratio, 2.0), as indicated. Sample sizes refer to the number of individuals analyzed per species.

F

Hbs of the high-altitude Ochotona princeps exhibit uniformly higher O2-affinities than those of the low-altitude O. collaris. (A) P50 values (mean ± SEM) for purified Hbs of O. princeps (n = 4 individuals) and O. collaris (n = 4 individuals), measured at pH 7.4 and 37 °C in the absence (stripped) and presence of allosteric effectors ([Cl−], 0.1 M; [HEPES], 0.1 M; DPG/Hb tetramer ratio, 2.0; [heme], 0.3 mM). (B) Log-transformed differences in P50 values of Hbs from O. princeps and O. collaris in the presence and absence of allosteric effectors. The Δlog-P50 values measure the extent to which Hb–O2 affinity is reduced in the presence of a given allosteric effector (Cl−, DPG, or both anions together). See table 1 for raw data.

F

Variation in Hb–O2 affinity (as indexed by P50) as a function of DPG concentration for purified Hbs of Ochotona princeps and O. collaris. Curves were fit using a double hyperbolic function. The Hbs of the two pika species exhibit no discernible difference in DPG binding constants, with an apparent dissociation equilibrium constant for the high-affinity DPG site of approximately 0.02 mM.

Sequence data from additional species of pikas enabled us to reconstruct the β-chain sequence of the princeps/collaris common ancestor (Anc1; fig. 4A). For each of the five amino acid sites, posterior probabilities of estimated ancestral states were as follows: β5G (1.000), β58A (0.984), β62A (0.987), β123T (0.997), and β126A (0.947). The ancestral sequence reconstruction revealed that three of the five β-chain substitutions that distinguish the Hbs of the two species occurred on the branch leading to the high-altitude O. princeps, β5Gly→Ala (βG5A), β62Ala→Thr (βA62T), and β126Ala→Val (βA126V), and the remaining two substitutions occurred on the branch leading to the low-altitude O. collaris, β58Ala→Pro (βA58P) and β123Thr→Ser (βT123S) (fig. 4A). As a first step in the protein-engineering experiments, we synthesized, purified, and functionally tested rHbs representing the wildtype Hbs of O. princeps and O. collaris and the reconstructed Hb of their common ancestor, Anc1 (fig. 4A). A triangulated comparison involving the reconstructed Anc1 and its modern-day descendants revealed that the elevated Hb–O2 affinity of O. princeps is a derived character state: P50(KCl+DPG) for the highland O. princeps rHb was 21.8% lower (i.e., O2 affinity was higher) relative to Anc1 (supplementary table S1, Supplementary Material online). Consistent with our experiments on the native Hbs, P50(KCl+DPG) for the highland O. princeps rHb was 14.4% lower relative to the lowland O. collaris rHb. In conjunction with the ancestral sequence reconstruction, these results indicate that the derived increase in O2 affinity of O. princeps Hb was attributable to the independent or joint effects of the three substitutions at sites β5, β62, and β126 (fig. 4A).

F

Inferred history of substitution at five sites that distinguish the HBB genes of the high-altitude Ochotona princeps and the low-altitude O. collaris. (A) Of the five amino acid substitutions that distinguish the Hbs of the two sister species, ancestral sequence reconstruction revealed that three occurred on the branch leading to O. princeps, and the remaining two occurred on the branch leading to O. collaris. “Anc1” represents the reconstructed HBB sequence of the princeps/collaris common ancestor. The phylogeny of pika species is taken from (Lanier and Olson 2009). (B) Graphical representation of sequence space, where each vertex of the cube represents a discrete three-site β5–β62–β126 genotype and each edge connects genotypes separated by a single mutational step. The low-affinity ancestral Anc1 and the high-affinity Hb of O. princeps are connected by six (=3!) possible forward trajectories through sequence space.

To complement experimental measures of the functional effects of each amino acid substitution, we used a codon-based maximum-likelihood approach to test for evidence of positive selection in an alignment of HBA and HBB sequences from all examined pika species (see Materials and Methods). In the case of HBB, likelihood ratio tests (LRTs) revealed significant variation among sites in the ratio of nonsynonymous-to-synonymous substitution rates (ω = dN/dS), and a branch-sites test (Zhang et al. 2005; Yang and dos Reis 2011) implicated positive selection in driving substitutions at both β62 and β126 in O. princeps (supplementary table S2, Supplementary Material online).

We used site-directed mutagenesis to synthesize the inferred ancestral genotype (Anc1), the derived, triple-mutant genotype of O. princeps, and each of the six possible mutational intermediates that connect them (fig. 4B). This combinatorially complete set of three-site genotypes enabled us to measure the additive and nonadditive effects of princeps-specific substitutions on each of eight ( = 23) possible genetic backgrounds. Of the six (=3!) possible forward trajectories that lead from Anc1 to the derived, triple-mutant β-globin genotype of O. princeps, our experiments revealed that only a single pathway yielded an incremental increase in Hb–O2 affinity at each successive step: βG5A (first), βA62T (second), and βA126V (third; fig. 5). Each of the three possible first mutational steps produced a significant increase in Hb–O2 affinity, as did one of six possible second steps (A126V on the “AAA” background) and one of three possible final steps (βA126V on the “ATA” background; fig. 5). However, each mutation’s effect on Hb–O2 affinity was highly dependent on the genetic background in which it occurred, as one of six possible second steps (βG5A) and two of three possible final steps (βG5A and βA62T) yielded significant “reductions” in Hb–O2 affinity. The βG5A and βA62T substitutions exhibited sign epistasis for Hb–O2 affinity: Both substitutions significantly increased affinity on the Anc1 background, but they both produced opposite effects when they occurred on backgrounds in which the other substitutions had already occurred (fig. 5; supplementary table S1, Supplementary Material online). Similar sign-epistatic effects of affinity-altering mutations have been documented in protein engineering studies of other vertebrate Hbs (Natarajan et al. 2013; Projecto-Garcia et al. 2013).

F

Trajectories of change in Hb–O2 affinity (indexed by P50(KCl+DPG)) for each of six forward mutational pathways leading from the ancestral three-site β5–β62–β126 “Anc1” genotype (GAA) to the derived, triple-mutant genotype of O. princeps (ATV). Only one of six possible forward pathways yields a monotonic increase in Hb–O2 affinity: G5A (1st), A62T (2nd), A126V (3rd). Error bars denote 95% confidence intervals.

In contrast to the sign-epistatic effects of βG5A and βA62T, the βA126V substitution had an affinity-enhancing effect on each of the four possible genetic backgrounds (fig. 5), although the magnitude of the affinity-enhancing effect was highly context-dependent: βA126V produced a 32.2% reduction in P50(KCl+DPG) when it occurred as the first substitution (on the Anc1 background) compared with a 5.8% reduction when it occurred as the final substitution (on a background in which βG5A and βA62T had already occurred). Relative to the Anc1 reference genotype, pairwise epistasis was consistently negative for the free energy of Hb–O2 binding (ΔG), indicating that the combined effects of affinity-altering mutations were always less than the sum of their individual effects on the Anc1 background (fig. 6).

F

Relative to the ancestral genotype (Anc1), pairwise epistasis for the free energy of Hb–O2 binding, ΔG, was consistently negative (indicating negative epistasis for Hb–O2 affinity). ΔG values for the double- and triple mutant genotypes were measured relative to Anc1 and were plotted against values expected under an additive model. Observed ΔG’s for the full set of double- and triple-mutant genotypes fell below the y = x diagonal, reflecting the fact that the combined effects of affinity-altering mutations were consistently less than the sum of their individual effects on the Anc1 background. Error bars denote 95% confidence intervals.

Structural Basis of Mutational Effects

Homology modeling and comparative data from human Hb mutants provided insights into the structural mechanisms responsible for the individual functional effects of HBB substitutions in the O. princeps lineage. The βG5A substitution is predicted to stabilize the A-helix of the β-chain by reducing backbone entropy and by reducing the solvent accessibility of apolar surface (Lopez-Llano et al. 2006). This stabilization of α-helical conformation can propagate to the α1β2 and α2β1 intersubunit contacts that mediate oxygenation-linked transitions in quaternary structure, thereby perturbing the allosteric equilibrium between the high-affinity (R-state) and low-affinity (T-state) conformations (Dumoulin et al. 1997, 1998; Bellelli et al. 2006).

The βA62T substitution is predicted to increase the rigidity of the E-helix because the γ-oxygen atom of β62Thr forms a polar contact with the β58Ala main-chain oxygen atom (fig. 7). This localized stabilization of E-helix secondary structure produces a subtle conformational change in the β-chain distal heme pocket, thereby reducing steric hindrance of β heme ligand-binding. Similarly, a nonpolar→polar βA62P mutation at the same site in human Hb (Hb Duarte) also produces an increased O2 affinity while preserving normal cooperativity and DPG sensitivity (Beutler et al. 1974; Ceccarelli et al. 2006). The structural mechanisms responsible for the observed interaction effects are not immediately clear, as sites β5, β62, and β126 are structurally remote from one another in the folded β-chain subunit (fig. 7), indicating that epistatic interactions do not involve direct steric contact between residue side chains or short-range steric interaction mediated by a third residue or bound ligand. Instead, the observed epistatic effects must stem from indirect, second-order perturbations, possibly involving dynamics of the allosteric transition in quaternary structure between different oxygenation states of the Hb tetramer.

F

Homology model of pika Hb showing detailed view of the β1 subunit, and the locations of three amino acid substitutions (βG5A, βA62T, and βA126V) that occurred in the Ochotona princeps lineage (A). The proximal histidine, β92 (which covalently binds the fifth coordination site of the heme iron), and the distal histidine, β63 (which stabilizes bound ligand at the sixth coordination site on the opposite side of the heme plane), are shown in blue. Counting from the N-terminal valine of the β-chain polypeptide, β5 is the second residue of the A helix, β62 is the sixth residue of the E-helix, and β126 is the fourth residue of the H helix. Inset (top right) shows the full Hb tetramer, with the focal β1 subunit depicted as a ribbon structure; the β2 subunit is shown in blue, and the α1 and α2 subunits are shown in green. (B) The βA62T substitution increases the rigidity of the E-helix because the γ-oxygen atom of β62Thr forms a polar contact with the β58Ala main-chain oxygen atom. (C) The methyl side-chain of the ancestral β62Ala forms no contact with β58Ala, resulting in a less stable α-helical conformation.

Pleiotropic Trade-Offs

Analysis of the full set of rHb mutants revealed that Hb–O2 affinity in the presence of allosteric effectors (P50(KCl+DPG)) was not significantly correlated with allosteric regulatory capacity, as measured by both Cl− sensitivity (r = −0.263, P = 0.530) and DPG sensitivity (r = −0.200, P = 0.636), the cooperativity of O2-binding, as measured by Hill’s cooperativity coefficient, n50 (r = 0.313, P = 0.450), autoxidation rate (r = 0.066, P = 0.877), or structural mobility, as measured by a computationally predicted index of energetic frustration (Jenik et al. 2012) (r = −0.451, P = 0.262). Functional experiments on the engineered rHb mutants revealed sign and magnitude epistasis for Hb–O2 affinity, but affinity-altering mutations produced no appreciable effects on other structural or functional properties that potentially trade-off with O2-affinity.

The substitutions that we examined are substitutions that actually occurred in the ancestry of modern-day O. princeps. It may be that the pleiotropic and epistatic effects of substitutions that actually did occur are qualitatively or quantitatively distinct from those of all possible affinity-altering substitutions that could have occurred (Draghi and Plotkin 2013; Gong and Bloom 2014). More general conclusions about how pleiotropy influences the fixation probabilities of mutations can be obtained by comparing the spectrum of substitutions that have actually contributed to evolutionary changes in protein function with the spectrum of spontaneous mutations that affect protein function (Streisfeld and Rausher 2011).

Evolutionary Implications of Sign Epistasis between Affinity-Altering Mutations

Insights into the form and prevalence of intramolecular epistasis can illuminate the extent to which realized pathways of protein evolution represent deterministic or historically contingent outcomes. Epistasis between mutant sites in the same protein can exert a deterministic influence on molecular evolution by constraining the number of selectively accessible mutational pathways to high-fitness genotypes, thereby enhancing the repeatability (and hence, predictability) of adaptive walks through sequence space (Weinreich et al. 2005, 2006; DePristo et al. 2007; Lozovsky et al. 2009; da Silva et al. 2010; Rokyta et al. 2011). However, epistasis can also exert a stochastic influence because if the phenotypic effects of mutations are conditional on the genetic background in which they occur, then the range of selectively accessible mutational pathways will be historically contingent on which particular mutations happen to occur first (Kvitek and Sherlock 2011; Salverda et al. 2011; Harms and Thornton 2014; Kryazhimskiy et al. 2014).

If the derived increase in Hb–O2 affinity in the O. princeps lineage evolved under the influence of positive directional selection (as suggested by results of the branch-sites test; supplementary table S2, Supplementary Material online), then our results indicate that any one of the three observed β-chain substitutions could have occurred as the first step in the adaptive walk. However, as the initial substitution can significantly alter the sign and/or magnitude of the phenotypic effect of each subsequent mutation that occurs on the newly modified background, the selection coefficient (and, hence, fixation probability) of each new mutation will be historically contingent on which affinity-enhancing mutation happened to be the first one to fix. This role for contingency also holds in cases where the epistatically interacting sites are simultaneously polymorphic.

Conclusions

The functional effects of affinity-altering mutations in pika Hb were highly dependent on the sequential order in which they occurred. Two of the historical β-globin substitutions in O. princeps increased or decreased Hb–O2 affinity depending on whether they occurred as the first step or the last step in the mutational pathway. Regardless of the temporal order in which the three β-chain substitutions actually occurred during the evolution of the increased Hb–O2 affinity in O. princeps, our results revealed no evidence for ineluctable pleiotropic trade-offs that could constrain mutational pathways of Hb evolution. These results suggest that the accessibility of alternative mutational pathways may be more strongly constrained by sign epistasis for positively selected biochemical phenotypes than by antagonistic pleiotropy. The caveat is that this inference is only valid for the circumscribed region of sequence space that we examined, and for the particular biochemical properties that we investigated. In contrast to the deleterious pleiotropic effects documented for spontaneous affinity-altering mutations in human Hb (Bellelli et al. 2006; Storz and Moriyama 2008; Varnado et al. 2013), the mutations that have contributed to evolutionary changes in Hb–O2 affinity in O. princeps produced fine-tuned changes in O2-affinity with minimal pleiotropic effects. These findings are consistent with theory suggesting that function-altering mutations with minimal pleiotropic effects will tend to make disproportionate contributions to phenotypic evolution (Otto 2004; Streisfeld and Rausher 2011).

Materials and Methods

Specimen Collection

For the survey of HBA and HBB sequence variation and the functional analyses of native Hbs, we collected seven O. princeps near the summit of Mt. Evans, Clear Creek Co., Colorado (39°15′24″N, 106°10′54″W; 4,262 m above sea level [a.s.l.]) and a total of six O. collaris from Eagle Summit, Yukon-Koyukuk Co., Alaska (65°29′8″N, 145°24′0″W; 949 m a.s.l.) and Rainbow Ridge, Fairbanks North Star Co., Alaska (66°18′32″N, 145°40′21″W; 1,119 m a.s.l.). Specimens of O. princeps were collected under permit 10TR2058 from the Colorado Division of Wildlife, and were deposited as voucher specimens in the Harold W. Manter Laboratory collection of the University of Nebraska State Museum (HWML; catalog nos. NP1195, 1201–1202). Specimens of O. collaris were collected under permit 11-057 from the Alaska Department of Fish and Game and were deposited as voucher specimens in the University of Alaska Museum of the North (UAMN; UAM111535, 111549, 111591, 111675, 111698, and 111847). Pikas were euthanized in accordance with guidelines approved by the University of Nebraska Institutional Animal Care and Use Committee (IACUC no. 519). Blood was collected from each animal through cardiac puncture and red cell fractions were snap-frozen in liquid nitrogen. Liver tissue samples were also collected from each individual as a source of genomic DNA. Liver samples were immediately frozen in liquid nitrogen and stored at −80 °C prior to DNA extraction.

To conduct a broader survey of intraspecific DNA polymorphism in the HBA and HBB genes, we augmented our sample of wild-caught pikas with frozen tissue samples from additional specimens of O. princeps and O. collaris from the University of California Museum of Vertebrate Zoology (MVZ 202371) and UAMN (UAM 100795, 102429, 102435). To sequence orthologous HBA and HBB genes in other pika species, we obtained tissue samples from field-collected specimens of O. daurica (n = 3) from Khan Taishir Mountain, Mongolia (46°14′23″N, 96°21′30″E; 1,965 m a.s.l.) and O. pallasi (n = 1) from Jargalant Khairkhan Mountain, Mongolia (47°44′21″N, 92°27′24″E; 1,724 m a.s.l.). Voucher specimens were exported under permit G.12.15 and were deposited in the HWML collection (NK223524–223526, 223869). We obtained frozen tissue samples of O. rufescens (MVZ 191920) and O. hyperborea (UAM 84368), and we retrieved O. curzoniae HBA and HBB sequences from GenBank (EF429202 and DQ839484, respectively).

DNA Sequencing

For each pika specimen, genomic DNA was extracted from liver tissue using the DNeasy Blood and Tissue kit according to the manufacturer's protocol (Qiagen, Valencia, CA). We polymerase chain reaction (PCR)-amplified the HBA and HBB genes of each of the examined pika species using the Expand High Fidelity PCR system (Roche Applied Science, Indianapolis, IN). For the survey of polymorphism in each of the two focal species, we sequenced both alleles of the HBA and HBB genes (for each gene, n = 16 and 18 alleles for O. princeps and O. collaris, respectively). Thermal cycler conditions were as follows: An initial denaturing step for 2 min, followed by 35 cycles of 94 °C for 30 s, 55–57 °C for 30 s, 72 °C for 2 min, and final extension of 72 °C for 10 min. Primers for the PCR-amplification of the HBA and HBB genes are provided in supplementary table S3, Supplementary Material online. PCR products were purified and sequenced bidirectionally. All sequences were deposited in GenBank under the following accession numbers: KC757706–757708, KC757676–KC757683, KC757705, KC757674, KC757675, KC757684, KC757685, KC757646, KC757688, KC757689, KC757648–757653, KC757690, KC757654, KC757655, KC757691-757696, KC757656–757661, KC757697–757704, KC757662–757673, KC757686, KC757687, and KC757647.

Molecular Evolution Analyses

To make inferences about patterns of molecular evolution and to reconstruct ancestral sequences we estimated separate phylogenies for the HBA and HBB genes of pikas, using the rabbit orthologs as outgroup sequences. We aligned the nucleotide sequences based on their conceptually translated amino acid sequences using MAFFT (version 7.182, Katoh and Standley 2013) and we then performed phylogeny reconstructions using maximum likelihood as implemented in Treefinder, ver March 2011 (Jobb et al. 2004). We evaluated support for the nodes with 1,000 bootstrap pseudoreplicates. We used the “propose model” tool of Treefinder to select the best-fitting models of nucleotide substitution for each codon position based on the Akaike information criterion with correction for small sample size. The estimated trees were used for all downstream analyses.

To reconstruct the HBB amino acid sequence of the princeps/collaris common ancestor, we used a maximum-likelihood approach (Yang et al. 1995) as implemented in PAML 4.8 (Yang 2007). To test for evidence of positive selection in pika HBA and HBB genes, we quantified variation in ω (=dN/dS, the ratio of the rate of nonsynonymous substitution per nonsynonymous site [dN] to the rate of synonymous substitution per synonymous site [dS] site) using a maximum-likelihood approach (Goldman and Yang 1994) implemented in the CODEML program, v. 4.8 (Yang 2007). In all cases, we used LRTs to compare nested sets of models (Yang 1998). We first compared models that allow ω to vary among codons in the alignments of HBA and HBB coding sequence (M0 vs. M3, M1a vs. M2a, M7 vs. M8). The latter two LRTs are tests of positive selection, as the alternative models (M2a and M8) allow a subset of sites to have ω > 1, in contrast to the null models (M1a and M7) that only include site classes with ω < 1. These analyses revealed statistically significant variation in ω among sites in the HBB alignment (supplementary table S2, Supplementary Material online), so we then implemented branch and branch-site analyses. The branch-site analyses permit inferences about changes in the selective regime in a subset of codons between foreground and background branches (Yang and Nielsen 2002; Yang et al. 2005; Yang and dos Reis 2011). In the alternative hypothesis (Model A), a fraction of sites in foreground branches are permitted to shift to a regime of positive selection (indicated by ω > 1); under the null model, in contrast, sites are only permitted to shift between neutral evolution and purifying selection. We implemented three different branch models: 1) A single-ω model, where all branches in the tree had the same ω, which corresponds to M0 from the sites analyses; 2) a two-ω model, involving a joint estimate of ω for branches subtending the two high-altitude species (O. curzionae and O. princeps) and a separate estimate of ω for the remaining “background” branches; and 3) a three-ω model with separate ω estimates for O. curzionae, O. princeps, and all background branches. We then conducted branch-site analyses of HBB sequences that focused specifically on O. princeps. We compared two sets of models, one in which the branch subtending the set of O. princeps sequences was set as the foreground branch, and one in which the entire set of O. princeps sequences (O. princeps clade) was set as foreground. We used the Bayes Empirical Bayes approach to calculate the posterior probability that a given site belongs to the site class with ω > 1 (i.e., that the substitution at the site in question was driven by positive selection in the foreground branch) (Yang et al. 2005).

Native Hb Purification

We used isoelectric focusing (PhastSystem, GE Healthcare Bio-Sciences, Piscataway, NJ) to characterize the Hb isoform composition of circulating red blood cells in both O. princeps and O. collaris. Native pika Hbs were purified by means of anion-exchange chromatography (HiTrap QHP, GE Healthcare, Piscataway, NJ), as described previously (Storz et al. 2012; Revsbech et al. 2013). The column was pre-equilibrated with 20 mM Tris (pH 8.4) and the sample was eluted using a linear gradient of 0–0.2 M NaCl. Samples were desalted by dialysis against 10 mM HEPES buffer (pH 7.4) and were concentrated using Millipore centrifugal filter units (30-kDa cutoff). Hb concentrations were determined from absorbance spectra using standard extinction coefficients at 577 and 540 nm.

Vector Construction and Site-Directed Mutagenesis

Nucleotide sequences of pika HBA and HBB genes were optimized with respect to Escherichia coli codon preferences, and were then synthesized by Eurofins (Huntsville, AL). Gene cassettes for the HBA and HBB genes and the methionine aminopeptidase (MAP) gene were tandemly cloned into the custom expression vector (pGM plasmid) described by Natarajan et al. (2011, 2013). In order to maximize efficiency in the posttranslational cleaving of N-terminal methionine residues from the nascent α- and β-chain polypeptides, an additional copy of the MAP gene was cloned into a coexpressed plasmid, pCO-MAP, with a kanamycin resistance gene. Site-directed mutagenesis was performed using the QuikChange II XL Site-Directed Mutagenesis kit from Stratagene (La Jolla, CA). Each engineered codon change in the globin cassette was verified by plasmid sequencing. Sequences for the mutagenic primers are provided in supplementary table S3, Supplementary Material online.

Expression and Purification of rHb

All rHb mutants were expressed in the JM109 (DE3) E. coli strain as described previously (Natarajan et al. 2011, 2013; Projecto-Garcia et al. 2013; Cheviron et al. 2014). Bacterial cells were subject to dual selection in an LB agar plate containing ampicillin and kanamycin, which ensures that transformed cells receive both the pGM and pCO-MAP plasmids. High-level expression of MAP enzyme proved critical for synthesizing fully functional rHbs. All pika rHb mutants were expressed in identical experimental conditions. Large-scale production was conducted in batches containing 1.5–2.0 l of TB medium. Cells were grown at 37 °C in an orbital shaker at 200 rpm until the absorbance reached an optical density of 0.6 at 600 nm. Expression of the HBA, HBB, and MAP genes was induced with 0.2 mM IPTG and the culture was supplemented with hemin (50 µg/ml) and glucose (20 g/l). The cells were then grown at 28 °C for 16 h in an orbital shaker at 200 rpm. The bacterial culture was saturated with CO for 15 min and the cells were harvested by centrifugation and stored at −80 °C.

As a preparative step for the protein purification, the cells were resuspended in lysis buffer (50 mM Tris, 1 mM ethylenediaminetetraacetic acid [EDTA], 0.5 mM DTT, pH 7.6) with lysozyme (1 mg/g wet cells) and were incubated in an ice bath for 30 min. The bacterial cells were sonicated, followed by addition of polyethyleneimine solution (0.5–1%) and centrifugation at 13,000 rpm (4 °C) for 45 min. The clarified supernatant was then subjected to overnight dialysis in Tris buffer (20 mM Tris–HCl with 0.5 mM EDTA, pH 8.3). The rHbs were then purified by two-step ion-exchange chromatography. We used prepacked Q-Sepharose columns (HiTrap QHP, 5 ml, 17-1154-01; GE Healthcare) equilibrated with Tris-buffer (20 mM Tris–HCl with 0.5 mM EDTA, pH 8.3). The samples were passed through the column and the rHb solutions were eluted against a linear gradient of 0–0.5 M NaCl. The eluted samples were desalted by overnight dialysis with phosphate buffer (10 mM sodium phosphate buffer, 0.5 mM EDTA, pH 6.8). Each dialyzed rHb sample was then passed through an SP-Sepharose column (HiTrap SPHP, 1 ml, 17-1151-01; GE Healthcare) with an elution buffer (20 mM sodium phosphate buffer, 0.5 mM EDTA, pH 8.3). Samples were concentrated and desalted by overnight dialysis against 10 mM HEPES buffer (pH 7.4) and were stored at −80 °C prior to the measurement of O2-equilibrium curves and autoxidation rates. The purity of each rHb preparation was confirmed by 20% sodium dodecyl sulfate-polyacrylamide gel electrophoresis and measures of the absorbance spectra of oxy, deoxy, and CO derivatives at 450–600 nm.

Experimental Measurements of Hb Function

Measures of Hb–O2 affinity were conducted using a modified diffusion chamber where changes in absorbance were recorded in conjunction with stepwise changes in the O2 tension of equilibration gases, as described previously (Weber 1992; Runck et al. 2010; Grispo et al. 2012; Storz et al. 2012; Revsbech et al. 2013). Measured O2-equilibrium curves were based on 4–6 saturation steps. Values of P50 (partial pressure of O2 [PO2] at half-saturation) and the cooperativity coefficient, n50, were calculated by using nonlinear regression to fit the O2-saturation data to the Hill equation (Y = PO2/[P50 + PO2], where Y is the fractional O2 saturation and n is the cooperativity coefficient). O2-binding equilibria of native and rHb solutions (0.3 mM heme) were measured at 37 °C in 0.1 M HEPES buffer, pH 7.4, in the absence (stripped) and presence of added allosteric effectors (0.1 M KCl, and/or DPG at 2.0-fold molar excess over tetrameric Hb). These are standard experimental conditions (Imai 1982; Mairbäurl and Weber 2012) that closely approximate intraerythrocytic effector concentrations in vivo: 100 mM Cl− approximates the naturally occurring concentration inside mammalian red blood cells (where K+ is the main counter ion, as opposed to Na+ in the plasma) and a DPG/Hb tetramer ratio = 2.0 falls within the physiological range for most eutherian mammals (Bunn 1971; Duhm 1971; Mairbäurl and Weber 2012; Tufts et al. 2013). Hb–O2 equilibria at physiological concentrations of Cl− ions and DPG ensure that the in vitro measurements are relevant to in vivo conditions (Mairbäurl and Weber 2012). Possible discrepancies in the measured P50 values for native Hbs and rHb mutants could stem from the fact that the native Hb samples are characterized by some degree of structural heterogeneity due to heterozygosity at α- and/or β-globin sites, whereas rHbs have an invariant amino acid composition. For this reason, all inferences are based on comparisons between native Hbs of the two species (table 1), or comparisons among the rHb mutants (supplementary table S1, Supplementary Material online).

To characterize the pleiotropic effects of affinity-altering mutations, we measured four properties that can potentially trade-off with Hb–O2 affinity: Structural mobility, cooperativity of O2-binding, allosteric regulatory capacity (the sensitivity of Hb–O2 affinity to the inhibitory effects of anionic effectors), and susceptibility to spontaneous heme oxidation (autoxidation rate; supplementary table S1, Supplementary Material online). Although there are numerous potential pleiotropic effects of amino acid mutations that could be measured, we chose to focus on four quantifiable properties that are of known importance and physiological relevance (Olson and Maillett 2005; Bonaventura et al. 2013; Varnado et al. 2013). To predict structural mobility of tetrameric Hb, we used the crystal structure of European hare Hb (Lepus europaeus; PDB ID, 3LQD) to build a structural model of pika Hb using SWISS-MODEL (Arnold et al. 2006); derivatives were modeled using MODELLER 9V11 (Sali and Blundell 1993). Molecular interactions were identified by PISA (Krissinel and Henrick 2007), and molecular graphics were generated using PyMol (Schrödinger, LLC). Conformational stress and ligand binding energy were computed using the Frustratometer (Jenik et al. 2012) and AutoDock (Morris et al. 2009) programs. As a measure of cooperativity, we used Hill’s coefficient (n50) in the presence of Cl− and DPG. As a measure of allosteric regulatory capacity, we measured the difference in log-transformed P50 values for stripped Hb in the presence and absence of KCl and DPG. Following standard protocols for vertebrate Hbs (Jensen 2001), we measured autoxidation rates of purified rHbs at 10 µM concentration in 50 mM Tris buffer pH 7.4 at 37 °C using an Agilent 8453 UV-Visible Spectrophotometer (Agilent Technologies, Santa Clara, CA). Absorbance was recorded every hour at 560, 576, and 630 nm over a 14-h period. For each experiment, measurements of autoxidation rates were derived from the average of absorbance traces from the three wavelengths.

Supplementary Material

Supplementary tables S1–S3 are available at Molecular Biology and Evolution online (http://www.mbe.oxfordjournals.org/).