RESEARCH ARTICLE

A Comparative Study Indicates Both Positive and Purifying Selection within Ryanodine Receptor (RyR) Genes, as well as Correlated Evolution PATRICK B. MCKAY AND CORTLAND K. GRISWOLD* Department of Integrative Biology, University of Guelph, Ontario, Canada

ABSTRACT

J. Exp. Zool. 321A:151–163, 2014

Ryanodine receptors are Ca2þ ion channels that allow Ca2þ to flow into the cytosol, from an internal store, in the form of transients. RyRs form a small gene family and vertebrates have three major isoforms, RyR1, RyR2, and RyR3, which are mixed and matched in different combinations in different tissues resulting in different Ca2þ transients for each tissue. In this study, we characterized the interspecies evolution of RyRs within vertebrates. First, we compared the nucleotide divergence of key gene regions including divergent regions (DRs), which are believed to be responsible for the functional divergence between RyRs, and mutation hot‐spot regions, which are responsible for RyR‐related pathologies. We found evidence that DRs undergo positive selection and mutation hot‐spot regions undergo purifying selection. Second, we estimated the extent of purifying selection for RyR1, RyR2, and RyR3 by estimating dN/dS ratios. We found all three genes to be under strong purifying selection, overall. This is consistent with RyRs being used in a diverse set of physiological contexts and therefore under potentially high pleiotropic constraint. Third, we tested for the correlated evolution of dN/dS ratios between RyR genes. We found that RyR2 and RyR3, and RyR3 and a skeletal form of dihydropyridine receptor (DHPR) have correlated rates of evolution. We propose that compensatory effects may explain their correlated evolution. We tested for compensatory function by simulating mutations via a physiological model of RyR function, but did not find evidence for compensation, which indicates that the correlation is likely a result of another process. J. Exp. Zool. 321A:151–163, 2014. © 2013 Wiley Periodicals, Inc. How to cite this article: McKay PB, Griswold CK. 2014. A comparative study indicates both positive and purifying selection within ryanodine receptor (RyR) genes, as well as correlated evolution. J. Exp. Zool. 321A:151–163.

Ryanodine receptors (RyRs) are large Ca2þ ion channels found within the endoplasmic/sarcoplasmic reticulum (ER) membrane. Physiologically, RyRs are important in a diverse set of processes, including muscle contraction, neurotransmitter release, and hormone secretion (Berridge et al., 2000, 2003). In vertebrates, there are three main types of RyR, each encoded by different genes: RyR1, RyR2, and RyR3. These genes form a gene family, arising by gene duplication, after which they diverged functionally (Sutko and Airey, '96; Kushnir et al., 2005). RyR1 and the ancestor to RyR2/RyR3 diverged first while RyR2 and RyR3 share a more recent common ancestor (Darbandi and Franck, 2009). A single type of RyR is used in different tissues. For instance, RyR1

Grant sponsor: Natural Sciences and Engineering Research Council (NSERC) —Discovery; grant sponsor: Canadian Foundation for Innovation; grant sponsor: Ontario Ministry of Research and Innovation.  Correspondence to: Cortland K. Griswold, Department of Integrative Biology, University of Guelph, 50 Stone Rd. E. Guelph, ON, Canada N1G 2W1. E‐mail: [email protected] Received 27 May 2013; Revised 15 September 2013; Accepted 28 October 2013 DOI: 10.1002/jez.1845 Published online 19 November 2013 in Wiley Online Library (wileyonlinelibrary.com).

© 2013 WILEY PERIODICALS, INC.

152 and RyR3 are used in skeletal muscle, neural and pancreatic tissues and RyR2 is used in cardiac, neural, and pancreatic tissues. Therefore, RyRs potentially affect multiple physiological processes and act pleiotropically. Thus, RyRs are a good model system to understand the physiological nature of pleiotropy and to investigate the consequences of pleiotropy on molecular evolution. The genes encoding RyRs occur on different chromosomes. For instance, in humans, RyR1 occurs on chromosome 19, RyR2 on chromosome 1 and RyR3 on chromosome 15. Each gene consists of about 5,000 amino acids and occupy between about 150 kB (RyR1) to 900 kB (RyR2) at the nucleotide level on their respective chromosomes. Two major classes of gene regions within RyR genes have been identified: divergent regions (DRs) and mutation hot‐spot regions. DRs were identified by comparing levels of divergence between RyR genes within the common rabbit [Oryctolagus cuniculus] (Hakamata et al., '92). In total, three DRs were found and labeled DR1, DR2, and DR3 (Sorrentino and Volpe, '93). Although RyR1 and RyR2 contain all three DRs, RyR3 lacks DR2. Studies have shown that these regions are accountable for many of the functional differences between RyR paralogs (Liu et al., 2004; Hamilton, 2005). For example, DR1 is responsible the inactivation of RyR1 and RyR3, but not RyR2 (Liu et al., 2004). DR2 is responsible for different sensitivities to Ca2þ activation (Hamilton, 2005). DR3 allows for the specialized binding of RyRs to the appropriate type of FK506 binding protein (FKBP), a regulatory RyR‐binding protein (Liu et al., 2004). Furthermore, DR3 contains a binding site within RyR1 for the dihydropyridine receptor (DHPR). DHPR is a Ca2þ ion channel found within the cell membrane. We propose that DRs are responsible for the functional divergence of individual RyRs between species, in addition to the functional divergence of RyR paralogs. Since nucleotide variability between species is indicative of adaptive evolution (Nielsen, 2005), we predict that nucleotide divergence within a DR region between species will be relatively high compared to the average divergence across the entire coding region of a gene. Mutation hot‐spot regions are gene regions within RyRs that contain a high concentration of mutations linked to human diseases including malignant hyperthermia and central core disease for RyR1 and catecholaminergic polymorphic ventricular tachycardia and arrhythmogenic right ventricular dysplasia for RyR2 (Yano et al., 2006; George et al., 2007; Betzenhauser and Marks, 2010; Mackrill, 2010; Lanner, 2012). Mutation hot spots within RyR3 have not yet been discovered. For both RyR1 and RyR2, one mutation hot‐spot region is located in the N‐terminus, one is located in the central region and two are located at the C‐ terminus of the gene, one of which overlaps with the transmembrane pore (George et al., 2007). It has been suggested that other vertebrates are susceptible to pathologies caused by mutations within these regions. For instance, Hirata et al. (2007) found that RyR1b mutations for zebrafish (Danio rerio) lead to a form of central core disease. Since these mutations can be highly J. Exp. Zool.

MCKAY AND GRISWOLD deleterious, we expect mutation hot‐spot regions to undergo strong purifying selection and for the interspecies nucleotide divergence for these regions to be relatively low. RyRs are involved in a process called calcium‐induced‐calcium release (CICR). During this process, a small influx of Ca2þ into the cytosol can induce a large amplitude Ca2þ transient or wave. The larger transient or wave of Ca2þ is brought about by the initial small‐amplitude influx of Ca2þ binding RyR channel receptors, inducing the channel to open and releasing the store of Ca2þ from the endo/sarcoplasmic reticulum. The nature of a Ca2þ transient or wave distinguishes different cell types and physiological processes. For instance, whereas in pancreatic tissue Ca2þ signaling involves protracted waves on a characteristic time scale of minutes (Sasaki et al., 2002), in cardiac tissue, signaling is characterized by short transients on a characteristic timescale of a second (e.g., Rice et al., '99; Sobie et al., 2002). Similarly, signaling in skeletal muscle is different from neural tissue, which is different from pancreatic and cardiac tissue. The diversity of Ca2þ signaling between tissues is caused, in part, by each RyR being functionally distinct and each tissue having different compositions of RyR channels (Berridge et al., 2000, 2003). The functional distinction between RyRs can be characterized, in part, by their propensity to be opened by Ca2þ. RyR1 and RyR2 are characterized by requiring a lower Ca2þ concentration to induce opening than RyR3. RyR1 and RyR2 are then distinguished from each other by RyR1 having a lower tendency to remain open compared to RyR2. Although RyR3 requires a higher threshold of Ca2þ to induce opening, it has a strong tendency to stay fully open for a short period of time before being induced to close by elevated Ca2þ in the cytosol (see Griswold [2011] for a summary of empirical work supporting these properties). Furthermore, RyR1 is distinguished from RyR2 and RyR3 in that it physically couples with DHPR in skeletal muscle (Paolini et al., 2004) and some neural tissue (Mouton et al., 2001; De Crescenzo et al., 2006). This physical coupling allows for channel opening in the absence of an influx of Ca2þ. These functional differences between RyRs allow RyRs to be mixed and matched between tissues, resulting in the appropriate phenotype for that tissue (Conti et al., '96; Sutko and Airey, '96). For example, RyR1 and RyR3, but not RyR2, are expressed in skeletal muscle tissue, RyR2 is expressed exclusively within cardiac tissue and RyR1, RyR2, and RyR3 are each expressed in neural and pancreatic tissues. Furthermore, distinct skeletal and smooth muscles express RyR1 and RyR3 in different proportions and distinct brain regions express RyR1, RyR2, and RyR3 in different proportions (Conti et al., '96; Mori et al., 2000; Mouton et al., 2001). Although RyRs can be mixed and matched between tissues to generate functional divergence, these tissues nevertheless share RyRs. Shared (pleiotropic) RyRs can cause correlations among tissues, when polymorphic (Griswold, 2011). For instance, mutations to RyR2 in mice lead to cardiac arrhythmias and epileptic

DIVERGENCE AND CORRELATED EVOLUTION OF RyRs seizures (Lenhart et al., 2008). One consequence of a gene acting pleiotropically is that it may be constrained from evolving (Arnold, '92; Hansen, 2003; Otto, 2004). For example, a mutation in RyR1 may be beneficial for skeletal muscle tissue, but deleterious for neural tissue causing RyR1 to evolve at a slower rate. Given that RyRs are pleiotropic, we expect to find evidence that RyRs undergo stronger purifying selection than other genes. Evidence for this would be exceptionally low dN/dS ratios for RyR genes compared to other genes. A related question is whether or not one RyR is under stronger pleiotropic constraint than the other two. For instance, RyR2 may be under stronger purifying selection than RyR1 and RyR3 as it is, to our knowledge, the only RyR subtype expressed in cardiac muscle tissue and a mutation that is beneficial for RyR2 in neural tissue may have a deleterious side effect in the heart. We therefore compare dN/dS ratios for RyR1, RyR2, and RyR3. Furthermore, we hypothesize that RyRs undergo correlated rates of evolution because of their integrated physiology generating tissue‐specific calcium transients or waves. We test this hypothesis by correlating dN/dS ratios between RyR gene pairs, as well as with the skeletal muscle copy of DHPR (encoded by the CACNA1S gene) across branch sections of a phylogeny. We use the dN/dS ratio as our evolutionary rate for two reasons. First it is a standardized measure of the rate of nonsynonymous substitution, thus capturing evolutionary processes subject to selection. Second, because it is a standardized measure of measure of evolutionary rate, it accounts for spurious processes such as phylogenetically conserved mutation rates or other hidden variables that may otherwise generate a correlation. Three evolutionary processes could cause correlated evolution between RyRs. The first process is nonadaptive and involves correlated sequence evolution between genes due to a shared hidden variable between genes. For instance, the genes share in common the effective population size of a lineage. If a lineage is characterized by a small effective population size, then all RyR genes may experience elevated rates of nearly neutral substitutions. In our analysis we control for this using a generalized least square approach that accounts for phylogeny. A second process is a pair of genes may simultaneously experience positive selection. A third process leading to a correlation is compensatory evolution. Given that RyRs act pleiotropically, compensatory evolution may be applicable. Generally, a compensatory mutation is a nonadaptive mutation that is subject to positive selection and restores a phenotype towards its original value after the occurrence of a deleterious mutation in another gene (Camps et al., 2007). In this study, we use a broader interpretation of compensatory mutations, where a mutation that is beneficial overall but has a deleterious pleiotropic effect in a particular tissue can be compensated by a secondary mutation in another gene that partially or completely restores the original phenotype in that tissue (Pavlicev and Wagner, 2012). This secondary mutation is said to be compensatory. For example, a

153 mutation in RyR2 may be beneficial overall due to its effects on heart function, yet have a deleterious effect on neural function. Given its overall benefits, the mutation may fix. Once fixed, however, there is the opportunity for RyR1 and RyR3, both of which are also found in neural tissue, to compensate for the change in physiology in neurons via mutation. These mutations are then expected to be selected for and be fixed within a population. We therefore expect that the rates of evolution between these genes would be correlated. In addition, because DHPR binds directly to RyR1 in skeletal muscle, we hypothesize that the skeletal muscle copy of DHPR and RyR1 has correlated rates of evolution. We test for a compensatory explanation for a correlation between genes using a model developed by Griswold (2011) of Ca2þ signaling in skeletal muscle, neural and heart tissue. The Griswold (2011) paper builds on models of Ca2þ signaling that focused on cardiac tissue by Rice et al. ('99) and Sobie et al. (2002), as well as more recent empirical work that characterizes the physiology of Ca2þ signaling. Key experimental results used to parameterize and verify the model include open probabilities for RyR1 (Chu et al., '93; Morrissette et al., 2000; Percival et al., 2004), RyR2 (Gyorke and Fill, '93; Gyorke et al., '94; Li and Chen, 2001), and RyR3 (Chen et al., '97). The model also accurately predicted that the addition of RyR3 increases the amplitude of Ca2þ transients, that there is a negative correlation between transient amplitude and the half‐length of a transient and that the amplitude of cardiac transients tends to be greater than skeletal muscle transients (e.g., Klein et al., '96; Weisleder et al., 2007). A more thorough explanation of the model is given the Supplementary File 4. Although the total number of RyRs remains constant, different tissue types in the model are characterized by different compositions of RyR channels. In the model there are parameters that govern the biophysical properties of RyRs, as well as DHPR, such as the rate at which a channel opens and the flow rate of Ca2þ through a channel. The different RyR channels share the same set of parameters, but differ in their respective values because they each have different biophysical properties (ref. Table 2, Griswold [2011]). The different biophysical properties of RyRs and DHPR and the different compositions of RyR channels between tissues result in different Ca2þ transient properties between tissues. In this analysis we mimic mutations in RyRs and DHPR by perturbing these parameters in the model and test whether a particular gene is best at compensating for perturbations in another gene.

MATERIALS AND METHODS Sequence Data Coding sequences for RyRs and DHPR were retrieved from GenBank (http://www.ncbi.nlm.nih.gov/genbank/) through NCBI (see Supplementary Table S1 for sequences used and accession numbers). All possible sequences within vertebrates were considered, but coding sequences with unusually large gaps were J. Exp. Zool.

154 suspected of being incomplete and excluded from the study. Sequences were aligned uusing ClustalX 2.1 (Larkin et al., 2007) with default parameters. Because a few sequences within each alignment had insertions that resulted in frameshift mutations, stop codons sometimes appeared within the aligned genes. We corrected for this by removing these insertions, as long as the sequence alignments were properly maintained. When this correction did not keep the sequences properly aligned, stop codons were accounted for by converting the nucleotides within stop codons to gaps. Nucleotide Divergence We tested for the best substitution model using MEGA5 (Tamura et al., 2011); Tamura and Nei's ('93) model, with nonuniform evolutionary rates and a fraction of sites that are invariable, was found to be the best. We accounted for the nonuniform evolutionary rates and the fraction of invariable sites by separating genes into exons, DR, mutation hot‐spot and transmembrane regions. We estimated and report Tamura–Nei (TN93) distances (Tamura and Nei, '93) for each region. In addition, we estimated Kimura 2‐parameter (K2) distances (Kimura, '80) and found them to be very similar to the TN93 values, on average (not reported). Both TN93 and K2 estimates were done using custom code in Mathematica (Wolfram Research, Inc., 2008). After aligning the sequences, we calculated pairwise distances between all species, while excluding gaps and ambiguous characters within the sequence alignments. We used resampling to obtain 95% confidence intervals of mean estimates. In particular, we sampled with replacement 100 pairwise differences 10,000 times, giving us 10,000 estimates of the mean TN93 distance for each region. Nucleotide positions spanning the DRs for all three genes were obtained from Liu et al. (2004). Only the DRs for RyR1 and RyR2 were given, so we used the corresponding homologous regions within RyR3 to calculate divergence values for DR1 and DR3 (RyR3 lacks DR2). Mutation hot‐spot regions span four distinct regions of RyR1 and RyR2 that have high concentrations of disease causing mutations. Start and end positions for these regions were obtained from Yano et al. (2006) and George et al. (2007). Estimating Rates of Evolution We preformed a maximum likelihood analysis with CODEML in the PAML 4.4 package to estimate the ratio of nonsynonymous changes per site (dN) to synonymous changes per site (dS; Yang, 2007). Multiple alignments were fitted to the F3  4 model of codon frequencies. We used the most recently updated phylogeny of vertebrates for our analysis (Meredith et al., 2011). A polytomy was allowed for branch sections for bovidae (Bos taurus), carnivora (Ailuropoda melanoleuca and Canis lupus), and equidae (Equus caballus) as their phylogeny is still disputed. Zebrafish was used as an outgroup for each phylogenetic tree. J. Exp. Zool.

MCKAY AND GRISWOLD When comparing levels of purifying selection between RyRs, we ran a model allowing each RyR to evolve at different rates. Using a log‐likelihood ratio test, we then compared this model to a null model where RyRs were not allowed to evolve at different rates. Our comparison of dN/dS ratios between genes is based on sequences from species that are shared between all RyR genes. To test for the correlated evolution of RyRs, we estimated correlations between dN/dS estimates for each pair of RyRs and between RyR and DHPR genes. For each correlation, we constructed a phylogenetic tree using all species that had sequences available for both genes (Fig. 1). As Figure 1 indicates, these correlations involved different subsets of species between gene pairs. We used this approach to maximize power to detect a correlation for a particular gene pair. Since we used different subsets of species, direct comparison of the magnitude of a relationship between species is not valid. Isoform RyR1a in zebrafish was used as an outgroup for RyR1 and isoform RyR2b in zebrafish was used as an outgroup for RyR2. For each correlation, we estimated dN/dS ratios for both genes for each branch of the phylogenetic tree using a free ratio model in PAML, which allows for independent dN/dS estimates for each branch. While discounting outgroups, we performed a linear regression on the log10 dN/dS estimates for each branch between RyR genes. Branches with dN/dS estimates of infinity, which are a result from a dS estimate of zero, and branches with dN/dS estimates of zero, which are a result of a dN estimate of zero, were also excluded from the study. We observed that each of the longer branches have dN/dS estimates less than 1 whereas some of the shorter branches have unreasonably high dN/dS estimates, possibly because shorter branch lengths have more sampling error as there are fewer codons undergoing mutations. Short branches with dN/dS estimate larger than 2 were therefore excluded. We fitted three regression models to the data and then compared how well they fit the data via the corrected Akaike information criterion (AICc), which accounts for finite sample sizes (Akaike, '74). We selected the model that gave a significantly more negative estimate. The first regression model was from Martins and Housworth (2002), which accounts for the shared common ancestry between species and unequal branch lengths. One can imagine a confounding variable—for example, effective population size (Ne)— resulting in a bias for dN/dS estimates between clades. For example, small mammals may have an Ne that affects dN/dS estimates a certain way, whereas large mammals may have an Ne that affects dN/dS estimates another way. This could result in large estimates for paired genes in one clade and lower estimates for paired genes in for another and lead to a positive correlation even if the evolution of both genes is not correlated. Under this model, the slope of the regression is estimated by applying a matrix to the data that accounts for shared common ancestry and shared branch lengths. We used the mean estimated divergence times, obtained from TimeTree, to calculate branch lengths (Hedges et al., 2006). The second model accounted for unequal branch lengths, but did

DIVERGENCE AND CORRELATED EVOLUTION OF RyRs

155

Figure 1. Phylogenetic trees used for comparing dN/dS ratios between RyR1 and RyR2 (a), RyR1 and RyR3 (b), RyR2 and RyR3 (c), RyR1 and DHPR (d), RyR2 and DHPR (e), and RyR3 and DHPR (f). Taxa used in this section include zebrafish (Danio rerio), bullfrog (Rana catesbeiana), lizard (Anolis carolinensis), turkey (Meleagris gallopavo), chicken (Gallus gallus), horse (Equus caballus), bull (Bos taurus), panda (Ailuropoda melanoleuca), dog (Canis lupus), boar (Sus scrofa), mouse (Mus musculus), rat (Rattus norvegicus), rabbit (Oryctolagus cuniculus), marmoset (Callithrix jacchus), rhesus (Macaca mulatta), orangutan (Pongo abelii), chimp (Pan troglodytes), and human (Homo sapiens). In this figure, the RyR1a gene is used for RyR1 and the RyR2b gene is used for RyR2 in zebrafish.

not account for shared common ancestry. Longer branches have smaller variances in estimates because there are a higher number of nucleotide substitutions between sequences, resulting in a low sampling error when estimating dN/dS. We accounted for this

variance by applying a transformation‐matrix to the data giving more weight to the points with longer branch lengths and assuming off‐diagonal elements are equal to zero. The third model involved an ordinary least squares (OLS) approach where we assumed equal J. Exp. Zool.

156

MCKAY AND GRISWOLD

variance in residuals and independence for all branches. Our AICc comparisons for each plot revealed that the third model fits the data significantly better than the other two (Akaike, '74). We therefore used the OLS model in our analysis. The significance of slopes and correlations was measured by conducting a two‐sided t‐test. Testing for Compensation We conducted simulations for seven types of neural tissue, approximating the composition of various areas of the brain found by Mori et al. (2000). We also conducted simulations for three types of skeletal muscle tissue, two of which are predominantly composed of RyR1 with low amounts of RyR3, as observed in nature (Conti et al., '96). The third simulation for skeletal muscle had equal amounts of RyR1 and RyR3. To measure the effects of mutations, we compared simulations with a single perturbed parameter and two perturbed parameters to a simulation without perturbed parameters and whose values are equal to those given in Tables 1 and 2 of Griswold (2011). Mutations within RyRs were simulated by either increasing or decreasing parameter values by 10%. Mutations within DHPR were simulated by increasing or decreasing the parameter that corresponds to the DHPR voltage to RyR1 channel opening constant by an order of magnitude. Recorded phenotypes include the average amplitude, duration, integrated concentration, and D50 of a calcium transient. The integrated concentration of a calcium transient integrates the amplitude of a transient over its duration and D50 is the duration of a transient between when its amplitude reaches one‐half of its maximum value at the front‐end of the transient to when it reaches one‐half of its maximum amplitude at the tail‐end of the transient. One hundred transient replicates were simulated per parameter set. Ninety‐five percent confidence internals of average estimates were obtained by resampling with replacement 100 times. To test for compensation, we used the following equation: S¼

jf 12  f 0 j jf 1  f 0 j

ð1Þ

where f12 is the mean phenotype for transients with two mutations, f1 is the mean phenotype for transients with one mutation, and f0 is the mean phenotype for transients with no

Table 2. Linear regression analysis comparing dN/dS ratios between RyR and DHPR gene pairs. Gene pairs

Slope (P‐value)

Correlation (P‐value)

DHPR–RyR1 DHPR–RyR2 DHPR–RyR3

0.0912 (0.5717) 0.0930 (0.8241) 0.7456 (0.0389)

0.1731 (0.5717) 0.0531 (0.8241) 0.5200 (0.0389)

Estimates of the slope and correlation for dN/dS estimates between RyR genes and the DHPR gene. Significant estimates are in bold.

mutations. Statistics greater than 1 suggest that the secondary mutations are decompensatory, statistics less than 1 suggest that the secondary mutations are compensatory and statistics equal to 1 suggest that the secondary mutations are neutral. After verifying for normality via the Kolmogorov–Smirnov test, we tested to see if test statistics were larger or smaller than 1 at the 5% significance level using a t‐test. We calculated statistics for the amplitude, duration, integrated concentration, and D50 for every possible set of mutations for all simulated tissues. An example of this calculation is found in Supplementary File S4. Once we calculated all test statistics, we counted the number of compensatory and decompensatory mutations for each RyR pair and calculated the corresponding proportions (Table 3). We then compared these proportions to see if any particular RyR pair has higher levels of compensatory function using a two‐sided t‐test. Additional analyses investigated whether RyRs compensate for different types of mutations (see Supplementary File 4).

RESULTS Nucleotide Divergence A comparison of nucleotide divergence between genes using a common set of species with gene sequences for all three RyR genes indicates that RyR2 has the most interspecies divergence, and RyR3 the least interspecies divergence. Specifically, the 95% confidence interval of Tamura–Nei (TN93) divergence for RyR1 is {0.183, 0.188}, for RyR2 is {0.219, 0.225}, and for RyR3 is {0.176, 0.181}.

Table 1. Linear regression analysis comparing dN/dS ratios between RyR gene pairs. Whole gene

Divergent regions

RyR gene pairs

Slope (P‐value)

Correlation (P‐value)

Slope (P‐value)

Correlation (P‐value)

RyR1–RyR2 RyR1–RyR3 RyR2–RyR3

0.5015 (0.1676) 0.1537 (0.0689) 0.5965 (0.0142)

0.4257 (0.1676) 0.5417 (0.0689) 0.6372 (0.0142)

0.6290 (0.1426) 0.5100 (0.1041) 0.4255 (0.2272)

0.5296 (0.1426) 0.5767 (0.1041) 0.4198 (0.2272)

Estimates of the slope and correlation for dN/dS estimates between whole RyR genes and divergent regions. Significant estimates are in bold.

J. Exp. Zool.

DIVERGENCE AND CORRELATED EVOLUTION OF RyRs

157

Table 3. Testing for compensation and decompensation between RyR and DHPR gene pairs.

Tissue

RyR pairs

Proportion that are compensatory

Proportion that are decompensatory

P‐values

Neural

RyR1–RyR3 RyR1–RyR2 RyR2–RyR3 Combined RyR1–RyR3 RyR1–DHPR RyR3–DHPR Combined

0.02995 0.04152 0.04375 0.03757 0.02083 0.07917 0.09583 0.08750

0.04268 0.06920 0.07634 0.06120 0.06750 0.36250 0.35000 0.35625

1.7526  104 3.8561 109 2.3585  1011 1.1856  1020 2.3917 1015 1.7076  1015 2.1328  1012 3.3440  1026

Skeletal

Proportion of mutations that are compensatory and decompensatory between RyR pairs in neural and skeletal tissues and between RyR and DHPR in skeletal muscle tissue. Proportions include mutations that affect transient amplitude, duration, integrated concentration, and D50. P values test for a difference between rates of compensation versus decompensation.

Figure 2 provides a picture of nucleotide divergence across exons within a gene for RyR1, RyR2, and RyR3. It also serves as a point of comparison when considering divergence within DR, mutation hot spot and transmembrane regions within a gene. The figure indicates a large amount of variation in divergence among exons within a gene. It is important to note that in Figure 2, different sets of species were used for each gene when estimating nucleotide diversity, so a comparison between genes within this figure is not valid. Within gene nucleotide divergence using a set of species that is shared between all genes (Supplementary Fig. S2) has a similar pattern as in Figure 2. Interspecies nucleotide divergence within DRs is significantly greater than the average divergence within a gene (Fig. 3). Mutation hot‐spot regions have divergence values that are significantly lower than the average divergence within a gene, except for the N‐terminal mutation hot‐spot region for RyR2. We found the transmembrane region for all three RyRs to have larger divergence values than the average divergence within a gene. The pattern of divergence between different regions is the same if sequences from a shared set of species between genes are used (Supplementary Fig. S3). Rates of Evolution of RyR Genes We found that a model that assumes different rates of evolution between RyR genes is better than a model that assumes equal rates between RyRs (P‐value ¼ 6.52  1019). The dN/dS ratio for RyR1 is 0.0223, for RyR2 it is 0.0176, and for RyR3 it is 0.0371. This analysis is based on the shared set of species between genes, so direct comparisons of dN/dS are valid. Testing for the Correlated Evolution of RyR Genes We found that dN/dS ratios for RyR2 and RyR3 are positively correlated across phylogenies (Table 1, Fig. 4). dN/dS ratios for RyR1 and RyR3 do not have a significant positive relationship,

although it is close to significance at the 5% significance level. There is no evidence that RyR1 and RyR2 undergo correlated evolution. Furthermore, there is no evidence that DRs undergo correlated evolution (Table 1). DHPR and RyR3 have a significant positive relationship in dN/dS ratios across lineages, whereas DHPR and RyR1 and DHPR and RyR2 do not (Table 2, Fig. 5). Testing for Compensation Between RyRs Using Simulations We tested for compensation by analyzing the proportions of compensatory mutations between RyR pairs. Generally, we did not find evidence for stronger compensation between particular RyR pairs for any transient property (Table 3). We tested for stronger compensation between particular RyR pairs by comparing the proportion of compensatory mutations between RyR pairs. For instance in the case of transient amplitude, we did not find the propensity to be compensatory to significantly vary between RyR pair [i.e., compensatory values given in Table 3 were not significantly different] (RyR1/RyR2 vs. RyR1/RyR3: P‐value ¼ 0.2547, RyR1/RyR2 vs. RyR2/RyR3: P‐value ¼ 0.0972, and RyR1/RyR3 vs. RyR2/RyR3: P‐value ¼ 0.6026). We also tested for compensatory function between RyR1 and DHPR and RyR3 and DHPR in skeletal muscle tissue but did not find significant evidence for stronger compensation between RyR1 and DHPR when compared to RyR3 and DHPR (P‐value ¼ 0.5278). These results are not consistent with the correlations between dN/dS ratios because there is no evidence that RyR3 is better at compensating for RyR2 and/or DHPR, or vice versa.

DISCUSSION Nucleotide Divergence of RyR Genes The high interspecies divergence found within each of the DRs suggests that these regions are under either relaxed purifying J. Exp. Zool.

158

MCKAY AND GRISWOLD

Figure 2. Bar charts depicting the average and 95% confidence intervals of the Tamura–Nei distance (DTN93) for all exons within RyR1 (a), RyR2 (b), and RyR3 (c), using a different set a species for each gene. The average TN93 distance across all species for the entire gene is represented by the horizontal line.

Figure 3. Bar charts depicting the average and 95% confidence intervals of the Tamura–Nei distance (DTN93) of divergent regions (DRs), mutation hot‐spot regions (MHRs), and transmembrane (TM) regions within RyR1 (a), RyR2 (b), and RyR3 (c), using a different set a species for each gene. The average TN93 distance across all species for the entire gene is represented by the horizontal line.

selection or greater positive selection. Given that prior studies indicate that DRs are responsible for functional divergence between RyR1, RyR2, and RyR3 (Liu et al., 2004; Hamilton, 2005), it is likely that the high interspecies divergence is driven by positive selection as opposed to relaxed selection. Our results provide further support that mutation hot‐spot regions are unique with RyR genes. Prior studies have found that

these regions have a relatively high density of disease‐associated mutations (Betzenhauser and Marks, 2010; Mackrill, 2010; Lanner, 2012). Our results indicate that mutation hot‐spot regions are under stronger purifying selection compared to other regions in a gene. The only mutation hot‐spot region that has higher divergence than average was the N‐terminus region of RyR2, which includes the first exon of each gene and is likely an artifact

J. Exp. Zool.

DIVERGENCE AND CORRELATED EVOLUTION OF RyRs

Figure 4. Plots of log10 dN/dS estimates between the entire genes of RyR1 and RyR2 (a), RyR2 and RyR3 (b), and RyR1 and RyR3 (c) with a simple linear regression line fitted. The comparison between RyR2 and RyR3 is the only one with a significant positive slope and correlation at the 5% significance level (Table 1).

of the sequence alignment. Generally, our results lend support to the idea that disease associated mutations arise out of functionally conserved regions of a gene. Future studies may examine intraspecies variability by calculating nucleotide divergence at the population level. When positive directional selection is occurring, the ratio of interspecies

159

Figure 5. Plots of log10 dN/dS estimates between DHPR and RyR1 (a), RyR2 (b), and RyR3 (c) with a simple linear regression line fitted. The comparison between DHPR and RyR3 is the only one with a significant positive slope and correlation at the 5% significance level (Table 2).

to intraspecies variability is increased, whereas when purifying selection is taking place, this ratio is reduced (Nielsen, 2005). In our study, the DRs, which appear to undergo stronger positive directional selection, are expected to have a higher ratio of interspecies to intraspecies variation than the mutation hot‐spot regions, which appear to undergo stronger purifying selection. J. Exp. Zool.

160 Evolutionary Rates of RyR Genes The dN/dS estimates indicate that RyRs undergo purifying selection, on average, which is not surprising given that very few genes have dN/dS values greater than one as most genes undergo purifying selection (Kimura and Ohta, '74). However, we suspect that the pleiotropy of RyRs make it undergo stronger purifying selection than other genes because mutations with beneficial phenotypes in one tissue type are likely to have deleterious side effects in others. A study involving yeast conducted by Salathé et al. (2006) found that the number of physiological processes that a gene takes part in is inversely correlated with its rate of evolution. When comparing our results with the results of this study, our dN/dS estimates for RyR1, RyR2, and RyR3 fall below the 95% confidence intervals of genes that are involved in five or fewer cellular processes. This comparison must be taken lightly, however, since we only examined animal genes within our study whereas Salathé et al. (2006) examined genes within fungi. When examining the literature, we found that our dN/dS estimates were lower than most dN/dS estimates for other animal genes. For example, the dN/dS ratio for all RyRs fell below the 95% confidence intervals of the mean dN/dS ratios for Allocentrotus fragilis (purple sea urchin) and Strongylocentrotus purpuratus (pink sea urchin) genes that were sampled by Oliver et al. (2010). We found that RyR3 undergoes the weakest purifying selection, whereas RyR2 undergoes the strongest. These results are interesting because RyR3 is the most recently derived RyR gene (it is likely RyR2 retains its ancestor's function) and physiologically RyR3 is typically the gene that is expressed at lowest frequency within a tissue and does not occur by itself within any tissue type, to our knowledge. In contrast, RyR2 occurs by itself in cardiac tissue and RyR1 is the predominant gene in adult skeletal muscle tissue. Because of the prevalence of RyR1 and RyR2 within these tissues, pleiotropic effects may be more severe in these tissues and mutations would therefore generally be subject to stronger purifying selection. The finding that RyR2 has higher interspecies divergence per site, on average, compared to RyR1 and RyR3, but a significantly lower dN/dS estimate than RyR1 and RyR3 is interesting. This is suggestive that RyR2 undergoes a higher rate of synonymous substitution than RyR1 and RyR3. This result is inconsistent with what is expected under background selection conditions in which it is expected that purifying selection at functional sites suppresses the fixation of linked neutral variation (Charlesworth et al., '95). Correlated Evolution Between Genes All significant positive relationships in our analysis of dN/dS ratios involved RyR3, be it with RyR2 or DHPR. The AIC analysis of alternate regression models indicates that hidden variables, such as effective population size (Counterman et al., 2004) that exhibit phylogenetic structure are not responsible for these positive relationships. Furthermore, our computer model does not suggest J. Exp. Zool.

MCKAY AND GRISWOLD that RyR3 is better at compensating for mutations in RyR2 or DHPR compared to other genes. Therefore, our results indirectly point to positive selection driving the correlation between RyR3 and RyR2, and RyR3 and DHPR. It should be noted that our model is restricted to the dynamical properties of Ca2þ transients and does not model direct molecular interactions between molecules. For instance, our model does not account for compensation that may occur between the RyR1 and DHPR as a result of their direct physical coupling. It is interesting to note that RyR2 and RyR3 are sister gene duplicates and have similar biophysical properties. In particular, their channels both tend to have a high propensity to remain open, once induced to open. There are important functional differences between RyR2 and RyR3, such as RyR3 inactivation by Ca2þ and RyR3 having a higher threshold to open, yet their functional similarities may result in them both responding evolutionary to a selective pressure. It is also surprising that this correlated evolution was not detected between the DRs of RyR2 and RyR3, since DRs are functionally important and we would expect positive selection to act on these regions. The correlated evolution between RyR3 and DHPR is surprising. A priori we expected that RyR1 and DHPR would undergo correlated evolution because they physically interact. In skeletal muscle tissue the two‐dimensional arrangement of RyR1, RyR3, and DHPR molecules appears to be important in modulating the physiology of muscle contraction (Yin et al., 2008). In particular, RyR1, RyR3, and DHPR form a lattice network with DHPR clustered in tetrads that are physically coupled to a single RyR1 molecule. Importantly, depolarization of DHPR not only directly induces opening of the channel of the RyR1 it is directly coupled to, but also indirectly induces opening of neighbor RyR channels, which may include RyR3. Therefore, there could be indirect selection that drives the correlated evolution of RyR3 and DHPR. There may also be an as yet unrecognized direct interaction between the two molecules. A factor that may drive a correlation between genes involved in calcium signaling is that our study included sequences from both ectotherms and endotherms. Ectotherms must maintain physiological functions over a larger temperature range. These larger temperature ranges can potentially lead to a larger range of intracellular environments for the organism in question. The genes involved in physiological processes must therefore be able to function properly for each of the numerous temperature‐ dependent environments, thus potentially constraining its ability to evolve. This is similar to a pleiotropic constraint. Consequently, RyRs within ectotherms may be under stronger purifying selection and have relatively lower dN/dS ratios than endotherms, thus resulting in a positive correlation. To the best of our knowledge, there are no studies that have compared the levels of purifying selection between endotherms and ectotherms. However, we found no evidence for this hypothesis as ectothermic lineages within our study had dN/dS ratios that were within the range of

DIVERGENCE AND CORRELATED EVOLUTION OF RyRs what was normal for RyRs. In addition, one may argue that ambient temperatures may affect dN/dS ratios. Studies have found that temperature is positively correlated with substitution rates within genes, but these studies have found no significant difference in dN/dS ratios between species living in different climates (Gillman et al., 2009; Wright et al., 2011). Another factor to consider concerns concerted evolution, where paralogous genes share greater homology than orthologous genes (Liao, '99; Chen et al., 2007). One can imagine a hypothetical scenario where RyR2 and RyR3 in humans are more similar to each other than RyR2 in humans and RyR2 in mice, even if the divergence of RyR2 and RyR3 occurred prior to the speciation event of humans and mice. This gene conversion may occur through a recombination event of two genes within a species. Although the functional constraints of RyRs make such a conversion difficult to imagine, this may have occurred for a small region within an RyR gene. Once this region is introduced into another RyR gene, the same mutations would be found within both genes, possibly driving their correlated evolution. Although we did not find evidence for the compensatory function of RyRs, it is still possible that compensation is important in the evolution of RyRs. For instance strong selection pressure on weak compensatory effects, which may not be significant enough to be detected by the computer model, may be responsible for the correlated evolution of RyRs. For instance, it is also possible that some of the simulated mutations are more likely to occur in nature than others and that if we looked exclusively at the most likely mutations at the level of the computer model, we would find evidence for compensation that is consistent with the correlated rates of evolution. It is possible that compensation may occur by changing the various concentrations of RyRs within tissues. For instance, a mutation in RyR1 that has a deleterious pleiotropic effect in skeletal muscle may be compensated by either increasing or decreasing the levels of RyR1 and/or RyR3 in skeletal muscle tissue. As was previously mentioned, different expression levels of each RyR are found across tissues (Conti et al., '96; Mori et al., 2000; Mouton et al., 2001). This lends support to the idea that there is the potential for compensation by changing expression levels of RyRs. Furthermore, Ca2þ signaling genes other than RyRs or DHPR may have evolved to compensate for RyRs. In particular, a future study may focus on the inositol triphosphate receptor (IP3R) gene family, a set of Ca2þ ion channels that have a similar sequence homology and function as RyRs (Yamamoto‐Hino et al., '94; Berridge et al., 2000, 2003). Similar to RyRs, there are three major isoforms of IP3Rs, each of which are encoded by their own genes. It is possible that compensatory mutations within these IP3Rs can restore cytoplasmic and ER Ca2þ concentrations to their original levels. A future study may test for the correlated evolution of IP3R isoforms, as well as test for the correlated evolution of individual IP3R and RyR isoforms. Other studies may include other Ca2þ signaling molecules that influence Ca2þ levels within the

161 cytoplasm such as SERCA pumps (Berridge et al., 2000, 2003) and nonpleiotropic genes. Compensatory mutations may also occur within the many accessory proteins that bind to and regulate RyR activity. In particular, it would be worth looking at genes that bind exclusively to RyRs. Pavlicev and Wagner (2012) argue that these “private genes” would be better candidates for compensation because mutations within these genes would result in fewer pleiotropic effects. Furthermore, one may expect that a particular RyR that compensates for a mutation at a different RyR would have a stronger functional resemblance to the other RyR once the compensatory mutation is selected for. This would make it more difficult for RyRs to diverge and specialize for their respective tissues. This argument therefore lends support to the prediction that if the compensation of RyRs occurs, it would occur in the accessory proteins or proteins that interact with RyRs and not RyRs themselves.

CONCLUSIONS Our study lends support to the assertion that regions within RyR genes that underlie functional differences between genes also underlie functional differences between species for the same gene. Our finding that mutation hot spots have comparatively low DTN93 distances indicates these are distinct regions undergoing stronger purifying selection. Our result that dN/dS ratios for RyR genes are exceptionally low is consistent with the observation that these genes are used in a diverse set of physiological processes and therefore potentially under strong selective constraint, overall. The finding that the evolution of RyR3 is correlated with RyR2 and DHPR indicates new lines of inquiry into positively selected processes that may cause the correlation.

ACKNOWLEDGMENTS We thank Jinzhong Fu, Brad Hanna and two anonymous reviewers for their comments. This study was funded by Natural Sciences and Engineering Research Council (NSERC)—Discovery, Canadian Foundation for Innovation, and Ontario Ministry of Research and Innovation.

LITERATURE CITED Akaike H. 1974. A new look at the statistical model identification. IEEE Trans Automat Contr 19:716–723. Arnold SJ. 1992. Constraints on phenotypic evolution. Am Nat 140:85–107. Berridge MJ, Lipp P, Bootman MD. 2000. The versatility and universality of calcium signaling. Nat Rev Mol Cell Biol 1:11–21. Berridge MJ, Bootman MD, Roderick HL. 2003. Calcium signaling: dynamics, homeostasis and remodeling. Nat Rev Mol Cell Biol 4:517–529. Betzenhauser MJ, Marks AR. 2010. Ryanodine receptor channelopathies. Pflugers Arch 460:467–480. Camps M, Herman A, Loh E, Loeb LA. 2007. Genetic constraints on protein evolution. Crit Rev Biochem Mol Biol 42:313–326. J. Exp. Zool.

162 Charlesworth D, Charlesworth B, Morgan MT. 1995. The pattern of neutral molecular variation under the background selection model. Genetics 141:1619–1632. Chen SRW, Li X, Ebisawa K, Zhang L. 1997. Functional characterization of the recombinant type 3 Ca2þ release channel (ryanodine receptor) expressed in HEK293 cells. J Biol Chem 272:24234–24246. Chen J, Cooper DN, Chuzhanova N, Ferec C, Patrinos GP. 2007. Gene conversion: mechanisms, evolution and human disease. Nat Rev Genet 8:762–775. Chu A, Fill M, Stefani E, Entman ML. 1993. Cytoplasmic Ca2þ does not inhibit the cardiac muscle saroplasmic reticulum ryanodine receptor Ca2þ channel, although Ca2þ‐induced Ca2þ inactivation of Ca2þ release is observed in native vesicles. J Membr Biol 135:49–59. Conti A, Gorza L, Sorrentino V. 1996. Differential distribution of ryanodine receptor type 3 (RyR3) gene product in mammalian skeletal muscles. Biochem J 316:19–23. Counterman B, Ortiz‐Barrientos D, Noor M. 2004. Using comparative genomic data to test for fast‐X evolution. Evolution 58:656–660. Darbandi S, Franck JPC. 2009. A comparative study of ryanodine receptor (RyR) gene expression levels in a basal ray‐finned fish, bichir (polypterus ornatipinnis) and the derived euteleost zebrafish (Danio rerio). Comp Biochem Physiol Biochem Mol Biol 154:443– 448. De Crescenzo V, Fogarty KE, Zhuge R, et al. 2006. Dihydrpyridine receptors and type 1 ryandine receptors constitute the molecular machinery for voltage‐induced Ca2þ release in nerve terminals. J Neurosci 26:7565–7574. George HG, Jundi H, Thomas NL, Fry DL, Lai A. 2007. Ryanodine receptors and ventricular arrhythmias: emerging trends in mutations, mechanisms and therapies. J Mol Cell Cardiol 4:34–50. Gillman LN, Keeling DJ, Howard AR, Wright SD. 2009. Latitude, elevation and the tempo of molecular evolution in mammals. Proc R Soc B 276:3353–3359. Griswold CK. 2011. A model of the physiological basis of a multivariate phenotype that is mediated by Ca2þ signaling and controlled by ryanodine receptor composition. J Theor Biol 282:14–22. Gyorke S, Fill M. 1993. Ryanodine receptor adaptation: control mechanism of Ca2þ‐induced Ca2þ release in heart. Science 60:807– 809. Gyorke S, Velez P, Suarez‐Isla B, Fill M. 1994. Activation of single cardiac and skeletal ryanodine receptor channels by flash photolysis of caged Ca2þ. Biophys J 66:1879–1886. Hakamata Y, Nakai J, Takeshima H, Imoto K. 1992. Primary structure and distribution of a novel ryanodine receptorcalcium release channel from rabbit brain. FEBS Lett 312:229–235. Hamilton SL. 2005. Ryanodine receptors. Cell Calcium 38:253–260. Hansen TF. 2003. Is necessary for evolvability? Remarks on the relationship between pleiotropy and evolvability. Biosystems 69:83– 94. Hedges S, Dudley J, Kumar S. 2006. TimeTree: a public knowledge‐base of divergence times among organisms. Bioinformatics 22:2971– 2972. J. Exp. Zool.

MCKAY AND GRISWOLD Hirata H, Watanabe T, Hatakeyama J, et al. 2007. Zebrafish relatively relaxed mutants have a ryanodine receptor defect, show slow swimming and provide a model of multi‐minicore disease. Development 134:2771–2781. Kimura M. 1980. A simple method for estimating rates of base substitutions through comparative studies of nucleotide divergence. J Mol Evol 16:111–120. Kimura M, Ohta T. 1974. On some principles governing molecular evolution. Proc Natl Acad Sci USA 71:2848–2852. Klein MG, Cheng H, Santana LF, et al. 1996. Two mechanisms of quantized Ca2þ‐release in skeletal‐muscle. Nature 379:455–458. Kushnir A, Mollah A, Wehrens X. 2005. Evolution of the ryanodine receptor gene family. Ryanodine receptors. Dev Cardio Med 254: 1–8. Lanner JT. 2012. Ryanodine receptor physiology and its role in disease. Adv Exp Med Biol 740:217–234. Larkin MA, Blackshields G, Brown NP, et al. 2007. Clustal W and clustal X version 2.0. Bioinformatics 23:2947–2948. Lenhart SE, Mongillo M, Bellinger A, et al. 2008. Leaky Ca2þ release channel/ryanodine receptor 2 causes seizures and sudden cardiac death in mice. J Clin Invest 118:2230–2245. Li P, Chen SRW. 2001. Molecular basis of Ca2þ activation of the mouse cardiac Ca2þ release channel (ryanodine receptor). J Gen Physiol 18:33–44. Liao D. 1999. Concerted evolution: molecular mechanism and biological implications. Am J Hum Genet 64:24–30. Liu Z, Zhang J, Wang R, Wayne Chen SR, Wagenknecht T. 2004. Location of divergent region 2 on the three‐dimensional structure of cardiac muscle ryanodine receptor/calcium release channel. J Mol Biol 338:533–545. Mackrill JJ. 2010. Ryanodine receptor calcium channels and their partners as drug targets. Biochem Pharmacol 79:1535–1543. Martins EP, Housworth EA. 2002. Phylogeny shape and the phylogenetic comparative method. Syst Biol 51:873–880. Meredith RW, Janecka JE, Gatesy J, et al. 2011. Impacts of the cretaceous terrestrial revolution and KPg extinction on mammal diversification. Science 334:521–524. Mori F, Fukaya M, Abe H, Wakabayashi K, Watanabe M. 2000. Developmental changes in expression of the three ryanodine receptor mRNAs in the mouse brain. Neurosci Lett 285:57–60. Morrissette J, Xu L, Nelson A, Meissner G, Block BA. 2000. Characterization of RyR1‐slow, a ryanodine receptor specific to slow‐twitch skeletal muscle. Am J Physiol Regul Integr Comp Physiol 279:1889–1898. Mouton J, Marty I, Villaz M, Feltz A, Mauet Y. 2001. Molecular interaction of the dihydropyridine receptors with type‐1 ryanodine receptors in rat brain. Biochem J 354:597–603. Nielsen R. 2005. Molecular signatures of natural selection. Annu Rev Genet 39:197–218. Oliver T, Garfield D, Manier M, et al. 2010. Whole‐genome positive selection and habitat‐driven evolution in a shallow and a deep‐sea urchin. Genome Biol Evol 2:800–814.

DIVERGENCE AND CORRELATED EVOLUTION OF RyRs Otto S. 2004. Two steps forward, one step back: the pleiotropic effects of favoured alleles. Proc Biol Sci 271:705–714. Paolini C, Fessenden JD, Pessah IN, Franzini‐Armstrong C. 2004. Evidence for conformational coupling between two calcium channels. Proc Natl Acad Sci USA 101:12748–12752. Pavlicev M, Wagner GP. 2012. A model of developmental evolution: selection, pleiotropy and compensation. Trends Ecol Evol 27:316– 322. Percival AL, Williams J, Kenyon JL, et al. 2004. Chicken skeletal muscle ryanodine receptor isoforms: ion channel properties. Biophys J 67:1834–1850. Rice JJ, Jafri MS, Winslow RL. 1999. Modeling gain and gradedness of Ca2þ release in the functional unit of the cardiac diadic space. Biophys J 77:1871–1884. Salathé M, Ackermann M, Bonhoeffer S. 2006. The effect of multifunctionality on the rate of evolution in yeast. Mol Biol Evol 23:721–722. Sasaki S, Nakagaki I, Kondo H. 2002. Involvement of the ryanodine‐ sensitive Ca2þ store in GLP‐1‐induced Ca2þ oscillations in insulin‐ secreting HIT cells. Eur J Physiol 445:342–351. Sobie EA, Dilly KW, dos Santos Cruz J, Lederer J, Jafri MS. 2002. Termination of cardiac Ca2þ sparks: an investigative mathematical model of calcium‐induced calcium release. Biophys J 83:59–78. Sorrentino V, Volpe P. 1993. Ryanodine receptors: how many, where and why? Trends Pharmacol Sci 14:98–103. Sutko JL, Airey JA. 1996. Ryanodine receptor Ca2þ release channels: does divergence in form equal divergence in function? Physiol Rev 76:1027–1071.

163 Tamura K, Nei M. 1993. Estimation of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans and chimpanzees. J Mol Biol Evol 10:512–526. Tamura K, Peterson D, Peterson N, et al. 2011. MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol 28:2731–2739. Weisleder N, Ferrante C, Hirata Y, et al. 2007. Systematic ablation of RyR3 altrs Ca2þ spark signaling in adult skeletal muscle. Cell Calcium 42:548–555. Wolfram Research, Inc. 2008. Mathematica, Version 7.0. Champaign, IL: Wolfram Research, Inc. Wright SD, Howard AR, Keeling DJ, McBride P, Gillman LN. 2011. Thermal energy and the rate of genetic evolution in marine fishes. Evol Ecol 25:525–530. Yamamoto‐Hino M, Sugiyama T, Hikichi K, et al. 1994. Cloning and characterization of human type 2 and type 3 inositol 1,4,5‐ trisphosphate receptors. Receptors Channels 2:9–22. Yang Z. 2007. PAML4: phylogenetic analysis by maximum likelihood. Mol Biol Evol 24:1586–1591. Yano M, Yamamoto T, Ikeda Y, Matsuzaki M. 2006. Mechanisms of disease: ryanodine receptor defects in heart failure and fatal arrhythmia. Nat Clin Pract Cardiovasc Med 3:43–52. Yin CC, D'Cruz LG, Lai FA. 2008. Ryanodine receptor arrays: not just a pretty pattern? Trends Cell Biol 18:149–156.

SUPPORTING INFORMATION Additional Supporting Information may be found in the online version of this article at the publisher's web‐site.

J. Exp. Zool.