Author Response: Functional Trade-Offs and Environmental Variation Shaped Ancient Trajectories in the Evolution of Dim-Light Vision

Article Figures and data Abstract eLife digest Introduction Results Discussion Materials and methods References Decision letter Author response Article and author information Metrics Abstract Trade-offs between protein stability and activity can restrict access to evolutionary trajectories, but widespread epistasis may facilitate indirect routes to adaptation. This may be enhanced by natural environmental variation, but in multicellular organisms this process is poorly understood. We investigated a paradoxical trajectory taken during the evolution of tetrapod dim-light vision, where in the rod visual pigment rhodopsin, E122 was fixed 350 million years ago, a residue associated with increased active-state (MII) stability but greatly diminished rod photosensitivity. Here, we demonstrate that high MII stability could have likely evolved without E122, but instead, selection appears to have entrenched E122 in tetrapods via epistatic interactions with nearby coevolving sites. In fishes by contrast, selection may have exploited these epistatic effects to explore alternative trajectories, but via indirect routes with low MII stability. Our results suggest that within tetrapods, E122 and high MII stability cannot be sacrificed—not even for improvements to rod photosensitivity. https://doi.org/10.7554/eLife.35957.001 eLife digest People can see in dim light because of cells at the back of the eye known as rods. These cells contain two key components: molecules called retinal, which are bound to proteins called rhodopsin. When light hits a rod cell, it kicks off a cascade of reactions beginning with the retinal molecule changing into an activated shape and ending with a nerve impulse travelling to the brain. The activated form of retinal is toxic, and as long as it remains bound to the rhodopsin protein it will not damage the rod or surrounding cells. The toxic retinal also cannot respond to light. It must be released from the protein and converted back to its original shape to restore dim light vision. As with all proteins, rhodopsin’s structure comprises a chain of building blocks called amino acids. Every land animal with a backbone has the same amino acid at position 122 in its rhodopsin. This amino acid, named E122, helps to stabilize the activated rhodopsin, slowing the release of the toxic retinal. Yet E122 also makes the rod cells less sensitive, resulting in poorer vision in dim light. In contrast, some fish do not have E122 but rather one of several different amino acids takes its place. What remains unclear is why all land animals have stuck with E122, and whether there were other options that evolution could have explored to overcome the trade-off between sensitivity and stability. By looking at the make-up of rhodopsins from many animals, Castiglione and Chang found other sites in the protein where the amino acid changed whenever position 122 changed. The amino acids at these so-called “coevolving sites” were then swapped into the version of rhodopsin that is found in cows, which had also been engineered to lack E122. These changes fully compensated for the destabilizing loss of E122 on activated rhodopsin but without sacrificing its sensitivity to light. Further experiments then confirmed that unless all amino acids were substituted at once, the activated rhodopsin was very unstable. Indeed, it was almost as unstable as mutated rhodopsins found in some human diseases. These findings suggest that, while there was in principle another solution available to land animals, the routes to it were closed off because they all came with an increased risk of eye disease. These findings highlight that rhodopsin likely plays a more important role in protecting humans and many other land animals against eye disease than previously assumed. More knowledge about this protective role may lead to new therapies for these conditions. Also, investigating similar evolutionary trade-offs could help to explain how and why different proteins work the way that they do today. https://doi.org/10.7554/eLife.35957.002 Introduction Nature-inspired strategies are increasingly recruited toward engineering objectives in protein design (Khersonsky and Fleishman, 2016; Jacobs et al., 2016; Goldenzweig and Fleishman, 2018, a central challenge of which is to successfully manipulate backbone structure to modulate stability without introducing undesirable pleiotropic effects on protein activity (Khersonsky and Fleishman, 2016; Goldenzweig and Fleishman, 2018; Starr and Thornton, 2017; Tokuriki and Tawfik, 2009). Engineering protein stability and activity requires an understanding of a protein’s sequence-function relationship, or landscape (Pál and Papp, 2017; Wu et al., 2016; Starr et al., 2017), where billions of possible pair-wise and third-order interactions can exist between amino acids (Starr and Thornton, 2017; Storz, 2016), and only a limited number of amino acid combinations will confer the function of interest (Wu et al., 2016; Starr et al., 2017; McMurrough et al., 2014; Mateu and Fersht, 1999; Tarvin et al., 2017). To understand the context-dependence of amino acid functional effects (also known as intramolecular epistasis [Starr et al., 2017; Storz, 2016; Echave et al., 2016]), approaches such as deep mutational scanning (Wu et al., 2016; Starr et al., 2017; Sailer and Harms, 2017) can explore a subset of sequence-function space formed in response to a limited set of artificial selection pressures (Starr and Thornton, 2017). By contrast, natural protein sequence variation reflects the range of protein function that evolved in response to changing ecological variables (Starr and Thornton, 2017; Pál and Papp, 2017; Ogbunugafor et al., 2016), where convergent ‘solutions’ for protein function and stability can be derived through the evolution of alternative protein sequences (McMurrough et al., 2014; Mateu and Fersht, 1999; Tarvin et al., 2017). This suggests that closer examination of natural sequence variation may reveal new blueprints for protein design. The dim-light visual pigment rhodopsin (RH1/RHO) is an excellent model for understanding how both ecological variables and biophysical pleiotropy may interact to determine the availability of functional evolutionary solutions for environmental challenges (Kojima et al., 2017; Gozem et al., 2012; Dungan and Chang, 2017; Castiglione et al., 2018). Spectral tuning mutations that shift the RH1 wavelength of maximum absorbance (λMAX) can adapt dim-light vision to a remarkable range of spectral conditions across aquatic and terrestrial visual ecologies (Hunt et al., 2001; Hauser and Chang, 2017a; Dungan et al., 2016). Recently, λMAX was revealed to exist within a complex series of epistasis-mediated trade-offs with the non-spectral functional properties of RH1 long understood as adaptations for dim-light (Gozem et al., 2012; Dungan and Chang, 2017; Castiglione et al., 2017; Hauser et al., 2017b). These include an elevated barrier to spontaneous thermal-activation, which minimizes rod dark noise and is promoted by blue-shifts in λMAX (Kojima et al., 2017; Gozem et al., 2012; Kefalov et al., 2003; Yue et al., 2017); and a slow decay of its light-activated conformation, which we refer to here as metarhodopsin-II (MII) for simplicity (Imai et al., 1997; Lamb et al., 2016; Kojima et al., 2014; Sommer et al., 2012; Schafer et al., 2016; Van Eps et al., 2017). The RH1 MII active conformation is associated with rapid and efficient activation of G-protein transducin (Gt) (Kojima et al., 2014; Sugawara et al., 2010), yet the reasons for its long-decay after Gt signaling remain unclear (Kefalov et al., 2003; Imai et al., 1997; Imai et al., 2007). To sustain vision, all-trans retinal (atRAL) chromophore must be released from MII after Gt signaling (Palczewski, 2006)—a process that depends on the conformational stability of the MII-active-state structure (Schafer et al., 2016; Schafer and Farrens, 2015). Cone opsins have low MII stability and therefore rapidly release atRAL (Imai et al., 1997; Chen et al., 2012a), where it is quickly recycled back into 11-cis retinal (11CR) through the cone visual (retinoid) cycle, enabling rapid regeneration of cone pigments for bright-light vision (Wang and Kefalov, 2011; Tsybovsky and Palczewski, 2015). Rods, in contrast, regenerate thousands of times slower than cones after bright-light exposure (Mata et al., 2002). Indeed, rod exposure to bright flashes of light leads to atRAL release that can outpace clearance by visual cycle enzymes (Sommer et al., 2014; Rózanowska and Sarna, 2005), thus leading to accumulation (Saari et al., 1998; Lee et al., 2010) and light-induced retinopathy through various modes of cellular toxicity involving oxidative stress (Maeda et al., 2009; Chen et al., 2012b). Interestingly, recent biochemical evidence suggests MII may play a role in retinal photoprotection by complexing with arrestin after Gt signaling to re-uptake and thus provide a sink for toxic atRAL after rod photobleaching (Sommer et al., 2014). This suggests the evolution of rhodopsin’s high conformational selectivity for toxic atRAL may be a functional specialization (Schafer et al., 2016; Schafer and Farrens, 2015), which could in turn reflect differences in retinoid metabolism between rods vs. cones (Wang and Kefalov, 2011; Tsybovsky and Palczewski, 2015; Imai et al., 2005). Consistent with the overlapping mechanisms of RH1 spectral and non-spectral functions via the highly constrained RH1 structure (Gozem et al., 2012; Yue et al., 2017), this biophysical pleiotropy likely necessitates costly trade-offs between the spectral and non-spectral functions of RH1 in natural systems (Dungan and Chang, 2017; Luk et al., 2016). By comparison, directed evolution and synthetic biology approaches have successfully engineered either spectral, or non-spectral aspects of rhodopsin function, but did not address trade-offs arising from shifts in function. It has thus been possible to shift the spectral absorbance of archaea and bacterial rhodopsins close to the limit of the visible spectrum (Herwig et al., 2017; McIsaac et al., 2014), and to engineer tetrapod rhodopsins with high thermal stability (Xie et al., 2003), constitutive activation (Deupi et al., 2012; Standfuss et al., 2011), and alternative chromophore-binding sites (Devine et al., 2013). However, it has not been investigated whether rod visual pigments with novel combinations of spectral and non-spectral functional properties can be engineered by manipulating the biophysical pleiotropy of RH1 otherwise exploited by natural selection. Site 122 (Bos taurus RH1 numbering) is a molecular determinant of both the spectral and non-spectral functional properties of rhodopsin and the cone opsins (Hunt et al., 2001; Yue et al., 2017; Imai et al., 1997; Imai et al., 2007; Yokoyama et al., 1999). Intriguingly, vertebrate visual pigment families show differences in which amino acid variants predominate at this site (Figure 1A), with I122 strongly conserved in the most ancestrally diverging cone opsins such as the long-wave sensitive opsins (LWS) (Lamb et al., 2007), whereas in the most derived opsin group, the rhodopsins (RH1), E122 predominates (Figure 1B,C) (Imai et al., 1997; Lamb et al., 2007; Imai et al., 2007; Carleton et al., 2005). E122 is a key component of an important hydrogen-bonding network with H211 that is known to stabilize the MII active-conformation (Choe et al., 2011). This stability increase is so dramatic that E122 is considered a functional determinant distinguishing rhodopsin from cone opsins (Figure 1B) (Imai et al., 1997; Lamb et al., 2016; Kojima et al., 2014). Paradoxically, by conferring this increase in MII stability, the evolution of E122 likely involved a costly fitness trade-off that diminished tetrapod rod photosensitivity (Yue et al., 2017), which can affect visual performance in animals (Kojima et al., 2017; Aho et al., 1988). Indeed, it is possible to improve tetrapod rod photoreceptor sensitivity by decreasing rod dark noise in vivo by replacing E122 with a cone opsin amino acid variant (COV; Figure 1A) at site 122, such as Q122, which predominates in RH2 cone opsins (Yue et al., 2017; Lin et al., 2017). The strict conservation of E122 in all tetrapod rhodopsins (Figure 1C, Table 1, Supplementary file 1) therefore suggests that during the evolution of tetrapod dim-light vision, natural selection may have prioritized MII stability (Figure 1B,D) at the expense of rod sensitivity. This apparent evolutionary trade-off is perplexing given that the low spontaneous thermal activation of rhodopsin (and therefore rod dark noise) is a functional hallmark of rhodopsin divergence from the cone opsins (Kojima et al., 2017; Gozem et al., 2012; Kefalov et al., 2003; Lamb et al., 2016). Figure 1 with 2 supplements see all Download asset Open asset Natural variation at site 122 determines rhodopsin function and stability. (A) Amino acid consensus residues at site 122 across vertebrate rod opsins (rhodopsin; RH1) and the cone opsins (long-wave (LWS), short-wave (SWS1 and SWS2) and middle-wave (RH2) sensitive). Modified from (Lamb et al., 2007). (B) Relative stability of the rod and cone opsin active-conformation (MII) in different vertebrates (Imai et al., 2005). (C) Schematic representation of naturally occurring cone opsin variants (COVs) and other amino acids across vertebrate RH1 (see Figure 1—figure supplements 1–2; Tables 1–2, Supplementary files 1–2). E122 is invariant in all Tetrapod RH1 genes sequenced to date. Natural deep-sea amino acid variants (Hunt et al., 2001; Yokoyama et al., 1999) are identified with an asterisk (*; Table 2). (D) Introduction of the ancestral cone opsin (LWS) variant I122 blue shifts tetrapod RH1 spectral absorbance and accelerates decay of the MII light-activated conformation. https://doi.org/10.7554/eLife.35957.003 Table 1 Variation at sites 119-122-123-124 in Tetrapods and Outgroup rh1. Sites with variation relative to the Vertebrate consensus (LEIA) are in bold and highlighted grey. Subterranean species are denoted (*). https://doi.org/10.7554/eLife.35957.006 SpeciesAccessionCommon name119122123124OutgroupsCallorhinchus miliiXP_007888679Elephant sharkLEIGOrectolobus ornatusAFS63882Ornate wobbegongLEVSLatimeria chalumnaeXP_005997879CoelacanthLQVANeoceratodus forsteriABS89278Australian lungfishFIIAMammalsDasypus novemcinctusXP_0044773039-banded armadillo*IEIAEptesicus fuscusXP_008150514Big brown batLEVAChrysochloris asiaticaXP_006868732Cape golden mole*MEIASorex araneusXP_004613289Common Shrew*LEVATupaia chinensisXP_006160726Tree ShrewLEVAIctidomys tridecemlineatusXP_00533384113-line ground squirrelLEVARattus norvegicusNP_254276brown ratLEIGSarcophilus harrisiiXP_003762497Tasmanian devilTEVAReptilesAlligator mississippiensisXM_006274155American alligatorLEVAAlligator sinensisXP_006039462Chinese alligatorLEVAAnolis carolinensisNP_001278316Carolina anoleLEMGPython bivittatusXP_007423324Burmese pythonLEMAAmphibiansAmbystoma tigrinumU36574Tiger salamanderMEIACynops pyrrhogasterBAB55452Jap. Fire belly newtLEIGXenopus tropicalisNP_001090803Western clawed frogLEMAXenopus laevisNP_001080517African clawed frogLEVA Why has tetrapod RH1 been constrained to this paradoxical compromise at site 122 for the last 350 million years? Interestingly, and in contrast to tetrapod rhodopsins, fish rhodopsins show variation at site 122, such as in the Coelacanth (Latimeria chalumnae), Lungfish (Neoceratodus forsteri), and deep-sea fish lineages, where COV (I, Q, M) and other residues at site 122 (V, D) are found (Figure 1C; Figure 1—figure supplement 1; Tables 1–2, Supplementary file 2) (Hunt et al., 2001; Yokoyama et al., 1999; Carleton et al., 2005). These substitutions have been shown to blue-shift λMAX by up to ~10 nm (Hunt et al., 2001; Yokoyama et al., 1999), and may improve dim-light sensitivity in poorly-lit aquatic environments (Yue et al., 2017). Strikingly, one of the largest freshwater groups—the Characiphysi (which includes piranhas, electric eels, and catfishes [Chen et al., 2013]) —has the COV I122 residue completely fixed (Figure 1—figure supplement 1, Supplementary file 3). In tetrapods by contrast, the red-shifting E122 mutation is strictly maintained, increasing MII stability (Imai et al., 1997) but greatly decreasing rod sensitivity (Yue et al., 2017). Why the strong constraints on high MII stability and E122 are relaxed only within certain aquatic visual ecologies, remains unknown. Table 2 Fish rh1 with variation at site 122 do not necessarily have variation at coevolving sites 119, 123, and 124. Sites with variation relative to the Vertebrate consensus (LEIA) are in bold and highlighted grey. https://doi.org/10.7554/eLife.35957.007 OrderSpeciesAccessionCommon name119122123124Ecology notes from FishBaseLepisosteiformesLepisosteus oculatusJN230969.1spotted garLMISFreshwater; brackish; demersal. (Ref. 2060)Atractosteus tropicusJN230970.1Tropical GarLMLSFreshwater; demersalOsteoglossiformesMormyrops anguilloidesJN230973.1Cornish JackTIIAFreshwater; demersal; potamodromous (Ref. 51243)Osteoglossum bicirrhosumKY026030.1Silver arowanaTIIAFreshwater; benthopelagic AlepocephalifromesAlepocephalus bicolorJN230974.1Bicolor slickheadLQIAMarine; bathydemersal; depth range 439–1080 m (Ref. 44023).Bathytroctes microlepisJN544540.1Smallscale smooth-headLDIAMarine; bathypelagic; depth range 0–4900 m (Ref. 58018)Conocara salmoneumJN412577.1Salmon smooth-headLQIAMarine; bathypelagic; depth range 2400–4500 m (Ref. 40643) GalaxiiformesGalaxias maculatusJN231000.1InangaLMIGMarine; freshwater; brackish; benthopelagic; catadromous (Ref. 51243).StomiatiformesArgyropelecus aculeatusJN412571.1Lovely HatchetfishHQIAMarine; bathypelagic; depth range 100–2056 m (Ref. 27311)Vinciguerria nimbariaJN412570.1Oceanic lightfishHQVAMarine; bathypelagic; depth range 20–5000 m (Ref. 4470) AteleopodiformesAteleopus japonicusKC442218.1Pacific Jellynose FishLMISMarine; bathydemersal; depth range 140–600 m (Ref. 44036). MyctophiformesBenthosema suborbitaleJN412576.1Smallfin lanternfishHQVGMarine; bathypelagic; oceanodromous; depth range 50–2500 m (Ref. 26165)Lampanyctus alatusJN412575.1Winged lanternfishHQVAMarine; bathypelagic; oceanodromous; depth range 40–1500 m (Ref. 26165)Neoscopelus microchirKC442224.1Shortfin neoscopelidLQIAMarine; bathypelagic; depth range 250–700 m (Ref. 4481)GadiiformesCoryphaenoides guentheriJN412578.1Gunther’s grenadierLVIAMarine; bathydemersal; depth range 831–2830 (Ref. 1371) BeryciformesMelamphaes suborbitalisJN231006.1Shoulderspine bigscaleLQIAMarine; brackish; bathypelagic; depth range 500–1000 m (Ref. 31511).HolocentriformesHolocentrus rufusKC442230.1Longspine squirrelfishLMISMarine; reef-associated; depth range 0–32 m (Ref. 3724).Myripristis murdjanKC442231.1Pinecone soldierfishLMIGMarine; reef-associated; depth range 1–50 m (Ref. 9710) ScombriformesAphanopus carboEU637938.1Black scabbardfishHQIGMarine; bathypelagic; oceanodromous (Ref 108735); 200–2300 m (Ref. 108733)Cubiceps gracilisEU637952.1Driftfish-QIAMarine; pelagic-oceanic; oceanodromous (Ref. 51243); In light of these ecological patterns, we questioned whether it was possible to synthesize an evolutionary alternative: a tetrapod RH1 that never lost COV at site 122 but still developed high MII stability. We reasoned that relative to tetrapods, the diversity and complexity of fish visual ecologies (Hunt et al., 2001; Hauser and Chang, 2017a) may have allowed selection the opportunity to explore the pleiotropic potential of site 122 through the evolution of novel structural interactions with nearby sites that could compensate for the destabilizing loss of the E122-H211 hydrogen bond. To identify these interactions, our goal was to use analyses of evolutionary rates to predict sites coevolving with site 122, and to investigate the functional consequences of coevolving sites with experimental site-directed mutagenesis studies. Ultimately, we used our analyses of natural variation as a guide to artificially engineer a tetrapod rhodopsin with increased MII stability, but within a non-E122 sequence background. We demonstrated that this synthetic alternative is possible, even if evolution did not proceed down this mechanistic trajectory toward a dim-light adapted visual pigment. Results Phylogenetic identification of an intramolecular coevolutionary network To better understand the selection pressures that may be constraining E122 to fixation during tetrapod evolution, we constructed a large vertebrate rhodopsin phylogenetic dataset (Figure 1—figure supplements 1 and 2, Supplementary files 1–2) and investigated the evolutionary history of site 122 using ancestral reconstruction (Materials and methods). We found that E122 (codon GAA; Figure 2A) has been fixed in tetrapod RH1 since the most recent common ancestor ~350 million years ago (MYA) (Hedges et al., 2015), where it appears along the ancestral branch leading to tetrapods (Figure 2A; Table 3) following the diversification from lungfishes (I122, codon ATA, Figure 2A; Supplementary file 1) and the coelacanth (Q122, codon CAA, Figure 2A; Supplementary file 1). This transition period in vertebrate evolution is characterized by extensive morphological modifications for vision within terrestrial environments, and likely included large increases in environmental light irradiance (MacIver et al., 2017; Warrant and Johnsen, 2013). Apart from the lungfishes and coelacanth, the high conservation of E122 in tetrapods is also reflected in other vertebrate rhodopsins (Figure 1, Figure 1—figure supplement 1; Tables 1–2; Supplementary files 1–2), but there are important exceptions within certain lineages of teleost fishes, such as the Characiphysi. Within this group, the COV residue I122 was introduced likely through E122I (codon ATC; Figure 2A), where I122 is now completely fixed across the extant Characiphysi (Supplementary file 3). Figure 2 with 1 supplement see all Download asset Open asset Local coevolutionary forces govern the evolution of site 122 differentially between tetrapods and fish (teleost) RH1. (A) Extant and reconstructed codon variation at site 122 (Materials and methods). Despite a variety of residues at site 122 across the Coelacanth (Q122), Lungfish (I122; Ceratodontiformes), and Tetrapods (E122), GAA codons encoding for E122 are nevertheless predicted as the ancestral state with high posterior probabilities (shown in parentheses). E122 (GAA/GAG) is also likely to have been present in the last common ancestor of Cypriniformes and the Characiphysi, although with low posterior probabilities and therefore high uncertainty. I122 codon ATC is fixed in all Characiphysi rhodopsin to our knowledge (Supplementary file 3). Approximate divergence times are from (Hedges et al., 2015). (B) Mutual information (MI) analyses (MISTIC [Simonetti et al., 2013]) reveal all sites coevolving with site 122 are within 6 Å. Significance thresholds were determined by reference to the highest MI z-score from all sites across analyses of randomized datasets (n = 150; z-score cut-off = 21.6), as previously described (Ashenberg and Laub, 2013). (C) Sites within this radius displayed decreased amino acid variation in tetrapod and characiphysi RH1, where E122 and I122 are fixed, respectively (asterisks). (D) In tetrapods and characiphysi RH1, reduction in amino acid variation (relative to teleosts) at positions within the 6 Å radius were driven by increases in purifying selection on non-synonymous codons. Statistically significant gene-wide increases in purifying selection (*) between lineages were detected by likelihood ratio tests of alternative (Clade model C [Bielawski and Yang, 2004]) and null (M2a_REL [Weadick and Chang, 2012]) model analyses of codon substitution rates (dN/dS) ((p<0.001); Tables 3–5). Sites estimated to be under this increase in purifying selection (*) were those identified in the divergent site class of the CmC model analyses through a Bayes empirical Bayes analysis as previously described (Castiglione et al., 2017. Site-specific dN/dS estimates are from M8 analyses on phylogenetically pruned datasets (Tables 8–10; Figure 1—figure supplement 2; Figure 2—figure supplement 1). https://doi.org/10.7554/eLife.35957.008 Table 3 Results of Clade Model C (CmC) analyses of vertebrate rh1 under various partitions. https://doi.org/10.7554/eLife.35957.010 Model and Foreground†ΔAIC‡lnLParametersNullP [df]ω0ω1ω2/ωd M2a_rel225.5−47185.370.02 (69%)1 (3%)0.20 (28%)N/A-CmC_Tetrapod Branch97.44−47119.330.20 (28%)1 (3%)0.02 (69%) Tetra Br: 0.00M2a_rel0.000 [1]CmC_Tetrapod4.92−47073.060.02 (67%)1 (3%)0.24 (30%) Tetra: 0.13M2a_rel0.000 [1]CmC_Teleost1.88−47071.540.02 (67%)1 (3%)0.14 (30%) Teleost: 0.24M2a_rel0.000 [1]CmC_Teleost vs Tetrapod0*−47069.600.02 (67%)1 (3%)0.17 (30%) Tetra: 0.13 Teleost: 0.24M2a_rel0.000 [2] †The foreground partition is listed after the underscore for the clade models and consists of either: the clade of Teleost fishes (Teleost); the clade Tetrapods (Tetrapod;Tetra) or branch leading to tetrapods (Tetrapod branch; Tetra Br); or the clades of both the teleost fishes and tetrapods as two separate foregrounds (Teleost vs Tetrapods). In any partitioning scheme, the entire clade was tested, and all non-foreground data are present in the background partition. ‡All ΔAIC values are calculated from the lowest AIC model. The best fit is shown with an asterisk (*). ωd is the divergent site class, which has a separate value for the foreground and background partitions. Significant p-values (α ≤0.05) are bolded. Degrees of freedom are given in square brackets after the p-values. Abbreviations—lnL, ln Likelihood; p, p-value; AIC, Akaike information criterion. Since fishes (Teleosts), unlike tetrapods, display amino acid variation at site 122 (Figure 1C), we hypothesized that compensatory mutations may be coevolving with site 122 across fish RH1. To test this hypothesis, we investigated across the entire transmembrane domain of rhodopsin (residues 53–302) for evidence of sites coevolving with site 122 within an alignment of Teleost RH1 (Materials and methods; Supplementary file 2). Using phylogenetically corrected mutual information (MI) analyses (MISTIC; [(Simonetti et al., 2013]) with z-score cut-off determined by analyses of randomized datasets (Ashenberg and Laub, 2013), we found significant evidence of coevolution with site 122 at several RH1 positions, all of which clustered within 6 Å of E122 (Figure 2B) in the MII crystal structure (Choe et al., 2011. This is within the range at which intramolecular forces such as Van der Waals and hydrophobic interactions between amino acids are thought to occur (Ivankov et al., 2014). It is known, however, that there is a tendency of covariation analyses such as MI to identify coevolving sites proximal to each other, which may in turn overlook more distal coevolving sites potentially indirectly interacting with site 122 (Talavera et al., 2015). Nevertheless, sites detected within this 6 Å radius (sites 119, 123) have been previously found capable of functionally compensating for human pathogenic mutations (e.g. A164V) disrupting the MII-stabilizing E122-H211 interaction (Stojanovic et al., 2003), suggesting that natural variation at coevolving sites within this radius could compensate for the functional effects of COV at site 122. We therefore decided to focus our investigations on identifying natural compensatory mutations at sites within this 6 Å radius. Relative to Teleost RH1 (where site 122 varies), we found that sites within this radius displayed decreased amino acid variation in Tetrapod and Characiphysi RH1, where E122 and I122 are fixed, respectively (asterisks, Figure 2C). This observation is consistent with an intramolecular evolutionary process known as entrenchment (Pollock et al., 2012; Goldstein and Pollock, 2017; Shah et al., 2015), where functionally favourable amino acid residues compensating for an original mutation tend to become fixed, thus mutually entrenching favourable amino acids at each position within the coevolving network. We therefore reasoned that if residues at nearby positions are indeed compensatory, then these sites should display a relative decrease in amino acid variation specifically in those vertebrate lineages where an amino acid has been fixed at site 122-- such as E122 in tetrapods and I122 in the Characiphysi. Furthermore, we hypothesized that decreases in amino acid variation observed in these lineages would be driven by an increase in purifying selection on non-synonymous codons, ultimately reflecting the entrenchment of compensatory amino acid residues by natural selection. We therefore employed codon-based phylogenetic likelihood methods to test for a relative increase of purifying selection at RH1 sites within 6 Å of site 122, within Tetrapod vs Teleosts, as well as in Characiphysi vs other Teleosts (Yang, 2007) (Materials and methods). Using likelihood ratio tests of alternative (Clade model C [(Bielawski and Yang, 2004]) and null (M2a_REL ([Weadick and Chang, 2012)]) model analyses of codon substitution rates (dN/dS) across the RH1 coding-sequence, we identified statistically significant evidence of gene-wide increases in purifying selection within Tetrapods (Table 3) and Characiphysi (Table 5) relative to teleosts ((p<0.001)). Sites estimated to be under this increase in purifying selection were those identified in the CmC divergent site class through a Bayes empirical Bayes analysis as previousl