Hundreds of Genes Experienced Convergent Shifts in Selective Pressure in Marine Mammals

Hundreds of Genes Experienced Rate Accelerations or Decelerations Associated with the Marine Environment

Using 59 placental mammal genomes, we calculated relative rates of evolution for all branches in amino acid alignments of 18,049 protein-coding gene trees (see Materials and Methods). Briefly, relative rates indicate whether that gene evolved at a faster or slower rate than the expectation after normalization for both the genome-wide average rate for that branch, and the gene-specific rate across all species. In the next step, we employed these relative rates to identify those genes that convergently shifted to higher (or lower) than average rates, specifically in the marine environment as represented by the branches leading to five marine species—bottlenose dolphin, killer whale, walrus, Weddell seal, and West Indian manatee—plus the branch leading to the common ancestor of the bottlenose dolphin and killer whale (fig. 1A). To illustrate this approach, the relative rates of the branches of the GNAT3 gene tree are plotted in figure 2A. GNAT3 encodes a G protein specific to taste and gustatory signaling. The relative rate of the cetacean ancestral branch (the branch leading to dolphins and whales) is positive and large, indicating that this gene accumulated more changes than we would have expected based on the average amount of divergence on that branch for all genes. In representing rates in this way, it is clear that the other marine mammal branches (red) also show higher rates compared to nonmarine branches (blue), which are centered at a relative rate of zero. The higher relative rates on marine branches in GNAT3 compared with all terrestrial branches is highly statistically significant, as measured using a nonparametric test (P = 8.23 × 10 ^− 5, Wilcoxon-rank sum test). The same analytical steps were repeated for all 18,049 genes to identify those with convergently higher or lower relative rates in marine lineages (table 1, supplementary tables S1 and S3, Supplementary Material online). We hereafter refer to genes with higher rates on marine branches as “marine-accelerated” and those with lower rates as “marine-decelerated.” This test does not require all marine branches to have positive relative rates, but because we employ a rank-based test, low P values indicate that a gene has consistently low or high rates in a majority of marine branches when compared with terrestrial branches.

We see in figure 1B that for genes accelerated or decelerated on marine branches, the genome-wide distribution of P values is markedly skewed downward, supporting the hypothesis of convergent shifts in selective pressure on hundreds of genes. In contrast, a control set of mammalian species that was matched for branch length and topology but lacked obvious phenotypic convergence—aardvark, alpaca, camel, microbat, and David’s myotis bat—was only modestly enriched for low P values. The marine and control P value distributions were significantly different for accelerated and decelerated genes (P = 2.46 × 10 ^− 8 and < 2.2 × 10 ^− 16, respectively, Kolmogorov–Smirnov test). Using the distribution of P values in our control species set we apply the q-value method for estimating the fraction of true null hypotheses and infer that genome-wide approximately 560 genes are convergently accelerated and 2,700 genes are decelerated in the marine species. In this analysis we did not designate the pinniped ancestral branch (that leading to walrus and seal) as marine, because the extent to which their last common ancestor was marine is unclear (Arnason et al. 2006; Rybczynski et al. 2009). This was a conservative decision, because if their ancestor were indeed marine, omitting it would only decrease the signal of convergence in our analytical framework. Nevertheless, when the analysis is done with the pinniped ancestral branch as marine, it did not greatly change the results of our study. The gene ranks from our original and the pinniped-ancestor analyses were highly correlated (R = 0.88, P value < 2.2e-16).

Marine-Accelerated Genes Are Enriched for Functions Consistent with Marine Adaptations

To ask which biological functions were enriched among marine-accelerated genes, we searched databases of gene pathways and mutant phenotypes (Smith and Eppig 2009; Liberzon et al. 2011). The convergently marine-accelerated genes contained statistically significant enrichment for functional categories that align well with current knowledge of marine mammal adaptations: sensory systems, muscle function, skin and connective tissue, lung function, and lipid metabolism (fig. 3 and table 1; supplementary tables S1 and S2, Supplementary Material online). The same functional enrichments were observed in the additional analysis with the pinniped ancestor as a marine branch. Marine-decelerated genes on the other hand contained a very different set of enriched functions, broadly described as DNA repair, chromosomal maintenance, immune response and apoptosis (supplementary tables S4 and S5, Supplementary Material online). The slower rate of change in these functions is consistent with increased constraint on somatic cell maintenance as would be required in these relatively long-lived and large-bodied mammals. In fact, this trend is illustrated by the additional large and long-lived species with slower rates in these genes (e.g., double-strand break repair gene XRCC4 is also highly constrained in elephant; figure 2B “XRCC4,” supplementary fig. S2, Supplementary Material online). Generally, the marine-decelerated set encompasses more genes, but their associations are less specific to the marine environment, which is also reflected in their less extreme P values (supplementary table S4 vs. Supplementary Data, Supplementary Material online). We then asked whether the marine-accelerated genes were associated with large body size as well. We identified the enriched functional categories for genes accelerated in a control set of six large-bodied, terrestrial species (alpaca, camel, elephant, horse, and rhino). Only one enriched functional category, “muscle contractility,” was shared between this terrestrial, large-bodied control, and the marine species (supplementary table S7, Supplementary Material online). None of the other significantly enriched functions for marine-accelerated genes were shared with this control.

Positive Selection Drove Marine-Acceleration in Skin and Lung Proteins

Marine-accelerated genes could have resulted from either bursts of adaptive evolution or relaxed constraint in response to the environment. We assigned each gene to these alternative selective regimes using tests of lineage-specific positive selection in CODEML. Specifically, we compared a set of three nested likelihood models: a neutral sites model (M1), the neutral branch-site model (BS1) and the branch-site model allowing positive selection (BS2) (Yang 2007). Rejection of M1 in favor of BS1 indicated relaxed constraint in marine branches, and further rejection of BS1 to prefer BS2 indicated positive selection in the marine species. A majority of the affected functions in marine-accelerated genes were under relaxed constraint, indicating a decline in their importance in the marine environment (table 2; supplementary table S1, Supplementary Material online). However, a high proportion of genes encoding structural components of cartilage and skin showed evidence of positive selection. Prominent examples of positively selected marine-accelerated genes included components of the cell junction-forming desmosome (e.g., PERP), cartilage collagen (COL9A2), and the epidermal cross-linking transglutaminase (TGM3). All of the marine-accelerated genes under positive selection also showed significant evidence of positive selection across the mammalian tree (supplementary table S1, Supplementary Material online), so their positive selection is not specific to marine branches. However, they did nevertheless show increased evolutionary rates on marine branches, suggesting these genes experienced a greater pressure to adapt in the marine environment (supplementary table S3, Supplementary Material online).

Marine-Accelerated Muscle and Sensory Proteins Were under Relaxed Constraint

The remaining marine-accelerated genes represent apparent cases of relaxed constraint, which include relaxation with continued functionality and complete loss of function. It is notable however that some of these genes may also have experienced short episodes of positive selection that are not detected in statistical models but which nonetheless led to rate acceleration (Anisimova et al. 2001). The four muscle-related genes seem to remain functional in most marine species, because two of them did not show genetic lesions in any marine species, β-myosin heavy chain (MYH7) and a muscle-specific glycolysis enzyme (PGAM2) (supplementary table S1 and Supplementary Data, Supplementary Material online alignments in supplementary data). The other muscle genes, nebulin and myoferlin, had early stop codons only in the dolphin sequence. The reason behind the relative rate increases for these muscle genes could be related to the relatively large body size of marine mammals, because a control analysis analyzing convergent rate changes for a set of large-bodied mammals (alpaca, camel, elephant, horse, and rhino) also showed significant rate increases for genes in this category (supplementary table S7, Supplementary Material online).

The strongest case of gene loss was seen in three digestion and absorption genes that, in addition to their marine-acceleration, showed obvious genetic lesions in 60% of marine species on average, compared with 7% in terrestrial species (permutation P value < 0.00001; supplementary table S6, Supplementary Material online). These genes were involved in transport in gastric surface mucus cells (CAPN8), salivary gland aldehyde dehydrogenase activity (ALDH3B2), and intestinal water permeability (aquaporin, AQP10). Many marine-accelerated genes were gustatory and olfactory genes encoding chemoreceptors, G-protein signaling proteins, and neuronal proteins specialized for taste and smell (table 1 and fig. 3; supplementary table S2, Supplementary Material online). Almost all of these olfactory and gustatory marine-accelerated genes showed evidence of relaxed constraint in marine species (supplementary table S1, Supplementary Material online), and a high number of genetic lesions, 29% in marine species versus 12% in terrestrial species (permutation P value < 0.00001). Taste-related genes showed strongly convergent rate increases in that two of the top three marine-accelerated genes were specific to gustatory signaling: G-protein subunit GNAT3 (1) and gustatory ion channel CALHM1 (3) (table 1). This is consistent with recent reports of taste receptor gene loss in toothed and baleen whales (Feng et al. 2014; Zhu et al. 2014). Similarly, we found that GNAT3 contained obvious genetic lesions in both the killer whale and bottlenose dolphin.

Uncharacterized Neuronal Genes Likely Contribute to Sensory Functions

There were also several poorly characterized genes that evolved convergently in marine species. Many of these are in gene families affecting axon guidance (UNC5C, EPHA8, FEZ1, ULK2), synaptic development and neurotransmission (CDH12, GRIN2C, GRIN3B), and other nervous system processes (TMEM25, OTOF, GLRA4, CLSTN2), yet they are not currently ascribed to organism-level neurological processes (supplementary table S1, Supplementary Material online). Given the strong signal of convergent relaxation of constraint in marine sensory systems (supplementary tables S1 and S6, Supplementary Material online), we propose that these genes function in olfactory and/or gustatory neurons in terrestrial mammals. These genes could also be involved in other sensory modalities as we observe acceleration in genes known to affect temperature (TRPV3) and pain sensation (ASIC4 and TACR2).

Alternatively, two of the marine lineages (cetaceans and pinnipeds) are vocal learners and convergence acceleration in neuronal genes may represent a shift in selection pressure due to that trait (Petkov and Jarvis 2012; Janik 2014; Jarvis et al. 2014; Reichmuth and Casey 2014). Although it is difficult to advance this hypothesis given that we do not know whether the ancestors of cetaceans or pinnipeds were vocal learners.

Marine-Accelerated Branches Do Not Contain an Excess of Site-Specific Convergent Amino Acid Changes

We tested the hypothesis that convergent or parallel amino acid changes at specific sites contributed to the reported marine-associated accelerations. We used the branch model of CODEML to reconstruct the substitutions on all branches in a mammalian phylogeny containing up to 49 placental mammals. We then tallied all convergent amino acid substitutions on marine branches for the top 45 marine-accelerated genes, including both strictly convergent and parallel changes. Strictly convergent refers to a scenario where the ancestors of two marine species had different amino acids at a given site that changed to the same amino acid on both marine branches. Parallel changes refer to two marine ancestors having the same amino acid at a given site which evolved to the same, new amino acid in both extant species. The proportion of all changes that was convergent on marine branches was generally low, with 42% of branches having none (supplementary fig. S3A, Supplementary Material online). Convergent changes on marine branches could have resulted from neutral processes or alternatively due to positive selection for the same amino acid variant. For the latter case, we would expect there to be more convergent amino acid changes in these genes for marine species compared to negative control species. We tallied convergent changes within a control set of species not expected to show convergence and chosen to match the topology and branch lengths of the marine species (supplementary fig. S3C, Supplementary Material online) (alpaca, camel and their ancestral branch, bushbaby, human, and aardvark). As an additional control, we tallied the convergent changes between the marine and control branches (supplementary fig. S3D, Supplementary Material online). The marine branch data set did not show an excess of convergent amino acid sites (mean proportion = 0.086) compared to the control data sets (0.088 and 0.052, for control branches and marine- vs. -control branches, respectively) (supplementary fig. S3, Supplementary Material online). Furthermore, the proportion of convergent changes on marine branches is small compared with the amount of excess changes that led to the acceleration in relative rate (mean proportion = 0.549; supplementary fig. S3B, Supplementary Material online). Overall, we found no evidence for adaptively driven site- and amino acid-specific convergence in marine-accelerated genes.

Comparison of Approaches to Identify Convergent Evolution

A previous study of molecular convergence in marine mammals employed the hypothesis that convergence would be seen as single amino acid sites changing to the same amino acid. They did not find an excess of such changes over background (Foote et al. 2015). We also asked if our marine-accelerated genes experienced such amino acid site convergence. We did not find an excess of convergent changes on marine branches when compared to negative control branches (supplementary fig. S3, Supplementary Material online). It has been reported that proteome-wide searches for convergent amino acids should have low power and that observed convergent changes can be sufficiently explained by neutral models that do not rely on adaptively driven convergence (Zhang and Kumar 1997; Thomas and Hahn 2015; Zou and Zhang 2015a). As an alternative approach, we introduced a novel analytic framework that focuses on those genes showing convergent acceleration or deceleration in their whole-gene evolutionary rates. This strategy represents a more course-grained perspective to identify common sets of genes and pathways that experienced a shift in their evolutionary pressures during phenotypic convergence. Our application of this gene-centric strategy to marine mammals uncovered a strong statistical excess of convergent shifts in evolutionary rate (fig. 1B). This suggests that convergent changes in selection at the gene level is relatively common compared to convergence at single amino acids, and that its signal is more readily detected. These two approaches operated at distinct levels but employed the same basic data set. Accordingly, we asked if they detected a shared set of genes. We compared our top 690 marine-accelerated genes with nominal P values below 0.05 to the 15 genes reported by Foote et al. in their table 1 as being both under positive selection in marine lineages and containing convergent amino acid changes. There was an overlap of two muscle-related genes, MYH7B and SMPX. While an overlap of two genes (13%) is at higher frequency than expected (3.8%), it was not entirely unlikely by chance (P = 0.11, binomial test). However, the comparison suffered from low power due to a low number of genes identified by Foote et al., so it is still plausible that the two approaches could be recovering similar targets in their quest for convergent evolution.

Foote et al. also reported genes under positive selection in one or multiple marine branches using the branch-site models of CODEML (their supplementary tables S3 and S9–S11, Supplementary Material online) (Yang 2007; Foote et al. 2015). However, the tests in Foote et al. did not exclude the presence of positive selection in terrestrial species. To illustrate, the gene ANPEP, which was featured in Foote et al., yields strong evidence for positive selection across all branches in the mammalian tree (CODEML sites models M1 vs. M2, P = 1.8 × 10⁻²¹). When all marine species are removed, ANPEP still shows strong evidence of positive selection, indicating that positive selection in ANPEP is not restricted to marine species (M1 vs. M2, P = 2.7 × 10⁻¹⁰; supplementary table S8, Supplementary Material online). Four of the five genes that passed Foote et al.’s correction for multiple testing (ANPEP, CD48, EMP1, and MUC1) show the same evidence of positive selection outside of marine species (supplementary table S8, Supplementary Material online). The exception was ZNF582, for which our alignment showed no evidence of positive selection in any model comparison, including the marine-species branch-site test. Our method, on the other hand, uses a different criterion that specifically tests for convergent acceleration or deceleration in marine branches when compared with terrestrial branches. The resulting gene set contains genes that convergently experienced positive selection or relaxed constraint (marine-accelerated) or additional constraint (marine-decelerated) in the marine environment. For these reasons, our studies report very different genes.

There are clear advantages and limitations to our gene-centric strategy. Generally, we found that a search for convergent rates in protein-coding genes can excel at revealing adaptations in physiological and structural genes. However, genes controlling convergent morphological adaptations were not identified. We expect this limitation is due to the fact that morphological changes typically involve regulatory rather than protein-coding changes (Carroll 2008); however, this suggests a potential extension to study rate changes in cis-regulatory regions as well. A notable advantage to our method is that it reveals convergently evolving biological processes even when nothing is known about adaptation a priori. Our approach applied no prior hypotheses about which genes or functions would be marine-accelerated or decelerated, yet it recovered genes significantly enriched for functions that are altered in marine mammals (fig. 3). Furthermore, it may be useful to assign probable functions to poorly characterized but convergently evolving genes, as we demonstrated for 11 neuronal genes likely participating in sensory perception. Overall, searching for convergent shifts in evolutionary rate represents a promising strategy to identify candidate genes underlying convergent phenotypic traits, as demonstrated in this and past studies (Goodman et al. 2009). The study of marine mammals is just the first application of our relative evolutionary rates approach. It can be readily applied to essentially any environmental or phenotypic variable with sufficient range in the mammalian phylogeny, or other taxonomic group for that matter.

Calculating Correlations of Protein-Coding Genes with the Marine Environment

Amino acid alignments of orthologous genes were derived from the “100-way” 100 vertebrate species genomic multiz alignment available at the UCSC genome browser (Blanchette et al. 2004), we studied those alignments with a minimum of 30 species. For each amino acid alignment we estimated branch lengths using the AAML program from the phylogenetic analysis using maximum likelihood (PAML) package (Yang 2007). Branch lengths were estimated under an empirical model of amino acid substitution rates with rate variability between sites modeled as a gamma distribution approximated with four discrete classes (for computational efficiency) and an additional class for invariable sites (AAML model “Empirical + F”) (Yang 1996; Whelan and Goldman 2001). Branch lengths were estimated by AAML on a published mammalian species tree topology (Murphy et al. 2004). Raw branch lengths were transformed into relative rates using a projection operator method (Sato et al. 2005). This method begins with each gene’s branch lengths in the form of a vector that we then scale to a vector length of 1. The next step is to transform these branch lengths in to relative rates. Relative rates quantify how much faster or slower this gene changed on a given branch after factoring out the divergence on that branch resulting from parameters affecting all genes (e.g., the time since speciation, effective population size, mutation rate). Removing these genome-wide effects will leave only the gene-specific relative rate for a given branch. Relative rates are calculated using a normalization vector that is calculated as the average amount of divergence on each branch across all scaled genes. Finally, we determine the relative rate of each branch as the residual of that branch after factoring out the normalization vector. A relative rate vector is thus defined as the component of the branch-length vector that is orthogonal to the normalization vector. These branch-specific relative rates were then used to perform a Wilcoxon rank sum test and correlation analysis over the binary variable of “terrestrial” or “marine” branches (fig. 1A). Marine branches are those leading to the manatee (Trichechus manatus latirostris), Weddell seal (Leptonychotes weddellii), walrus (Odobenus rosmarus divergens), killer whale (Orcinus orca), bottlenose dolphin (Tursiops truncatus), and the cetacean ancestral branch leading to the killer whale and dolphin. Control branches for P value enrichments in figure 1 were those leading to the aardvark (Orycteropus afer afer), alpaca (Vicugna pacos), camel (Camelus ferus), microbat (Myotis lucifugus), and David’s myotis bat (Myotis davidii), as well as the ancestral branch leading to alpaca and camel. Large-body control branches were those leading to alpaca (V.pacos), camel (C.ferus), elephant (Loxodonta africana), horse (Equus caballus), and white rhino (Ceratotherium simum).

Functional Enrichment Analysis

Functional information for marine-accelerated and marine-decelerated genes in table 1 was taken from the Uniprot and RefSeq databases, and from literature cited directly (Pruitt et al. 2007; Uniprot Consortium 2007). Computational tests for functional enrichment were performed using the hypergeometric test with the background set of genes restricted to genes that were tested for marine convergence and had at least one annotation in the corresponding annotation file. Correction for multiple testing was performed using false discovery rate (FDR) q-values (Storey 2002). We used two sources of annotations for computational tests for functional enrichment, the “canonical pathways” from MSigDB (Liberzon et al. 2011) and mammalian phenotypes from MGI (Smith and Eppig 2009). The mammalian phenotype annotations were compiled by associating gene symbols listed in the genotype name to the reported phenotypes and all their ancestors in the mammalian phenotype ontology (https://bioportal.bioontology.org/ontologies/MP, last accessed June 6, 2016). Enriched categories are displayed in a network in figure 3. In this network we required a minimum of 2-fold enrichment within the 500 most marine-accelerated genes and an FDR q-value below 5%. Edges were drawn between categories if they shared a minimum of 20% of their accelerated gene content as calculated using the smaller of the two categories. Accordingly, a value of 100% results when all of a category’s accelerated genes are found in the other category. The network in figure 3 was visually laid out by an edge-weighted force-directed algorithm in Cytoscape (Killcoyne et al. 2009).

FDR Analysis

FDR analysis was performed on probabilities resulting from the Wilcoxon rank sum test. FDRs and the proportion of tests conforming to the null hypothesis were estimated by a modified version of the “qvalue” function in the corresponding “qvalue” R module (Storey and Tibshirani 2003; Storey 2015). The q-value method first estimates π₀ (the proportion of tests for which the null hypothesis holds) by comparing the cumulative distribution of observed P values to a uniform distribution. The π₀ value is subsequently used to calculate a FDR for each P value threshold. We modified the π₀ calculation by replacing expected uniform quantiles with the observed null distribution quantiles calculated with the control species set. This approach parallels the empirical permutation-based FDR calculations that are often employed in genomic data analysis which are typically more conservative as they preserve dependence in the data structure. In our case this produced FDR estimates that are slightly more conservative because the P value distribution from our control species was not perfectly uniform. The slight departure from uniformity in the control species was seen as a slight enrichment of small P values due to nonindependence arising from phylogenic relationships.

Phylogenetic Models of Selective Pressure

Of the marine-accelerated genes, the top 115 were subjected to phylogenetic models of codon evolution in order to test for significant evidence of relaxation of constraint or positive selection over the marine mammal branches. We used CODEML of the PAML package to run three models: the branch-site neutral model (BS1), the branch-site selection model (BS2), and a neutral sites model (M1) (Yang 2007). Likelihood ratio tests (LRT) between BS1 and its nested null model M1 were used to assign significance of relaxation of constraint on marine mammal branches specifically. LRTs between BS2 and its null BS1 were used to infer positive selection on the marine mammal branches (Zhang et al. 2005). Probabilities were assigned for each of these two LRTs using the chi-square distribution with 1 degree of freedom. In a similar way, mammal-wide positive selection was inferred using the M8 (selection) and M8A (neutral) sites models and their respective LRT (Swanson et al. 2004). Significance for this LRT was assigned with the chi-square distribution with 1 degree of freedom, which is a conservative assumption (Wong et al. 2004).