1. Computational and Systems Biology
  2. Genetics and Genomics
Download icon

Pan-mammalian analysis of molecular constraints underlying extended lifespan

  1. Amanda Kowalczyk  Is a corresponding author
  2. Raghavendran Partha
  3. Nathan L Clark
  4. Maria Chikina  Is a corresponding author
  1. Joint Carnegie Mellon University-University of Pittsburgh PhD Program in Computational Biology, United States
  2. University of Pittsburgh, United States
  3. University of Utah, United States
Research Article
  • Cited 0
  • Views 983
  • Annotations
Cite this article as: eLife 2020;9:e51089 doi: 10.7554/eLife.51089

Abstract

Although lifespan in mammals varies over 100-fold, the precise evolutionary mechanisms underlying variation in longevity remain unknown. Species-specific genetic changes have been observed in long-lived species including the naked mole-rat, bats, and the bowhead whale, but these adaptations do not generalize to other mammals. We present a novel method to identify associations between rates of protein evolution and continuous phenotypes across the entire mammalian phylogeny. Unlike previous analyses that focused on individual species, we treat absolute and relative longevity as quantitative traits and demonstrate that these lifespan traits affect the evolutionary constraint on hundreds of genes. Specifically, we find that genes related to cell cycle, DNA repair, cell death, the IGF1 pathway, and immunity are under increased evolutionary constraint in large and long-lived mammals. For mammals exceptionally long-lived for their body size, we find increased constraint in inflammation, DNA repair, and NFKB-related pathways. Strikingly, these pathways have considerable overlap with those that have been previously reported to have potentially adaptive changes in single-species studies, and thus would be expected to show decreased constraint in our analysis. This unexpected finding of increased constraint in many longevity-associated pathways underscores the power of our quantitative approach to detect patterns that generalize across the mammalian phylogeny.

Introduction

Humans age in the sense that an individual’s probability of dying increases as a function of time lived. Interestingly, this trend is not true of all species but is true of mammals generally (Jones et al., 2014). Numerous hypotheses of mammal-specific aging exist, including the antagonistic pleiotropy hypothesis (Williams, 1957) and the mutation accumulation hypothesis (Medawar, 1952), both of which refer to changes within an individual throughout its lifetime that result in aging. The antagonistic pleiotropy hypothesis postulates that genes that are beneficial to early life become detrimental later in life. Such genes are retained because they increase early-life reproductive output and thus increase fitness. The mutation accumulation hypothesis predicts that a gradual accumulation of errors in DNA sequence as a result of repeated replication during a lifetime’s worth of cell divisions will lead to a gradual breakdown of functionality. Support for both hypotheses has been found in individual species. Recent work has identified disease-related SNPs in age-related genes that are beneficial in early life and detrimental in later life in humans, thus indicating selective pressures associated with gene evolution related to aging and supporting the antagonistic pleiotropy hypothesis (Rodríguez et al., 2017). A study of SNVs in the human brain found that number of mutations was positively correlated with age and that mutations were at loci associated with age-related disease, thus supporting the mutation accumulation hypothesis (Lodato et al., 2018). However, in human populations, variability in the aging phenotype is limited and many confounding biological changes are correlated with aging, making it difficult to pinpoint specific biological processes that are causal and thus amenable to manipulation.

On the other hand, lifespan varies dramatically (>100 fold) across mammals (Nowak, 1999), making comparative genomics a fruitful avenue for aging research. Numerous studies have investigated the genomic features of mammals with extreme lifespan such as bats (Foley et al., 2018; Seim et al., 2013), naked mole-rats (Kim et al., 2011), whales (Keane et al., 2015), and elephants (Sulak et al., 2016) to identify potential causative genetic changes. In Myotis, the longest-lived bat genus, species show lack of telomere shortening and corresponding expression changes in telomere maintenance and DNA repair genes (Foley et al., 2018). Comparative genomics studies have also suggested that changes in the insulin growth factor one pathway may enable increased lifespan and cancer resistance in Myotis brandtii (Brandt’s bat) (Seim et al., 2013). Similarly, naked mole-rats show differential regulation of genes associated with macromolecule degradation, mitochondrial function, and TERT, a gene associated with telomere maintenance, as well as changes to genes related to tumor suppression (Kim et al., 2011). In primates, genes associated with cardiovascular function, coagulation, and healing have been demonstrated to show evolutionary correlations with lifespan (Muntané et al., 2018). Sequencing of the bowhead whale genome revealed species-specific changes in DNA repair, cell cycle, and aging genes (Keane et al., 2015). In elephants, a striking increase in TP53 copy number has been linked to increased cancer resistance enabling longer lifespan (Sulak et al., 2016). Despite compelling results, these single- and limited-species studies have limitations. The species studied differ from their nearest sequenced relatives in multiple physiological traits as well as millions of nucleotides. Thus, while single-species studies have yielded some credible candidates for genes associated with increased lifespan, it is difficult to know to what extent these represent insights into the universal mechanisms of lifespan regulation rather than species-specific adaptation or coincidental neutral changes. In this study, we develop new methodology to evaluate the relationship between the evolutionary constraint of genes and pathways and quantitative lifespan traits in an unbiased, genome-wide, pan-mammalian analysis.

The wide range of lifespans across the mammalian phylogeny (Figure 1A) provides the ideal dataset to investigate lifespan from a comparative genomics perspective. Because independent changes in lifespan occurred repeatedly in the mammalian species tree, lifespan can be viewed as a convergent trait. Molecular features that correlate with convergent changes in lifespan therefore may also occur repeatedly across a variety of organisms. In our study we use protein evolutionary rates quantified as the number of amino acid substitutions on a phylogenetic branch to infer convergent rate shifts associated with lifespan traits across the mammalian phylogeny.

Figure 1 with 11 supplements see all
Many genes have evolutionary rates correlated with longevity phenotypes as demonstrated by analysis with RERconverge.

(A) A subset of species used for this analysis alongside their maximum longevity values. Lifespan varies widely across mammals independent of phylogeny. (B) Mammal body size and maximum lifespan values for 61 species. Lines represent the 3L phenotype and the ELL phenotype (also see Figure 1—figure supplement 1). (C) RERconverge pipeline to find correlation between relative evolutionary rates of genes and change in lifespan phenotypes. (D and E) Distribution of p-values from correlations between evolutionary rates of genes and change in the 3L and ELL phenotypes indicate an enrichment of significant correlations (also see Figure 1—figure supplement 11).

Evolutionary rates are useful for linking phenotypes to genes because they reflect evolutionary constraint experienced by a genetic element (Zhang and Yang, 2015). In the absence of diversifying selection, genetic elements that support a specific trait are expected to be more constrained in species where the trait has a larger contribution to fitness. In agreement with this expectation, multiple studies have shown that a genetic element providing a function less important for a given species is under less constraint and hence exhibits a faster evolutionary rate (Clark et al., 2013; Janiak et al., 2018; Roscito et al., 2018; Wertheim et al., 2015). Reciprocally, when an element becomes relatively more important, its rate is expected to slow. Thus, in cases of phenotypic convergence, rates can be exploited to reveal important genes associated with the phenotype, such as changes to muscle and skin genes associated with the mammalian transition to a marine environment (Chikina et al., 2016) and loss of constraint of vision-related genetic elements in subterranean mammals (Partha et al., 2017; Prudent et al., 2016). Rate shifts can thus provide an evolutionary perspective on the contribution of genes, non-coding elements, or pathways to phenotypes of interest (Hiller et al., 2012). Here we report the genome-wide, pan-mammalian correlations between evolutionary rates of genes and lifespan phenotypes.

Results

In mammals, lifespan is strongly positively correlated with adult body size such that the largest mammals (whales) are longest-lived and the smallest mammals (small rodents) are shortest lived (Figure 1B). However, if lifespan is corrected for body size, species including bats, the naked mole-rat, and some primates are clearly exceptionally long-lived given their body sizes. The term ‘longevity’ is applied to both phenotypes and previous studies have focused on small numbers of species at both phenotypic extremes to find genetic and physiological explanations for the extended phenotypes (Foley et al., 2018; Keane et al., 2015; Kim et al., 2011; Seim et al., 2013; Sulak et al., 2016). In our study we explicitly distinguish two different extended longevity traits: the ‘long-lived large-bodied’ trait (3L) and the ‘exceptionally long-lived given body size’ trait (ELL). Using maximum lifespan and body size data (Tacutu et al., 2018) we define the 3L and ELL phenotypes (Supplementary file 1) to be the first and second principal components of body size and maximum lifespan (Figure 1B, Figure 1—figure supplement 1). The resulting trait values are orthogonal with respect to each other, ensuring that they can be analyzed independently.

Having defined the 3L and ELL phenotypes, we compute the association between these phenotypes and protein-specific relative evolutionary rates (RERs) using the RERconverge method (Kowalczyk et al., 2018; Partha et al., 2018) (Figure 1C and Materials and methods). Relative evolutionary rates quantify the deviation in evolutionary rate of a protein along a specific phylogenetic branch from proteome-wide expectations. Negative RERs indicate fewer substitutions than expected due to increased constraint. Positive RERs correspond to more substitutions than expected, which could arise due to relaxation of constraint or positive selection.

After computing correlations between all protein relative evolutionary rates and the 3L and ELL phenotypes, we find an excess of low p-values (Figure 1D and Figure 1E). In order to evaluate how much of the signal is due to true association with the phenotypes, we used a phylogenetically restricted permutation strategy (termed ‘permulations’, see Materials and methods) to generate null correlation statistics. We find that the real and permuted p-value distributions are indeed different, indicating that a large fraction of genes are correlated with lifespan phenotypes. We quantify the fraction of non-null genes using the π1 method (Storey, 2003). Using our permulation p-values as the null distribution, the fraction of true positives was inferred to be ~15% (π1 = 0.153) for the 3L phenotype and ~7.5% (π1 = 0.075) for the ELL phenotype. Using the theoretical uniform null distribution, the corresponding values are π1=0.108 and π1 = 0.021, respectively. Regardless of the null distribution choice, our analysis clearly demonstrates a significant molecular signal for gene evolutionary rates correlated with lifespan phenotypes, with a considerably higher number of associations for the 3L phenotype.

Our analysis investigates both positive and negative correlations between evolutionary rates of genes and changes in lifespan phenotypes. Positive correlations represent genes with faster evolutionary rates in species with high 3L and ELL phenotype values relative to species with low phenotype values. Conversely, negative correlations represent genes with slower evolutionary rates in species with high 3L and ELL phenotype values relative to species with low phenotype values. Note that correlation directionality is relative – faster evolution in low phenotype values corresponds to slower evolution in high phenotype values and vice versa. However, the choice has important consequences for biological interpretation and necessitates an in-depth discussion.

RERs reflect the amount of constraint on the genetic element. A decrease in RER implies greater evolutionary constraint and a greater contribution of that gene to fitness. An increase in RER, on the other hand, has two nearly opposite interpretations. An increased RER could arise because of a relaxation of evolutionary constraint driven by a reduced contribution to fitness. Alternatively, an increased RER could arise due to positive (also termed directional) selection, which implies that that the gene is actively undergoing directed evolution and thus could be contributing to trait-related innovation. However, since there is no default evolutionary rate for protein coding sequence, genes evolving slower in long-lived species could just as easily be interpreted as evolving faster in short-lived species. In our study, we have chosen to interpret the rate changes with reference to the effect in long-lived species (Figure 2). The ancestral mammal is believed to have been small and short-lived, and thus large values of 3L and ELL are derived traits. Consequently, our interpretation reduces to assuming that change in phenotype and shifts in evolutionary rate coincide.

During selective pressure shifts that drive phenotypic changes, the genetic evolutionary rate landscape shifts for many genes.

Imagine a simplified case where Gene 1 and Gene 2 start at the same evolutionary rate (solid red and blue lines, respectively). Trait-driving genes like Gene 1 enter a transient state of increased evolutionary rates to drive trait change. Complementary genes that support the trait like Gene 2 experience increased purifying selection to allow the trait to persist as it is established. True rates (solid red and blue lines) are not measurable because they represent rates at inaccessible ancestral evolutionary timepoints – only accessible extant sequences can be used to access average rates over time (dashed red and blue lines, which represent positive and negative RERs). Therefore, transient periods of positive selection are less readily able to be captured by RERconverge than sustained purifying selection due to the greater magnitude of their captured rate shifts.

With this interpretation, negative correlations imply that a gene is more important in species with large values of 3L and ELL, while positive correlations imply either a relaxation of constraint or positive selection in species with large values of 3L and ELL. While such positively correlated genes would be of great interest if they indeed represent molecular innovations underlying evolution of extended lifespan, we find relatively fewer such genes, no evidence of positive selection among them (see Supplementary file 6), and no enriched pathways associated with them.

We thus focus our analysis on the negatively correlated genes, which we interpret as being under increased purifying selection in species with high longevity (3L or ELL) values. It is theoretically possible for negative correlations to be caused by accelerated evolution in small and short-lived species because of directional selection associated with development of low longevity values. However, this possibility can be rules out via branch-site models for positive selection using the low longevity species as foreground, as these show little evidence for positive selection (see Supplementary file 6). Together, these analyses support a single interpretation of the main rate convergence signal as a decrease in evolutionary rate, and thus an increase in purifying selection, experienced by species with large 3L and ELL values. It is important to emphasize that these genes are unlikely to have contributed to molecular innovation that lead to the establishment of 3L and ELL traits, but rather these represent existing biological systems that become especially important after the traits are established (see Discussion).

While we observe a clear excess of genes at low p-values, we focus on pathway enrichment analysis which both demonstrates a stronger signal and facilitates interpreting our results in the context of existing knowledge. We investigate enriched pathways for both 3L and ELL phenotypes using a rank-based method (Supplementary file 2). After performing standard multiple-hypothesis testing corrections on the empirical p-values from permulations, there remains considerable pathway-level signal underlying the 3L and ELL traits. Both the gene-level and pathway-level results were highly robust to species removal, which indicates the biological pathways revealed here are important for longevity across mammals and are not restricted to specific species (see Materials and methods, Figure 1—figure supplement 2, Supplementary file 4, and Supplementary file 5).

Among pathways under increased constraint in 3L species, we find a striking abundance of pathways related to cancer control. Those pathways can be organized into the broad categories of ‘cell cycle control’, ‘cell death’, and ‘innate and adaptive immunity’, and they also include other cancer-related pathways such as p53 regulation and telomere maintenance (Figure 3A). We likewise see a significant enrichment in cancer-related genes more broadly. We compared correlation statistics between known cancer genes from Bailey et al. (2018) and found that tumor suppressor genes had significantly lower correlation statistics (rho) than all genes (Wilcoxon rank-sum p-value=4.114e-8) while oncogenes had no significant difference in correlation statistics compared to all genes (Wilcoxon rank-sum p-value=0.3745) (Figure 3—figure supplement 1). This may indicate preferential purifying selection on tumor suppressor genes in large, long-lived species for lower cancer incidence. Considering that these findings, our 3L results can be naturally interpreted in the context of Peto’s Paradox (Peto, 2016). The paradox reasons as follows: if all cells have a similar probability of undergoing a malignant transformation, organisms with more cells should have a greater risk of developing cancer. However, empirical cancer rates do not vary with body size (Peto, 2016), which implies that larger animals harbor mechanisms to suppress cancer rates. Top 3L constrained pathways are associated with multiple cancer control mechanisms, including DNA repair, cell cycle control, cell death, and immune function (Figure 3B and Figure 3C). A normal cell’s transformation to malignancy involves failure of all these processes, and our analysis suggests that 3L animals are invested in the maintenance of each of their associated pathways through increased purifying selection. Based on enrichment and permulation results, we can infer that cell cycle fidelity, an early step in cancer development, is most important over evolutionary time scales for 3L species. Further, there is no evidence for enrichment of pathways associated with metastasis and angiogenesis, later steps in cancer development. This finding suggests that large, long-lived species have experienced increased selective pressure to protect pathways involved in early cancer stages but not later stages, perhaps because the most severely negative fitness impacts of cancer are felt earlier in its development. Species-specific cancer control mechanisms have been identified in individual species, such as increased TP53 copy number in elephants (Sulak et al., 2016), but we show here that investment is cancer control is key to longevity across the entire mammalian phylogeny because top enriched pathways for the 3L phenotype do not depend on a handful of species (Figure 1—figure supplement 3 and Figure 1—figure supplement 2).

Figure 3 with 1 supplement see all
Pathways that evolve slower in long-lived, large-bodied mammals are related to control of cancer.

(A) Significantly enriched pathways under increased constraint in species with larger values of the 3L phenotype. Each dot represents a pathway, and the size and color of the dot represents the negative log of the rank-sum enrichment statistic. Width of lines connecting pathways represent the number of genes the pathways have in common. (B and C) Pathways under increased constraint in 3L species play various roles in cancer control. Pathways associated with early stages of cancer development (DNA repair, cell cycle control, cell death, and immune functions) are significantly enriched, while pathways for later stages of cancer development (angiogenesis and metastasis) are not enriched. In C), each dot represents a pathway. Yellow dots have significant permulation p-values while black dots do not. Note that dots for ‘All’ pathways are excluded for the sake of clarity.

An additional pathway that shows a strong signal of increased constraint with the 3L phenotype is the insulin-like growth factor (IGF1) signaling pathway (Figure 4), which deserves special consideration. Perturbations of IGF1 signaling result in changes in lifespan and body size in diverse organisms (Johnson et al., 2013; Kimura et al., 1997; Stout et al., 2013), which suggests that the IGF1 pathway may be a source of innovation underlying the evolution of the 3L trait. However, we find that IGF1 pathway genes in fact evolve more slowly and are thus under increased purifying selection in large, long-lived species. The magnitude of this signal is quite striking; IGF1 and IGF1R are ranked 1014 and 94 respectively for increased constraint and several other pathway members are near the top (see Figure 3A). A possible explanation is that the IGF1 pathway plays an important role in cancer control (Larsson et al., 2005), and this cancer-related signal of constraint dominates any adaptive signal related to lifespan. There are also further reasons to believe that the IGF1 pathway is not the main source of the 3L trait. Across the mammalian phylogeny, lifespan is strongly correlated with body size, but genetic perturbations in the IGF pathway result in longer-lived individuals that are of smaller size (Holzenberger et al., 2003; Sutter et al., 2007) thus decoupling the two traits. This strongly suggests that changes in the IGF1 pathway are unlikely to drive the natural evolution of the large, long-lived phenotype, which is established through a different, yet-unknown mechanism.

The IGF1 signaling pathway is under increased evolutionary constraint in long-lived, large-bodied mammals.

(A) IGF1 signaling pathway genes are significantly correlated with change in the 3L phenotype. (B) The IGF1 signaling pathway is significantly enriched for increased evolutionary constraint in large, long-lived species. The barcode indicates ranks of genes in the pathway within the list of all pathway-annotated genes. The worm indicates enrichment as calculated by a tricube moving average, a type of moving average in which values near the end of the sliding window are down weighted to reduce the effect of extreme values in any given window. The dashed horizontal line marks the null value indicating no enrichment. (C) The IGF1 signaling pathway contains many genes whose evolutionary rates are negatively correlated with the 3L phenotype. Shading indicates the Rho-signed negative log p-value for the correlation. Genes in white are not included in the IGF1 pathway annotation used to calculate pathway enrichment statistics, but they are included in the diagram for sake of completeness. The GNB2L1 gene (gray) is in the IGF1 pathway annotation, but correlation statistics were not calculated for that gene because too few branches in the gene tree met the minimum branch length cut-off.

For the ELL phenotype, we find a smaller, more focused set of enriched genes and pathways. Of those significantly enriched constrained pathways for ELL, we see some overlap with functional groups represented in results from the 3L phenotype, notably immune-related and DNA repair pathways (Figure 5A). However, although the functional groups are the same, the pathways contained within them differ between the two phenotypes (File S2, File S3). In particular, the only significantly constrained DNA repair pathways for the 3L phenotype involve Fanconi’s anemia, while the ELL phenotype shows significantly constrained DNA repair pathways for a variety of repair functions (Figure 5). Such pathways stand out not only because of the connection between DNA repair and cancer control, but also because of the observed relationship between DNA repair and aging independent of cancer incidence. This relationship can be demonstrated experimentally by creating double-stranded DNA breaks in laboratory mice to induce an aging phenotype (White et al., 2015). There is also evidence that DNA damage causes dysregulation of the cellular chromatin state and thus can contribute to aging even in post-mitotic cells (Oberdoerffer et al., 2008; Shi and Oberdoerffer, 2012).

DNA repair pathways are under increased evolutionary constraint in mammals that are exceptionally long-lived given their body size.

(A) Significantly enriched pathways under increased constraint in species with larger values of the ELL phenotype. Bar height indicates the negative log permulation p-value for each pathway, and the color of bars indicates the pathway enrichment statistic. (B) DNA repair pathways are more significantly enriched for increased evolutionary constraint in species with large values of the ELL phenotype than species with large values of the 3L phenotype. The barcodes indicate ranks of genes in the pathways within the list of all pathway-annotated genes. The worms indicate enrichment as calculated by a tricube moving average, a type of moving average in which values near the end of the sliding window are down weighted to reduce the effect of extreme values in any given window. The dashed horizontal lines mark the null value indicating no enrichment.

In addition to DNA repair-related pathways, we also noted pathways related to NFKB signaling for which overexpression in downstream targets has been associated with aging. Experimental evidence suggests a connection between NFKB signaling and DNA repair through sirtuins, a chromatin regulator family that has already been implicated in lifespan control (Howitz et al., 2003; Mao et al., 2012). Sirtuins mediate DNA damage-induced dysregulation and are also responsible for silencing NFKB-regulated genes, thus connecting the two processes (Salminen et al., 2008). Overall, our analysis strongly suggests that fidelity in DNA repair and NFKB signaling contributes to the fitness of ELL species, indicating that these pathways may be a fruitful avenue for aging research and intervention.

Discussion

We employed an evolutionary rates-based method as an unbiased, genome-wide, pan-mammalian scan to identify genes and pathways that evolve significantly slower in long-lived large-bodied species (3L) and species that are exceptionally long-lived given their body size (ELL). Pathways related to cancer control, including cell cycle, DNA repair, cell death, and immunity, evolve significantly slower in 3L species, which suggests that cancer resistance is an important functionality the enable to evolution of large and long-lived species. Alternatively, a broader set of DNA repair genes and a more focused set of immune genes related to NFKB signaling evolve significantly slower in ELL species, both of which may be linked to effective DNA repair in order to preserve chromatin state.

Our analysis differs from previous efforts in both methodology and results. Firstly, we do not consider individual amino acid changes within a protein as the unit of convergence, but rather we calculate the overall evolutionary rate of a protein on each branch of a phylogeny and consider these evolutionary rates as the unit of convergence. Secondly, unlike previous studies that focused on lineage- or species-specific changes, we look for correlations between evolutionary rates and quantitative life history traits across the entire phylogeny. This pan-mammalian approach allows us to generalize our findings to describe evolutionary trends throughout all mammals. We also draw a careful distinction between absolute lifespan and relative lifespan, which allows us to make district conclusions about the 3L and the ELL phenotypes.

The most important distinction between our work and previous work is that prior studies began with the assumption that some genes must be under positive selection in association with evolution of extended lifespan. From our unbiased analysis, we find that in fact the strongest signal of evolutionary convergence is that of increased constraint on certain genes and gene families in long-lived species (see Figure 3—figure supplement 1). While some of the pathways have been reported previously (such as cell cycle control, DNA repair, telomerase repair, and IGF1 signaling) our finding is actually the opposite of (though not necessarily contradictory to) the positive selection signal that has been emphasized in previous work, which from the perspective of evolutionary rate is decreased constraint. In fact, we find comparatively few genes and no significant pathway enrichments for the opposite trend – faster relative rates in longer lived species – which would correspond to genes potentially under positive selection in longer-lived species. We therefore focus our discussion on genes and pathways evolving slower in species with large values for the 3L and ELL traits. We considered these genetic elements to be important to the evolution of the 3L and ELL traits because they have been protected from accumulating deleterious mutations and therefore evolve slower.

Despite considerable previous work on genes under positive selection in species with extreme lifespan traits, our analysis did not recover any positively selected pathways in the pan-mammalian analysis. However, for any new trait to arise, some corresponding genetic changes must occur. It is possible that many of the molecular innovations that lead to increased lifespan are species-specific and thus would not be detected by our method. However, it is also true that our relative evolutionary rates method is much more suited to detecting the kind of long-term convergent increases in purifying selection that we observe. Figure 2 represents a simple schematic history of the evolution of a trait and corresponding rate changes. During the establishment of a new trait, some genes experience a brief period of positive selection which generates the molecular innovation to create the new trait, after which the derived sequence will again be subject to purifying selection. Using extant species as a genomic reference for evolutionary history, we can only infer the average rate over the entire history and the brief period of positive selection results in a small, potentially undetectable signal. On the other hand, we can hypothesize another class of genes which become more important as a trait is established and thus experience continuously increased purifying selection after the trait exists. These are the trait ‘enabling’ as opposed to the trait ‘establishing’ genes. Because the rate change in the case of trait-enabling genes is permanent, it has a greater impact on the inferred average rate observable from extant data. Our analyses strongly support that in the case of extended longevity, such ‘enabling’ genes not only exist but are also convergent across independent trait change events.

We have interpreted the rate correlation from the prospective of species with large phenotype-values, thereby assuming phenotype changes and rate changes coincide. The implicit choice of default rate does not affect the final gene ranking on which our analysis is based, but it does have one important consequence. While decreases in rate are easily interpretable as increased constraint, rate increase can be due to relaxation of constraint or positive selection. We address this asymmetry simply by testing both possibilities. We test the relative accelerated genes for positive selection in the long-lived species and the relatively decelerated genes for positive selection in the short-lived species. We find no evidence of systematic positive selection in either direction.

We interpret our strongest pathway enrichment signals as indicating that cancer control is important for enabling evolution of the 3L phenotype, which is in agreement with previous work that has found similar changes in expression levels of cell cycle and immune function genes in both cancerous and aging cells (Chatsirisupachai et al., 2019). However, there may be alternative explanations. Specifically, differential pathogen pressure is expected to affect some of the same pathways. The connection to immune pathways is clear, and viruses are known to co-opt the cell cycle. TP53 and cell death are likewise important mechanisms by which cells control viral infection. All of these pathways may thus be under increased constraint in species more likely to experience viral infections. However, we believe these explanations are unlikely for several reasons. Firstly, we do not see an enrichment of genes that are specialized for virus control, such as MX1, OAS1, DHX58, and genes activated by type I interferon. Secondly, the interaction between pathogen pressure and evolutionary rates has been studied extensively and it is found that it typically drives adaptive changes in pathogen interacting protein as pathogens are often unique to a species (Kosiol et al., 2008; Sackton et al., 2007; Shultz and Sackton, 2019). This is at odds with the increased constraint of immune genes that we observe. Finally, the probability of cancer increases with age and there is theoretical and empirical evidence that importance of cancer resistance increases with body-size and lifespan. On the other hand, for humans the probability of dying from many infections decreases sharply in early life (though increases again at post-reproductive ages) (Palmer et al., 2018) suggesting that infection resistance should not be preferentially selected for in long-lived species.

Overall, the genes and pathways we identified whose constraint is negatively correlated with the 3L and ELL traits encompass functionalities important for the evolution of extended lifespan, and they therefore represent candidate genes and pathways for further experimental exploration. Importantly, such genes were uncovered using an unbiased, genome-wide pan-mammalian scan. As such, these results point to keys to exceptional longevity that are not specific to one or a handful of species, but that are universal across mammals.

Materials and methods

RERconverge measures relationships between relative evolutionary rates of orthologous genes and phenotype values based on a set of gene trees whose branch lengths represent protein evolutionary rates and a set of phenotype values, either binary or continuous. The RERconverge package is freely-available for use at https://github.com/nclark-lab/RERconverge along with walkthroughs for beginner users (Kowalczyk et al., 2018). The original RERconverge method was successfully used to find genes undergoing convergent evolutionary rate shifts in marine mammals (Chikina et al., 2016) and subterranean mammals (Partha et al., 2017), and recent statistical improvements have made RERconverge more robust and have given it more power to detect such rate shifts (Partha et al., 2018). Here, the RERconverge methods were extended to use on continuous traits.

To perform the RERconverge analysis, we first used PAML (Yang, 2007) to create maximum likelihood gene trees whose branch lengths represent evolutionary rates by means of number of amino acid substitutions. Gene trees were generated for 19,149 amino acid alignments of orthologous genes from the 62 mammal species from the UCSC 100-way alignment (Blanchette et al., 2004; Harris, 2007; Kent et al., 2002) using the hg19 genome, all of which are available at http://genome.ucsc.edu/. Trees were read into R using the readTrees() function from RERconverge. This step also estimated a master tree with branch lengths that represented the average branch lengths across many gene trees. For all further steps, only 61 species were used as listed in Supplementary file 1; the cape golden mole was excluded because longevity data were not available for it. Relative evolutionary rates corrected for genome-wide evolutionary rates for each species as well as corrected for branch length heteroskedasticity using a weighted regression approach (Partha et al., 2018) were calculated using getAllResiduals() with weight = T, scale = T, cutoff = 0.001, and useSpecies = species names for the 61 mammals species.

Relative evolutionary rates (RERs) quantify the deviation in evolutionary rate of a protein along a specific phylogenetic branch from proteome-wide expectations. Negative RERs indicate fewer substitutions than expected due to increased constraint. Positive RERs correspond to more substitutions than expected, which could arise due to relaxation of constraint or positive selection. To correlate RERs with lifespan phenotypes, we use phenotypic change along phylogenetic branches computed from maximum likelihood ancestral state reconstruction (Revell, 2012). This transformation is equivalent to phylogenetically independent contrasts (Felsenstein, 1985) and thus removes phylogenetic dependence from the phenotype values.

We used the char2Paths() function on each phenotype with trees read in using readTrees() to create a phenotype vector that represented the predicted difference in phenotype values between each species and its ancestor. The char2Paths() function first uses maximum likelihood estimation through the fastAnc() function from phytools (Revell, 2012) and then subtracts the values between connected pairs of species to calculate a phylogenetically-independent measure of the change in the phenotype over evolutionary time. Finally, the getAllCor() function was run with each phenotype and the RER matrix with method = ‘p’, min.pos = 0, and winsorizeRER = 3 to perform Pearson correlations with RERs winsorized such that the two most extreme RER values are set to the value of the third most extreme RER. See Figure 1—figure supplement 4 for evolutionary rate-phenotype scatterplots with high correlation values.

The gene lists produced from the RERconverge correlation analyses were used to calculate pathway enrichments for all canonical pathways from mSigDB (Liberzon et al., 2011) and Mouse Genome Informatics (MGI) functional annotations (Blake et al., 2003; Eppig et al., 2015). For each pathway, a Wilcoxon Rank-Sum statistic was calculated to compare the sign of Rho times the negative log p-value for correlations of genes in the pathway to the same measure for all genes included in a pathway annotation using getStat() from RERconverge to calculate the sign of Rho times the negative log p-value for correlations for each gene and fastwilcoxGMTall() to quickly calculate an approximation of the Wilcoxon Rank-Sum statistic.

Permulation analysis (Phylogenetically-Restricted Permutations)

Request a detailed protocol

In addition to calculating theoretical enrichment statistics, we developed a novel phylogenetically-restricted permutation strategy dubbed ‘permulations’ to calculate empirical p-values for each pathway enrichment statistic. Permulations are a combination of permutations and simulations in reference to the strategy for generating null phenotype values to use to generate null statistics. To perform permulations, first phenotype values (3L and ELL) were simulated for each species using phylogenetic simulations, and then original phenotype values were reassigned to the species based on the rank of the simulated values (see Figure 1—figure supplement 5). Simulations were performed using the geiger library in R using a Brownian motion approach and the average evolutionary rate tree created by RERconverge (Harmon et al., 2008). After creating new phenotype values, RERconverge analyses were performed to calculate correlations between evolutionary rates of genes and phenotype values, and enrichment statistics were calculated using these gene results. The process of creating new phenotype values and calculating enrichment statistics was repeated 1000 times for each phenotype and the empirical p-value for pathway enrichment was calculated as the proportion of times the permulation enrichment statistic was greater than the enrichment statistic calculated using real phenotype values.

Permulations are an attractive option to calculate empirical p-values because they use fabricated phenotypes that are independent of RERs but respect the underlying phylogenetic relationships between species (due to use of phylogenetic simulations) while maintaining the original range of the data (due to rank-based assignment of true phenotype values). These types of analyses are a common tool in genomics because they allow for control of the dependence across tests caused by intercorrelations among genes. Permulations are highly analogous to permutations commonly performed in differential expression analysis, and we use them for the same purpose, namely to calculate empirical pathway-level statistics (Subramanian et al., 2005). This step is critical for evaluating pathway enrichment results because many pathways have elevated RER correlations even when conditioned on the phenotype of interest. A pathway with high RER correlation among its genes is more likely to show up as enriched when using a test that assumes gene independence, such as the Wilcoxon Rank-Sum test. Permulations allow us to generate an empirical null distribution for pathway enrichment statistics to correct for interdependence among gene ranks. Indeed, we find multiple pathways that show significant enrichment using a Wilcoxon Rank-Sum test but insignificant empirical p-values using permulations. We also show that permulation p-values, which take phylogeny into account, are more conservative than permutation p-values, which ignore phylogenetic dependence, and equally as conservative as the phylogenetic simulation p-values (see Figure 1—figure supplement 6).

Positive selection tests

Request a detailed protocol

In genes significantly correlated with the 3L and the ELL phenotypes, we investigated evidence for relaxation of constraint and positive selection on trait-defining foreground branches using phylogenetic models of codon evolution. We did so using a representative subset of the full mammalian phylogeny (see Figure 1—figure supplement 7). Trait-defining foreground branches were specified independently for the two phenotypes along both axes of trait values – positive and negative. In total we have four sets of trait-defining foreground branches, namely positive 3L, negative 3L, positive ELL, and negative ELL. Figure 1—figure supplement 8 shows the phylogeny of species used for this analysis and the four sets of trait-defining foreground branches. Species were selected for the representative subset based on three criteria: 1) all species were placental mammals (the monotreme and marsupials were excluded), 2) species with highest and lowest 3L and ELL phenotype values were included and used as foreground species in their respective tests, and 3) closely-related non-foreground species were included. Only placental mammals were used because the inference of positive selection can be confounded in the non-placental clade due to the long divergence time between the clades and hence saturation of synonymous sites for many genes. Foreground species were selected based on 3L and ELL values. A subset of non-foreground species was selected to correct for over- and under-sampling in particular clades while maintaining valid ‘outgroups’ for foreground species. This subset also eliminated species that were essentially duplicates in terms of tests for positive selection. For example, the primate clade is sampled much more than other clades, that is more primate genomes are available. Additionally, because most of these species are non-foreground species, their presence has little impact on positive selection tests. To correct for the oversampling, we used five out of twelve primates in our subtree, three of which were foreground species (human, marmoset, and squirrel monkey) and two of which were outgroups (rhesus as an outgroup to human and bushbaby as an outgroup to the whole clade). We used similar logic to select outgroup species throughout the rest of the phylogeny.

We inferred the significance of relaxation of constraint on each foreground branch set using likelihood ratio tests (LRT) between Branch-site Neutral (BS Neutral) and its nested null model M1 (sites neutral model) in PAML (Yang, 2007). Similarly, we performed LRTs between branch-site selection model (BS Alt Mod) and its null BS Neutral were used to infer positive selection on the foreground branches. Probabilities for each of these two LRTs were estimated using the chi-square distribution with 1 degree of freedom. We additionally inferred significance of mammal-wide relaxation of constraint and positive selection using the LRTs between M8A (Neutral model) vs M0 (null model) and M8 (positive selection model) vs M8A respectively. Prior to performing the mammal-wide tests in genes corresponding to each of the four foreground branch sets, we removed the corresponding foreground branches, allowing us to obtain unbiased estimates for significance of relaxation of constraint and positive selection from only the background mammalian branches.

Species robustness through subtree analysis

Request a detailed protocol

We were interested in assessing the sensitivity of our results to the choice of mammalian species used in our analyses. The mammalian genomes available represent a subset of not only all extant mammals, but all mammal species that have ever existed. Therefore, the genomes used in this analysis represent an incomplete, and perhaps even a biased representation of mammal species. Since we would like to extend our conclusions to pertain to mammals in general, we sought to quantify the effects of our incomplete data on our gene and pathway results.

To do this, we created subsets of our data that contained fewer species. These subsets had 10, 20, 30, 40, 50, 60, 70, and 80 percent of species randomly removed (6, 12, 18, 24, 30, 36, 42, and 48 species removed out of 61 total, respectively) with ten random subsets created for each species removal level. We then ran the standard RERconverge analysis on these subsets to acquire correlation statistics representing the relationship between evolutionary rate of each gene across species and longevity phenotypes in those species. We also ran enrichment analyses on the subset gene results to acquire pathway enrichment statistics. After performing these analyses, we calculated correlations between results from our full dataset and results from the subsets, where results were quantified as the negative log p-value times the sign of the statistic.

We further tested the sensitivity of our results to specific species presence/absence by performing targeted species removal of species groups with which we expected to have other phenotypes confounded with body size, namely marine mammals, and species that may have non-convergent genetic mechanisms for lifespan extension, namely bats and the naked mole rat. We created two new data subsets that only contained non-marine species (dolphin, manatee, killer whale, walrus, and Weddell seal were removed) and only non-bat and non-naked-mole-rat species (megabat, black flying-fox, microbat, David’s myotis bat, big brown bat, and naked mole-rat were removed). We then performed the standard RERconverge analysis to find correlations between the evolutionary rates of genes across species and longevity phenotypes in those species. We also performed pathway enrichment analyses on the gene results.

We found good correlation of both enrichment and correlation statistics between subset results and full dataset results based on negative log p-values times the sign of the statistic (Figure 1—figure supplement 3). For data subsets in which a proportion of species were randomly eliminated, there was upwards of a median 60% correlation between results even with half all species removed from our analyses, which indicates that randomly removing species did not significantly change our results. These findings are quantified numerically as an among-subtree variance in Supplementary file 4 and Supplementary file 5. We also compared full dataset results to results from targeted data subsets without marine mammals and without bats and the naked mole-rat by quantifying correlations between their negative log p-values times the sign of the statistic for gene correlations and pathway enrichment statistics (Figure 1—figure supplement 8). From these comparisons, we found a strong relationship between full dataset results and targeted subset results, which indicates that presence or absence of potentially problematic species such as marine mammals, bats, and the naked mole-rat is not strongly impacting our results. Quantitative representations of these findings are available in Supplementary file 4 and Supplementary file 5 Together, these findings indicate that our results are not species, clade, or species subgroup-specific, but instead represent trends across all mammals.

Alternate tree topology analysis

Request a detailed protocol

In addition to determining whether individual or groups of species were disproportionately driving our conclusions, we were also interested in verifying that the fixed tree topology used for analyses was not affecting results. To do so, we reran analyses using alternate plausible topologies representing different ancestral relationships among species (trees shown in Figure 1—figure supplement 9). Differences in the alternate tree topologies represent points of potential incomplete lineage sorting that lead to uncertainty, they and are thus reasonable alternatives compared to the original tree used for analyses.

RERconverge was run using gene trees generated using the alternate topologies and results were compared to results using the original tree. As shown in Figure 1—figure supplement 10, there is a strong correlation between all sets of results, which indicates that uncertainties in tree topology do not strongly affect results. This is true for very similar topologies (Robinson-Foulds distance 6) and fairly different topologies (Robinson-Foulds distance 22).

Phylogenetic trees

Request a detailed protocol

Full Tree (61 mammals):

((((((((((((Human:0.005957477577,Chimp:0.006721826689):0.001382639829,Gorilla:0.007765177171):0.005572327638,Orangutan:0.0164503644):0.002187630666,Gibbon:0.01770384793):0.007043113559,(Green_monkey:0.007693724903,((Crab-eating_macaque:0.001292320552,Rhesus:0.00713015786):0.002951690224,Baboon:0.005199240711):0.002049749893):0.01566263562):0.0135408115,(Marmoset:0.02474184521,Squirrel_monkey:0.02096868307):0.02784675729):0.04299750653,Bushbaby:0.108738222):0.01379370868,((((((Guinea_pig:0.09048639907,(Chinchilla:0.05332953299,Brush-tailed_rat:0.08476954109):0.01287861561):0.02118937782,Naked_mole-rat:0.08588673524):0.07432515556,Squirrel:0.08896424642):0.006291577528,((((Chinese_hamster:0.04084640027,Golden_hamster:0.04456203524):0.02314125062,Prairie_vole:0.06932402649):0.01947113467,(Mouse:0.05273642272,Rat:0.05576007402):0.04435347588):0.08380065137,Lesser_Egyptian_jerboa:0.1438649666):0.04270536633):0.01663675397,(Pika:0.1256544445,Rabbit:0.07131655591):0.06535533418):0.009050428462,Chinese_tree_shrew:0.1191189141):0.003894252213):0.01425600689,(((((Panda:0.03854019703,((Weddell_seal:0.02002160645,Pacific_walrus:0.02064385875):0.01734764946,Ferret:0.04613997497):0.002879093616):0.009005888384,Dog:0.05339127565):0.01185166857,Cat:0.05020331605):0.03285617057,((((((Cow:0.02168740723,((Domestic_goat:0.01157093136,Sheep:0.01246322594):0.0049716126,Tibetan_antelope:0.01522587482):0.01465511149):0.0662523666,(Killer_whale:0.006371664911,Dolphin:0.01086552617):0.06014682602):0.01216198069,Pig:0.0796745271):0.006785823323,(Bactrian_camel:0.01240650215,Alpaca:0.01096629635):0.06374554586):0.02551888691,(White_rhinoceros:0.04977357056,Horse:0.061454379):0.02510111297):0.00331214686,((Big_brown_bat:0.03248546656,(Davids_Myotis_bat:0.02344332842,Microbat:0.01567729315):0.02193849809):0.09455328094,(Black_flying-fox:0.005833353548,Megabat:0.01611220178):0.07567400302):0.02385546003):0.002057771224):0.004845253848,(Star-nosed_mole:0.1239823369,(Hedgehog:0.1696142244,Shrew:0.1934205791):0.02079474546):0.0235875333):0.01477733374):0.01316436193,(((((Cape_golden_mole:0.1017903453,Tenrec:0.1749615473):0.01592632003,Cape_elephant_shrew:0.1516860647):0.006610995228,Aardvark:0.08326528894):0.008243787904,(Elephant:0.06812658238,Manatee:0.06198982615):0.0224994529):0.03384011363,Armadillo:0.1342602666):0.005989703247):0.2206952867,((Wallaby:0.1270943532,Tasmanian_devil:0.09944141622):0.02717055443,Opossum:0.1181200712):0.1802966572):0,Platypus:0.4322118716);

Branch Site Tree (34 placental mammals):

((((((Human:0.02214318927,Rhesus:0.0277942336):0.0135408115,(Marmoset:0.02474184521,Squirrel_monkey:0.02096868307):0.02784675729):0.04299750653,Bushbaby:0.108738222):0.01379370868,(((Chinchilla:0.08739752642,Naked_mole-rat:0.08588673524):0.07432515556,Squirrel:0.08896424642):0.006291577528,((((Chinese_hamster:0.04084640027,Golden_hamster:0.04456203524):0.02314125062,Prairie_vole:0.06932402649):0.01947113467,(Mouse:0.05273642272,Rat:0.05576007402):0.04435347588):0.08380065137,Lesser_Egyptian_jerboa:0.1438649666):0.04270536633):0.02958143464):0.01425600689,((((Weddell_seal:0.02002160645,Pacific_walrus:0.02064385875):0.04108430003,Cat:0.05020331605):0.03285617057,(((((Cow:0.02168740723,(Sheep:0.01743483854,Tibetan_antelope:0.01522587482):0.01465511149):0.0662523666,Killer_whale:0.06651849094):0.01216198069,Pig:0.0796745271):0.03230471024,Horse:0.08655549197):0.00331214686,((Big_brown_bat:0.03248546656,(Davids_Myotis_bat:0.02344332842,Microbat:0.01567729315):0.02193849809):0.09455328094,Black_flying-fox:0.08150735657):0.02385546003):0.002057771224):0.004845253848,(Star-nosed_mole:0.1239823369,(Hedgehog:0.1696142244,Shrew:0.1934205791):0.02079474546):0.0235875333):0.01477733374):0.01316436193,((Cape_elephant_shrew:0.1665408478,(Elephant:0.06812658238,Manatee:0.06198982615):0.0224994529):0.03384011363,Armadillo:0.1342602666):0.005989703247);

Data and materials availability

Request a detailed protocol

All data and code are publicly available, provided in the main text or supplementary materials, or available through the RERconverge package on GitHub at https://github.com/nclark-lab/RERconverge. (Kowalczyk et al., 2019; copy archived at https://github.com/elifesciences-publications/RERconverge).

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
    Improved Pairwise Alignment of Genomic DNA
    1. RS Harris
    (2007)
    Pennsylvania State University.
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
    An Unsolved Problem in Biology
    1. PB Medawar
    (1952)
    HK Lewis, editors. London: Scientific Research.
  31. 31
  32. 32
  33. 33
  34. 34
    Walker’s Mammals of the World (6th ed)
    1. RM Nowak
    (1999)
    Baltimore: The Johns Hopkins University Press.
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
    Chromatin dynamics in DNA double-strand break repair
    1. L Shi
    2. P Oberdoerffer
    (2012)
    Biochimica Et Biophysica Acta (BBA) - Gene Regulatory Mechanisms 1819:811–819.
    https://doi.org/10.1016/j.bbagrm.2012.01.002
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59

Decision letter

  1. Matt Kaeberlein
    Reviewing Editor; University of Washington, United States
  2. Diethard Tautz
    Senior Editor; Max-Planck Institute for Evolutionary Biology, Germany

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

This study provides important insights into the evolution of lifespan across a range of mammalian species and identifies several genes under increased evolutionary constraint in long-lived species.

Decision letter after peer review:

Thank you for submitting your article "Pan-mammalian analysis of molecular constraints underlying extended lifespan" for consideration by eLife. Your article has been reviewed by two peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Diethard Tautz as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

Kowalczyk and colleagues report on a comparative genomics approach to understand potential mechanisms of lifespan determination in mammals. They consider protein evolutionary rates across a wide range of mammalian species and identify signatures suggesting that genes related to cell cycle, DNA repair, cell death, the IGF1 pathway, and immunity are under evolutionary constraint in long-lived mammalian species. The unbiased approach is a particular strength of this study, as much of the prior work in this area has involved targeted studies aimed at querying specific mechanisms/pathways in long-lived species. Although the conclusions are a bit short on new insight into potential conserved mechanisms of longevity determination, they generally make biological sense and add further support for the necessity of robust anti-cancer mechanisms to evolve in parallel with long lifespans in large mammals.

Essential revisions:

1) Parts of the selection analysis should be expanded to consider shorter lifespans. A first assumption seems to be that longer lifespans tend to be achieved through positive selection; while shorter lifespans tend to be achieved through relaxed purifying selection. This fits with Mutation Accumulation and Antagonistic Pleiotropy theories in the sense that shorter lifespans would evolve due to the relative lack of efficient purifying selection at older ages; while one would need to fix advantageous mutations to lengthen lifespan. This seems correct in general, even if, to the best of our knowledge, no one has shown this to always be the case. A second assumption seems to be that once long lifespans have been achieved, purifying selection will keep them this way and, hence, they will come with would come with constrained rates of protein evolution in similar pathways. That is what they set out to detect, since it would be indicative of convergent evolution at the molecular level.

However, there is the additional possibility that once a species has achieved a shorter lifespan (for instance, due to a trade-off to increase early-age reproduction) purifying selection keeps it this way. It may be the case, thus, that increased constraint would also be detected in similar pathways in species with lower 3L and ELL values. We think this study is particularly suited to address this point with some rigor. So please, look at "more constrained in short lived" pathways and tell us what you see.

2) The multiple testing correction needs a better explanation. Reading the Materials and methods, it would seem that "permulations" only apply to pathways, "…to calculate empirical values for each pathway enrichment statistic", only pathway enrichment is mentioned there). However, it seems (from the first paragraph in the Results section and from Figures 1C and 1D) that the first thing done with "permulations" is to see if they have an excess of significant p-values in the set of individual gene correlations. Is that so? We think clarifying this is critical since, of course, any distribution of so many correlations between RERs and phenotypes will have extremes and significant p-values. Incidentally, we would suggest including conventional QQ-plots, at least as supplementary information.

3) Related to the previous point, we are fine with using RERs, 3L and ELL. However, we are not fully convinced that "permulations" are so conservative. First, using the "average evolutionary rate tree" is likely to generate mismatches with closely related species where incomplete lineage sorting or gene flow have been relevant phenomena. In these cases, the distribution of observed correlation p-values may be different from the simulated (permulated) distributions. Second, we know that rates of evolution (of both, molecules and phenotypes) do not follow an exact clock.

https://doi.org/10.7554/eLife.51089.sa1

Author response

Essential revisions:

1) Parts of the selection analysis should be expanded to consider shorter lifespans. A first assumption seems to be that longer lifespans tend to be achieved through positive selection; while shorter lifespans tend to be achieved through relaxed purifying selection. This fits with Mutation Accumulation and Antagonistic Pleiotropy theories in the sense that shorter lifespans would evolve due to the relative lack of efficient purifying selection at older ages; while one would need to fix advantageous mutations to lengthen lifespan. This seems correct in general, even if, to the best of our knowledge, no one has shown this to always be the case. A second assumption seems to be that once long lifespans have been achieved, purifying selection will keep them this way and, hence, they will come with would come with constrained rates of protein evolution in similar pathways. That is what they set out to detect, since it would be indicative of convergent evolution at the molecular level.

However, there is the additional possibility that once a species has achieved a shorter lifespan (for instance, due to a trade-off to increase early-age reproduction) purifying selection keeps it this way. It may be the case, thus, that increased constraint would also be detected in similar pathways in species with lower 3L and ELL values. We think this study is particularly suited to address this point with some rigor. So please, look at "more constrained in short lived" pathways and tell us what you see.

We appreciate the reviewers’ thoroughness and extensive knowledge of aging-related theories, and we have spent considerable time discussing these ideas. One key finding in this study was that levels of purifying selection, or conservation, across species with a range of 3L and ELL phenotypes represent relative gene importance. Our broad hypothesis based on findings in this work is that latent positive selection at some point in evolutionary history drove establishment of the 3L and ELL phenotypes, while selective pressure shifts associated with trait establishment caused sustained increased purifying selection of genes supporting the newly-established trait. Since RERconverge detects average rate changes over evolutionary time, we are unable to detect periods of brief positive selection, while we are much more readily able to detect sustained purifying selection. A visual of this hypothesis is shown in Figure 1—figure supplement 1. Results from both RERconverge and branch-site models to test for positive selection were unable to detect instances of positive selection support this model. It is also possible that instances of positive selection were non-convergent, in which case they would neither be detectable to RERconverge nor relevant to generalizable conclusions about the origin of extended lifespan in mammals. We performed branch-site model tests for positive selection using species with both large and small values of both 3L and ELL phenotypes as foreground species to test genes accelerated in species with both large and small values of both 3L and ELL phenotypes, thus accounting for both correlation directions. None of the branch-site model tests detected meaningful positive selection in any foreground species, and broadly, RERconverge detected more significant negative correlations (genes evolving faster in small, short-lived species and slower in large, long-lived species) than positive correlations (genes evolving slower in small, short-lived species and faster in large, long-lived species), with no significantly enriched pathways evolving faster in species with large values of the 3L and ELL phenotypes. Therefore, although it is theoretically possible that genes could be under increased purifying selection to increase phenotypes related to evolution of small size and short lifespan, we do not observe such genes in our analysis.

Further, since long lifespan and large body size are inferred to be derived traits (the ancestral mammal is believed to have been small and short lived), evolutionary states of species with large values of 3L and ELL are likely representative of changes related to evolution of those traits, while evolutionary states of species with small values of 3L and ELL are likely representative of maintenance of ancestral states. This may explain why we do not observe these positive correlations in our analyses.

Mutation Accumulation and Antagonistic Pleiotropy theories broadly apply to lifespan on the scale of an individual rather than changes in lifespan over evolutionary time. These theories state, respectively, that accumulation of mutations throughout an individual’s lifetime eventually leads to breakdown of essential functions that drive senescence and eventual age-related death, and that functions that benefit early lifetime survival and fitness are actually detrimental in old age. We do not believe that potential reduced mutation rates in old age in large and long-lived mammals are due to effective increased purifying selection in these species, but instead that increased purifying selection on targeted genes, such as those for DNA repair and cell cycle control, facilitates improved functionalities of these genes throughout an individual’s lifespan and thus increased conservation indirectly contributes to existing hypotheses regarding aging.

We have updated the Results section as follows to further clarify that we addressed both correlation directionalities and considered both possible directions of positive selection:

“Our analysis investigates both positive and negative correlations between evolutionary rates of genes and changes in lifespan phenotypes. […] Both the gene-level and pathway-level results were highly robust to species removal (see Materials and methods, Figure 1—figure supplement 6, Supplementary file 4, and Supplementary file 5).”

We have also added Figure 1—figure supplement 1 to clarify why we detect more genetic signal for pathways under increased purifying selection than positive selection.

2) The multiple testing correction needs a better explanation. Reading the Materials and methods, it would seem that "permulations" only apply to pathways, "…to calculate empirical values for each pathway enrichment statistic", only pathway enrichment is mentioned there). However, it seems (from the first paragraph in the Results section and from Figures 1C and 1D) that the first thing done with "permulations" is to see if they have an excess of significant p-values in the set of individual gene correlations. Is that so? We think clarifying this is critical since, of course, any distribution of so many correlations between RERs and phenotypes will have extremes and significant p-values. Incidentally, we would suggest including conventional QQ-plots, at least as supplementary information.

We agree with reviewers that multiple hypothesis testing correction is a key part of our analyses since performing nearly 20,000 tests for genes will lead to false positives. Reviewers are also correct that permulation p-values were only calculated for pathways, although null permulation p-values were generated both for genes and pathways, as they are necessary in order to calculate permulation p-values (permulation p-values are the proportion of null permulation statistics more extreme than an observed statistic for a particular test, be it gene-level or pathway level). At the gene level, we perform classic multiple hypothesis testing correcting using the Benjamini-Hochberg correction. We compare this multiple hypothesis testing correction, which uses a standard uniform distribution to represent the assumed proportion of false positive p-values, to multiple hypothesis testing correction using the true empirical null hypothesis for this dataset as calculated using permulation null p-values, or p-values generated by shuffling phenotype values with respect to the underlying phylogeny. Based on this empirical null distribution generated using null permulation p-values, we find that our theoretical adjusted p-values are actually more conservative than we would expect because the distribution of empirical null p-values is slightly non-uniform and in fact slopes downward at the left end of the p-value distribution (Figures 1D and 1E).

At the pathway level, we do calculate empirical p-values on a per-pathway basis. These p-values represent the proportion of null permulation enrichment statistics that are more extreme than the observed pathway enrichment statistics, akin to bootstrap p-values. Again, these permulation p-values are multiple hypothesis testing corrected using a standard Benjamini-Hochberg correction.

Permutation analyses are a common tool in genomics because they allow one to control for the dependence across tests. They are ubiquitously used in differential expression analysis to calculate empirical FDRs for gene-level test and empirical p-values for pathway-level test. In the simplest case we are concerned with comparing two groups with different clinical phenotypes. The null hypothesis is that the phenotype differences are not reflected in the data and thus with respect to the data the two groups are interchangeable. Consequently, permutations of group labels can be used to generate an empirical null distribution while controlling for correlations across genes (which remain unperturbed). This empirical null can then be used to calculate empirical FDRs for gene level tests. Correlation among genes also must be considered for evaluating pathway level analysis. Many tests for pathway enrichment assume that genes are independent when conditioned on the phenotype, which is clearly false. Genes that are co-expressed tend to have similar differential expression statistics. Permutations can be used to generate per-pathway null test statistic distribution by controlling for correlations among pathway genes. This distribution can then be used to calculate an empirical p-value for pathway enrichment, which will typically be larger than those derived from a parametric test that assumes independence.

Our permutations pipeline follows closely what would typically be done in a gene expression context. However, in our case we cannot simply permute the species values because even under the null hypothesis these are not independent. To extend the differential expression analogy, if some patients were profiled more than once any valid permutation would need to place those samples in the same group. In our case, the species are related through the phylogeny, so we perform phylogenetically constrained permutations that ensure that the species that are phylogenetically closer are also more likely to receive similar phenotype values. As expected, this results in more conservative p-values than the unconstrained permutations (Figure 1—figure supplement 9).

The main text has been altered as follows to further clarify the role of permulations in hypothesis testing:

“Permulations are an attractive option to calculate empirical p-values because they use fabricated phenotypes that are independent of RERs but respect the underlying phylogenetic relationships between species (due to use of phylogenetic simulations) while maintaining the original range of the data (due to rank-based assignment of true phenotype values). […] We also show that permulation p-values, which take phylogeny into account, are more conservative than permutation p-values, which ignore phylogenetic dependence, and equally as conservative as the phylogenetic simulation p-values (see Figure 1—figure supplement 9).”

We have also added Q-Q plots as a supplementary figure as requested (Figure 1—figure supplement 10).

3) Related to the previous point, we are fine with using RERs, 3L and ELL. However, we are not fully convinced that "permulations" are so conservative. First, using the "average evolutionary rate tree" is likely to generate mismatches with closely related species where incomplete lineage sorting or gene flow have been relevant phenomena. In these cases, the distribution of observed correlation p-values may be different from the simulated (permulated) distributions. Second, we know that rates of evolution (of both, molecules and phenotypes) do not follow an exact clock.

We agree with reviewers that rates of evolution do not follow an exact clock, and we fully take advantage of that in order to perform our analyses to detect evolutionary rate shifts associated with phenotypic changes. Also note that we take into account broad differences in evolutionary rates on the species level by normalizing for differences in species evolutionary rates caused by factors such as generation times by quantifying evolution as a rate relative to species average evolutionary rate (thus relative evolutionary rates as RERs).

We additionally fully subscribe to the notion that a fixed phylogeny is a model, and although all models are wrong, some are useful and some are better than others. We strived to create the best phylogeny possible by incorporating information from two highly-regarded studies that built phylogenies as well as our own literature review of highly-contested relationships among species. Although most of the species we considered are so distantly related that we don’t expect many completely incorrect relationships in extant species, there undoubtedly was incomplete lineage sorting on internal branches of the phylogeny. In our analyses, this circumstance would be reflected by extremely short branch lengths representing little evolutionary change that would be filtered out prior to statistical analyses, thus making the analysis robust to phylogenetic inaccuracies. We have demonstrated this robustness by performing our analyses using two different phylogenies, one with a Robinson-Fould distance of 6 and one with a distance of 22 from the phylogeny used in the original study (Figure 1—figure supplement 11). These distances represent the number of “swaps” among branches that occurred in the phylogenetic topology. As highlighted in the supplementary figure, locations of swaps are at points of uncertainty, such as the relationship between the Chinese tree shrew and primates, that likely exist because of phenomena such as incomplete lineage sorting that cause confusion among species at the molecular level. There is a strong correlation between both gene-level and pathway-level statistics between the Meredith+ topology and the alternative topologies (Figure 1—figure supplement 12), which indicates that the reported results are robust to changes in tree topology representing topological uncertainties and instances of potential incomplete lineage sorting.

We have included Figure 1—figure supplements 11 and 12 to demonstrate the tree topologies tested and the correlations between results from different tree topologies.

We have also added the following text to the manuscript to describe the alternate tree topology analyses:

“Alternate Tree Topology Analysis

In addition to determining whether individual or groups of species were disproportionately driving our conclusions, we were also interested in verifying that the fixed tree topology used for analyses was not affecting results. […] This is true for very similar topologies (Robinson-Foulds distance 6) and fairly different topologies (Robinson-Foulds distance 22).”

In regards to the relationship between permulation p-values and theoretical p-values, we expect and observe differences between the two values at the level of individual pathways. Some pathways that exhibit strong enrichment according to theoretical p-values have non-significant permulation p-values, likely because the permulation p-values account for confounders related to the structure of the data that we as statisticians do not anticipate. The beauty of permulation p-values is that they represent the true null expectation for the given dataset – under the given phylogeny and set of statistical tests, we can know what the null expectation is for the statistics because we have generated thousands of “fake” (permulation) phenotype values that should not have significant statistics (because they are literally fake phenotypes). For these “fake” phenotypes, we do not expect to observe significant correlations with genes and subsequent significant pathway enrichment except in false positive cases, which may or may not occur at a rate of about 0.05 (α). If, for a particular pathway, false positive cases happen at random more often than we would expect (greater than α=0.05 of the time), that is quantified by the permulation p-values – in that case, more null permulation p-values would be significant, more null permulation statistics would be more extreme than the observed statistic, and thus the calculated empirical permulation p-value for the pathway would be higher (less significant) because we know that the probability of getting a false positive significant p-value is higher than the typical null expectation.

https://doi.org/10.7554/eLife.51089.sa2

Article and author information

Author details

  1. Amanda Kowalczyk

    1. Joint Carnegie Mellon University-University of Pittsburgh PhD Program in Computational Biology, Pittsburgh, United States
    2. Department of Computational and Systems Biology, University of Pittsburgh, Pittsburgh, United States
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Investigation, Visualization, Methodology
    For correspondence
    kowaae22@pitt.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9061-1336
  2. Raghavendran Partha

    1. Joint Carnegie Mellon University-University of Pittsburgh PhD Program in Computational Biology, Pittsburgh, United States
    2. Department of Computational and Systems Biology, University of Pittsburgh, Pittsburgh, United States
    Contribution
    Conceptualization, Resources, Data curation, Software, Formal analysis, Supervision, Funding acquisition, Visualization, Methodology, Project administration
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7900-4375
  3. Nathan L Clark

    1. Department of Computational and Systems Biology, University of Pittsburgh, Pittsburgh, United States
    2. Pittsburgh Center for Evolutionary Biology and Medicine, University of Pittsburgh, Pittsburgh, United States
    3. Department of Human Genetics, University of Utah, Salt Lake City, United States
    Contribution
    Conceptualization, Resources, Supervision, Funding acquisition, Project administration
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0006-8374
  4. Maria Chikina

    Department of Computational and Systems Biology, University of Pittsburgh, Pittsburgh, United States
    Contribution
    Conceptualization, Resources, Data curation, Software, Supervision, Funding acquisition, Visualization, Methodology, Project administration
    For correspondence
    mchikina@pitt.edu
    Competing interests
    No competing interests declared

Funding

National Institutes of Health (R01HG009299)

  • Nathan L Clark
  • Maria Chikina

National Institutes of Health (U54 HG008540)

  • Nathan L Clark
  • Maria Chikina

Howard Hughes Medical Institute (T32 EB009403)

  • Amanda Kowalczyk

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We would like to thank Dr. Andreas Pfenning and Dr. Dennis Kostka for helpful discussion and feedback, as well as all members of the Clark and Chikina labs. 

Senior Editor

  1. Diethard Tautz, Max-Planck Institute for Evolutionary Biology, Germany

Reviewing Editor

  1. Matt Kaeberlein, University of Washington, United States

Publication history

  1. Received: August 14, 2019
  2. Accepted: January 14, 2020
  3. Version of Record published: February 11, 2020 (version 1)

Copyright

© 2020, Kowalczyk et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 983
    Page views
  • 123
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Computational and Systems Biology
    2. Immunology and Inflammation
    Bisrat J Debebe et al.
    Research Article
    1. Computational and Systems Biology
    2. Evolutionary Biology
    Rahil Taujale et al.
    Research Article