1. Evolutionary Biology
  2. Genetics and Genomics

Universal and taxon-specific trends in protein sequences as a function of age

  1. Jennifer E James  Is a corresponding author
  2. Sara M Willis
  3. Paul G Nelson
  4. Catherine Weibel
  5. Luke J Kosinski
  6. Joanna Masel  Is a corresponding author
  1. Department of Ecology and Evolutionary Biology, University of Arizona, United States
  2. Department of Physics, University of Arizona, United States
  3. Department of Mathematics, University of Arizona, United States
  4. Department of Molecular and Cellular Biology, University of Arizona, United States
https://doi.org/10.7554/eLife.57347
Share this article
Cite this article
  1. Jennifer E James
  2. Sara M Willis
  3. Paul G Nelson
  4. Catherine Weibel
  5. Luke J Kosinski
  6. Joanna Masel
(2021)
Universal and taxon-specific trends in protein sequences as a function of age
eLife 10:e57347.
https://doi.org/10.7554/eLife.57347
  2. Download BibTeX
  3. Download .RIS

Abstract

Extant protein-coding sequences span a huge range of ages, from those that emerged only recently to those present in the last universal common ancestor. Because evolution has had less time to act on young sequences, there might be ‘phylostratigraphy’ trends in any properties that evolve slowly with age. A long-term reduction in hydrophobicity and hydrophobic clustering was found in previous, taxonomically restricted studies. Here we perform integrated phylostratigraphy across 435 fully sequenced species, using sensitive HMM methods to detect protein domain homology. We find that the reduction in hydrophobic clustering is universal across lineages. However, only young animal domains have a tendency to have higher structural disorder. Among ancient domains, trends in amino acid composition reflect the order of recruitment into the genetic code, suggesting that the composition of the contemporary descendants of ancient sequences reflects amino acid availability during the earliest stages of life, when these sequences first emerged.

Introduction

New protein-coding genes emerge over time, either through a combination of rearrangement, duplication, and divergence from ancestrally coding sequences, or de novo from previously non-coding DNA (for a review see Van Oss and Carvunis, 2019, Long et al., 2003). Although genes clearly originated de novo during the emergence of life, de novo gene birth after this unique event was once believed to be extremely rare (Jacob, 1977; Keese and Gibbs, 1992; Zuckerkandl, 1975). Today, de novo gene birth is increasingly acknowledged to be real and important, with well-documented examples confirmed in diverse taxa (Baalsrud et al., 2018; Cai et al., 2008; Khalturin et al., 2009; Milde et al., 2009; Zile et al., 2020), despite the technical challenges associated with correctly identifying de novo genes (McLysaght and Hurst, 2016). Compared to genes born early in the history of life, recently emerged de novo genes have had little time to evolve and adapt. If their properties converge only slowly to those of ancient genes, this creates trends in gene properties with age (Neme and Tautz, 2013).

It has been claimed that younger genes encode shorter and faster evolving proteins (Lipman et al., 2002; Toll-Riera et al., 2012) that are more disordered, containing fewer hydrophobic amino acids (Wilson et al., 2017), which are more clustered together along the primary sequence (Foy et al., 2019). Foy et al., 2019 suggest that the two latter trends, in hydrophobicity and in its clustering, together indicate that evolution is slow to find the sophisticated folding strategies employed by older proteins, which minimize their propensity to form aggregates within cells despite relatively high overall hydrophobicity.

The universality of some of these trends has, however, been questioned. In particular, some authors have found that in some taxa, such as Saccharomyces cerevisiae (Carvunis et al., 2012; Vakirlis et al., 2018), young genes are in fact less disordered than older genes. Random sequences with high %GC tend to encode polypeptides with higher intrinsic structural disorder (ISD) (Ángyán et al., 2012). Differences in genomic GC content between species might therefore shape disorder at the time of gene birth differently (Basile et al., 2016; Van Oss and Carvunis, 2019). Because GC content evolves over relatively short timescales, GC content could explain differences in de novo genes among species, but less so among broad taxonomic groups.

Alternatively, the low disorder of young ‘protogenes’ may be an artifact of including sequences that do not meet the evolutionary definition of functionality, which requires that the loss of the gene be deleterious, that is, have a selection coefficient of less than zero (Graur et al., 2013). In a species whose %GC gives rise to low structural disorder in non-functional protogenes, the pooling of high-disorder functional and low-disorder non-functional genes in varying proportions can create a spurious trend with age. Indeed, the direction of the disorder trend in S. cerevisiae was reversed when dubious gene candidates were removed from the data of Carvunis et al., 2012 (Wilson et al., 2017).

Nevertheless, with phylostratigraphy so far studying a limited set of relatively well-annotated focal species, it remains possible that trends in protein properties as a function of age are not universal, and instead depend on taxonomic group. The forces that shape gene birth could differ across lineages because the first imperative to ‘do no harm’ operates differently as a function of differences in proteostasis machinery. Alternatively, different amino acids might be abundant, and hence cheap to use, in different taxa. We therefore set out to test whether trends in protein properties are the same vs. different across sequences of different ages and across different species, and which forces might underlie any differences.

Identifying phylostratigraphy trends as a function of gene age requires age estimation. Gene ages are based on the date of the most basal node in the phylogeny of lineages containing homologs (Domazet-Loso et al., 2007). But when sequences are highly divergent, programs used to detect homologs such as BLASTp (Altschul et al., 1990) are prone to false negatives, and thus underestimate gene age (Elhaik et al., 2006; Jain et al., 2019; McLysaght and Hurst, 2016; Moyers and Zhang, 2015; Wolfe, 2004). This is particularly problematic when studying protein properties such as length, evolutionary rate, and degree of conserved structure, because these properties themselves directly impact our ability to detect sequence similarity. Simulation studies have shown that these biases in homology detection have the potential to drive spurious phylostratigraphic trends (Elhaik et al., 2006; Moyers and Zhang, 2015; Moyers and Zhang, 2016). However, the number of ‘true’ young genes identified by phylostratigraphy studies may dwarf the numbers of genes that are falsely identified as young due to errors in homology detection, making that bias insufficient to explain previously observed trends (Domazet-Lošo et al., 2017). Statistically correcting for length and evolutionary rate did not remove observed trends in ISD (Wilson et al., 2017) or clustering (Foy et al., 2019).

Simulation studies suggest that reducing false negatives will reduce the strength of trends driven by homology detection bias (Moyers and Zhang, 2018). Therefore, if we improve the sensitivity of homology detection methods, this is predicted to result in weaker trends in protein properties with age if previously reported trends in protein properties, such as disorder and clustering, were spurious. In contrast, if previously reported trends are not artifacts, we predict that with increased sensitivity they will get stronger. This is because many false negatives may be random errors, which reduce power in detecting trends (McLysaght and Hurst, 2016), rather than systematic errors, which can create spurious trends.

Whatever the issues with homology detection, abandoning it is not a viable option. Because homologs share an evolutionary history, they are not statistically independent. But many –omic studies treat genes as independent datapoints and are thus flawed due to pseudoreplication. Phylostratigraphy can solve this problem by making each datapoint represent the independent evolutionary origin of a protein sequence, avoiding phylogenetic confounding (Felsenstein, 1985; Thornton and DeSalle, 2000). It is therefore critical that we improve upon the BLAST-based homology detection methods used within phylostratigraphy, rather than abandon the effort.

Fortunately, more recently developed methods, such as PSI-BLAST (Altschul and Koonin, 1998) and HMMER3 (Finn et al., 2011), are more sensitive to distant homology, while still being computationally efficient. Both are multiple sequence methods, which have long been known to be more sensitive to distant homologs than pairwise methods (Park et al., 1998). PSI-BLAST is a heuristic method that constructs a position-specific score matrix from the alignment of known homologs, which can then be used to iteratively search a protein database for additional homologs. HMMER instead constructs a profile HMM – a probabilistic model of the alignment – and is thus based on a formal statistical framework, yielding superior treatment of indels (Eddy, 2011). Profile HMM based methods perform significantly better than sequence based methods in comparative studies on diverse types of sequence data, with little difference in the rate of false positives (Freyhult et al., 2010; Madera and Gough, 2002).

Unfortunately, because both methods are based on iterative searches, both are vulnerable to model corruption, where a single false positive hit to a non-homologous sequence will pull in many more non-homologous sequences, potentially snowballing over multiple iterations (Pearson et al., 2017). Automated pipelines for whole-genome analysis are technically challenging, and HMMER is generally used with manual supervision of alignments.

Here we increase the sensitivity of homology detection during phylostratigraphy by using the manually curated pfam database, which was constructed using a HMMER3 pipeline. For each pfam, an HMM was built from high quality seed alignments, homologs were then found, and matches above a curated threshold were realigned back to the HMM (El-Gebali et al., 2019).

Pfams are intended to correspond to protein domains, which are structural units, capable of folding independently (Holm and Sander, 1994; for a review discussing domain definitions and identification, see Ponting and Russell, 2002). These domains are often considered to be the ‘true’ units of homology, with full proteins made up of one or more domains in a variety of different combinations (Bagowski et al., 2010; Chothia et al., 2003; Koonin et al., 2000; Moore et al., 2008). The pfam database has developed over time and now curates sets of clearly homologous sequences that are no longer required to correspond to compact folds. We continue, however, to use the term ‘domain’ in conjunction with pfams.

Here we apply HMMER3-based phylostratigraphy to a broad range of taxa. A key methodological innovation is to use pfam domains rather than BLASTp as our unit of homology. A second innovation is that instead of studying only sequences that are present in one or a few focal species, and using other species only to assign age of origin, we make full use of 435 fully sequenced species, effectively treating most of them as focal species. This not only increases our power, but also lets us evaluate how general trends are across different phylogenetic groups. We use only those pfams that are present in at least two species, a quality filter that ensures the exclusion of non-genic contaminants from our data set. We also analyze patterns of protein evolution in full genes, dating these by the oldest pfams that they contain, making results comparable to earlier phylostratigraphies based on the oldest short segment BLASTp hits (Albà and Castresana, 2005; Domazet-Loso et al., 2007). Additionally, to avoid spurious trends most likely to be due to homology detection bias, we avoid the protein metrics that most strongly affect homology detection (i.e. evolutionary rate and length), and focus on trends in properties expected a priori to have less effect on homology detection, such as the hydrophobicity and the degree of clustering of hydrophobic amino acid residues.

Our large data set gives us unprecedented power to answer two key questions. First, which trends in long-term evolution are supported rather than being due to artifacts? Here our increased power enables us to detect finer-scale trends, such as changes in the frequencies of each of the 20 amino acids, rather than simply summary statistics such as ISD. Trends in ISD and clustering get stronger with improved methods, consistent with being biological rather than due to artifacts. Second, are trends in long-term evolution universal across our tree, or specific to a particular group? To answer this, we compare trends among pfams that were born in animals, trends among pfams that were born in plants, and trends among pfams that were born prior to the origin of eukaryotes. We find that while trends in clustering are general across our data set, ISD and amino acid composition trends are specific to particular lineages – the latter compatible with specific biology rather than general artifacts of homology detection.

Results

In this study, we analyzed all protein-coding genes containing at least one pfam annotated in the genomes of 435 species, with a total of 8209 pfams identified (see Materials and methods). We used a variety of quality filters to exclude pfams that may be due to contaminants, annotation errors or horizontal gene transfer (see Materials and methods). Eukaryotic species were included in our data set if they were marked ‘Complete’ by GOLD (Mukherjee et al., 2019), and also present in TimeTree (Hedges et al., 2006), a phylogenetic database that we used to assign evolutionary ages to pfams. A phylogenetic tree summarizing the species used in this analysis is given in Figure 1—figure supplement 1.

Where not otherwise indicated, throughout our main results we focus on non-transmembrane pfams. It is a priori likely that transmembrane pfams experience different selective pressures as they evolve, given the very different biophysical characteristics of a lipid bilayer vs. the cytosol. Pooling is therefore not appropriate unless similarities have first been established.

To avoid pseudoreplication caused by phylogenetic confounding, we take the average of each protein property across all homologs of a pfam, and then treat each homologous set as a single data point in subsequent analysis (see Materials and methods). Phylostratigraphy assigns each such set of homologous pfams to an age class, dated using timetree. We then use linear regression to obtain slope, which provides our estimate of the effect size of the relationship between sequence properties and age. In order to compare our results with those of previous studies, we also conducted our analyses on full genes.

Trends in ISD

Young genes (Willis and Masel, 2018; Wilson et al., 2017; Foy et al., 2019; Mukherjee et al., 2015) and domains (Bornberg-Bauer and Albà, 2013; Buljan and Bateman, 2009; Ekman and Elofsson, 2010; Moore and Bornberg-Bauer, 2012) have been reported to have high ISD, although some have claimed this depends on taxon (Vakirlis et al., 2018). We use IUPred to estimate ISD from amino acid sequence alone, allowing estimation across large-scale genomic data sets that contain many sequences of undetermined structure (Mészáros et al., 2018). We confirm that across our entire data set, ISD is higher in young pfams (Figure 1, linear model: R2 = 0.13, p = 3 × 10−245) and genes (Figure 1—figure supplement 2A, linear model for all genes: R2 = 0.057, p = 1 × 10−229).

Figure 1 with 4 supplements see all
Download asset Open asset
Young domains have high intrinsic structural disorder (ISD), but this trend is driven exclusively by recent animal domains.

Results are for non-transmembrane pfams. (A) The brown linear regression was calculated for recent animal pfams (slope = −0.062, R2 = 0.0097, p = 6 × 10−7, n (number of pfams) = 2456), green for recent plant pfams (slope = 0.056, p = 0.1, n = 183), and black line over ancient pfams in all lineages (slope = −4.2 × 10−4, p = 0.9, n = 3102). Data underlying animal and plant curves is visualized separately in Figure 1—figure supplement 3. Slopes represent the decrease in average IUPred2 score, that is, the predicted propensity of the average amino acid to be disordered, per billion years. Each data point consists of the average across all instances of homologous pfams, across all species in which it occurs. Phylostratigraphy assigns these to age classes, dated using timetree. To help visualize the data, every age class is represented in the figure by a weighted box plot, where the width of the plot indicates the number of pfams in that age class. The median is shown in red, with the boxes representing upper and lower quartiles (the 75th and 25th percentile), and the whiskers indicating 9 and 91 quantiles. For age classes with only a single pfam, values are presented as small red dots. For clarity of presentation our plots do not show outliers, although we note that these are included in our linear regression models. (B) Phylostratigraphy slopes for pfams calculated over different subsets of the data are plotted with their 95% confidence intervals. The point colors correspond to the regression slopes in (A).

Our improved methodology increased the steepness (i.e. effect size) of the relationship between gene ISD and gene age in mouse genes above that previously estimated by Foy et al., 2019 (Figure 1—figure supplement 2B), from −0.028 to −0.050; note that mouse genes have a steeper slope than our results across all taxa. However, we note that age explains a relatively small proportion of the variance in ISD, with the remainder presumably driven by a combination of function, random biological variation, and measurement error (IUPred being an imperfect proxy for ISD or whatever other correlated biophysical property truly underlies the trend).

The slope of the relationship between ISD and age is 1.8-fold steeper for pfams (Figure 1A) than for full genes (Figure 1—figure supplement 2A); this is unsurprising, because whole genes are made up of combinations of sequences of different ages. The fact that protein domains are a more fundamental evolutionary unit than genes (Bornberg-Bauer et al., 2005; Moore et al., 2008), which is reflected in stronger relationships with pfam age, demonstrates an advantage of using a pfam-based phylostratigraphy approach. While this study focuses on pfam domains, we note that results are similar, but with smaller effect sizes, when calculated over genes.

To investigate how the ISD trend varies over time, we recalculated phylostratigraphy slopes over age-restricted subsets of our data. We also compared animal vs. plant lineages, given that these are the two kingdoms for which we have the most data (343 animal genomes and 87 plant genomes). The trend in ISD is not consistent (Figure 1B, Figure 1—figure supplement 3), but is instead driven by recent animal pfams (in which we include all pfams that emerged after the divergence of animals/fungi from plants, 1496 MYA, that are now present in animals).

There is no significant change in ISD over ‘ancient’ pfams (those that emerged prior to the last eukaryotic common ancestor [LECA], which is estimated to have existed around 2101 MYA). This does not mean that ISD is static; for example, the same domain in more exquisitely adapted vertebrates has higher ISD (Weibel, 2020).

There is also no change in ISD with age over recent plant pfams (i.e. all pfams found in plants that are younger than 1496 MY old). Note that we have relatively few plant-specific pfams in our data set compared to animal-specific pfams (Figure 1—figure supplement 3). This is likely due to differences in genome quality and annotation between the groups, and the corresponding lack of available plant genomes that meet our quality standards. However, it is also possible that these numbers reflect a biological reality that animal proteomes are richer in domains (Ramírez-Sánchez et al., 2016). We also note that animal- and plant-specific pfams have different age distributions due to differences in the respective phylogenies, and thus the set of possible divergence times in these groups.

Recent animal pfams have higher mean disorder (0.36) than recent plant pfams (0.26) (Welch’s t-test p = 4 × 10−5, figure 1a shows how this result depends on age within the recent pfam categories), with the latter still higher than ancient pfams (0.21). The difference between animals and plants does not reflect differences in the birth process alone; even in ancient pfams that are shared by plants and animals, mean ISD is 0.27 in animals vs. 0.24 in plants, a smaller but still significant difference (Wilcoxon signed rank test [a paired, non-parametric test] on difference p = 0.005). These results are compatible with animal domains experiencing more selection for high ISD than plants; plant domains have consistently less disorder with the difference being more pronounced in young, lineage-specific domains. This may be because plants produce aggregation inhibiting molecules (Velander et al., 2017), thus reducing the risk that protein aggregation poses to plant cells (see Discussion).

Trends in amino acid frequencies

IUPred scoring of ISD primarily reflects amino acid composition, with hydrophilicity being a major determinant of disorder (Dosztányi et al., 2005b; Wilson et al., 2017). Our larger data set has sufficient power to investigate age trends in the frequencies of each of the 20 amino acids individually, rather than just trends in the single IUPred summary statistic. This can reveal which amino acids drive our ISD results, as well as other potentially interesting patterns in amino acid occurrence.

Trends in amino acid frequencies with age among ancient pfams are essentially identical whether they are assessed using only plant data or only animal data (Figure 2A, Spearman’s ρ = 0.94, p = 6 × 10−6, unweighted Pearson’s R = 0.98, p = 8 × 10−14). In subsequent analyses beyond Figure 2, we therefore pool data across all species in our data set to calculate the phylostratigraphy slopes among ancient pfams with higher resolution.

Figure 2 with 1 supplement see all
Download asset Open asset
Trends in amino acid frequencies as a function of age differ across lineages.

Results are shown for non-transmembrane pfams. Phylostratigraphy slopes are in units of the change in percentage points of an amino acid per billion years. ‘Ancient’ refers to pfams older than 2101 MY, assessed in only plant or only animal instances, whereas ‘plant’ and ‘animal’ slopes are calculated over pfams appearing after the divergence between the animal/fungi and plant lineages, 1496 MY, again assessed in only plant or only animal instances. Lines indicate the standard errors on the slopes. Points are color-coded by their hydrophobicity, as measured by 1-mean relative solvent accessibility (RSA) (Tien et al., 2013), such that buried, hydrophobic amino acids are dark purple, and exposed, hydrophilic amino acids are yellow. RSA scores are based on the water-accessible surface area of amino acids in a training set of proteins of known structure. Amino acid shapes indicate whether they are essential in animals. In (A), y=x is shown as a dashed line, in all other plots dashed lines for x=0 and y=0 are shown for clarity. In (A), Spearman’s ρ = 0.94, p = 6 × 10−6. Spearman correlations for (B), (C), and (D) are not significant (p = 0.8, 0.8, and 0.2, respectively).

In contrast, recent animal pfams and recent plant pfams show different trends in amino acid frequencies with age, with slopes of similar magnitude (Welch’s t-test on absolute slope values p = 0.3), but no correlation in value (Figure 2B, p = 0.8). Recent frequency trends are mostly unrelated to ancient trends, with no relationship for animals (Figure 2C, p = 0.8), and a weak correlation for plants that is not significant (Figure 2D, p = 0.2) and is entirely driven by cysteine.

In plants, cysteine has the steepest phylostratigraphy slope. Miseta and Csutora, 2000 previously reported that %cysteine content decreased with ‘complexity’, from mammals, to plants, to a set of single celled organisms that included both photosynthesizing and non-photosynthesizing bacteria and archaea. Over all pfams, our results agree with the findings of Miseta and Csutora, in that cysteine is slightly more common in animals (2.5%) than in plants (2.2%). But surprisingly, young, plant-specific pfams have 3.5% cysteine, as opposed to 2.6% in animal-specific pfams, and 2.1% in ancient pfams.

Photosynthesis likely creates greater oxidative stress in plants than in animals, which is likely to disrupt disulfide bonds. But in the training data used by IUPred, cysteine residues are mostly in stable disulfide bonds (Mészáros et al., 2018), so IUPred scores cysteines as highly order promoting. If unpaired cysteines are common in plants, this might obscure a trend in ISD. However, recalculating IUPred scores from sequences from which we had excised cysteine resulted in very little change in our ISD results (Figure 1—figure supplement 4).

In young animal pfams, the two amino acids that are most enriched are serine and proline. These are two of the amino acids previously identified as most responsible for differences between eukaryotes and prokaryotes, due in particular to the greater quantity of ‘linker’ regions (protein sequences between pfams) of eukaryotic proteins (Basile et al., 2019). Our results suggest instead that there is nothing special about linker regions. Linker regions might simply be young sequences that happen not to have been annotated as a pfam domain. We find that not just linkers, but also young annotated pfam domains in animals have a high percentage composition of serine and proline. This is perhaps unsurprising; the simple categorization of sequences as either ‘pfam’ or ‘linker’ in genome annotation is artificial and might introduce a variety of biases.

Because animal species with more tyrosine kinases have less tyrosine, it has been argued that selection against deleterious tyrosine phosphorylation has driven a decline in tyrosine content in metazoa (Tan et al., 2009). However, evidence that selection drives tyrosine loss is limited (Pandya et al., 2015), and it has been argued that changes in GC content might explain trends in tyrosine (Su et al., 2011, although see Tan et al., 2011). Examining this question with a phylostratigraphy approach for the first time, we do not see evidence for tyrosine loss in animals; tyrosine has a relatively shallow phylostratigraphy slope that is indistinguishable from 0 (Supplementary file 1). This agrees with the findings of Pandya et al., 2015, who found no significant difference between the frequency spectra of alleles that removed or created tyrosine, suggesting that tyrosine is neither strongly selected for nor against, with trends in tyrosine being modest compared with those of other amino acids. There is, however, more tyrosine in younger ancient pfams than in the oldest ancient pfams, compatible with an ancient rather than recent process of tyrosine loss (see Supplementary file 1).

To assess whether animal-specific trends also reflect amino acid availability, we examined essential amino acids (i.e. amino acids that plants but not animals can synthesize [Guedes et al., 2011]), expecting them to be rarer in young animal domains. While this is the case (with six of nine of their phylostratigraphy slopes being positive; Figure 2), the difference is not significant (Binomial test with 50% chance of a positive slope, one-tailed p = 0.3).

Instead, the ISD score returned by IUPred seems to best summarize animal-specific trends in overall sequence hydrophobicity. We do not find differences among amino acids in hydrophobicity scores such as relative solvent accessibility (RSA) (Tien et al., 2013). While we do see a tendency for hydrophobic amino acids to have positive phylostratigraphy slopes among animal pfams and negative slopes among plant pfams (amino acids are colored by RSA in Figures 2 and 3), this correlation between RSA and phylostratigraphy slope is not significant for any of the lineage- or phylostrata-specific subsets of the data set shown in Figure 2 (see Figure 2—figure supplement 1).

Figure 3 with 2 supplements see all
Download asset Open asset
Ancient domains exhibit similar amino acid trends, whether transmembrane or non-transmembrane.

Phylostratigraphy slopes are in units of percentage point change in composition per billion years. Taxonomic and temporal sub sets of the data are as the same as in Figure 2. Lines indicate the standard errors on the slopes. Points are color-coded by their hydrophobicity, as measured by 1-relative solvent accessibility (RSA) (Tien et al., 2013), such that buried, hydrophobic amino acids are dark blue, and exposed, hydrophilic amino acids are yellow. In all plots, y=x is shown as a dashed line. Congruence is strong for ancient pfams in (A) (Spearman’s ρ = 0.83, p = 5 × 10−6), weakly significant for animal-specific pfams in (B) (p = 0.06) where transmembrane trends tend to be weaker and leucine is an outlier, and not detectable in our weakly powered plant-specific data set in (C) (p = 0.3).

What is more, if hydrophobicity was the main determinant of phylostratigraphy slopes, we would expect amino acid composition to evolve differently in a hydrophobic membrane environment than in the cytosol. However, amino acid slopes are highly correlated between transmembrane and non-transmembrane ancient pfams (Figure 3A). Clearly, trends among ancient domains are not primarily driven by hydrophobicity. Amino acid slopes are also weakly, albeit not significantly, correlated between transmembrane and non-transmembrane pfams in animals (Figure 3B). Breaking this animal correlation is leucine, a clear outlier in that it is enriched in young transmembrane pfams but depleted in young non-transmembrane pfams. Leucine also shows the same pattern to some degree among ancient pfams. There is very little power to detect a correlation in plants, should it exist (Figure 3C).

There has been concern that phylostratigraphy slopes could be driven by homology detection bias (Moyers and Zhang, 2015, Moyers and Zhang, 2016; Wilson et al., 2017). In this case, amino acids that are more changeable, as assessed by the changeability scores of Tourasse and Li, 2000, making homology more difficult to detect, should be over-represented in young genes. However, we find if anything the opposite: more changeable amino acids are (slightly) enriched rather than depleted among the oldest of ancient pfams (Figure 3—figure supplement 1, Spearman’s ρ = 0.50, p = 0.02 for nontransmembrane pfams, and ρ = 0.40, p = 0.08 for transmembrane ancient pfams). Our results are therefore not driven by homology detection bias even at extremely long timescales. There is no correlation between recent animal (p = 0.57) or recent plant (p = 0.7) phylostratigraphy slopes and amino acid changeability, and therefore no evidence that these results are affected by homology detection bias either. Furthermore, we note that homology detection bias is expected to create similar patterns for all taxa. The fact that the strength and even direction of amino acid trends can be different for different taxa is itself also evidence against a strong role for homology detection bias.

Given the striking similarity in ancient trends between transmembrane and non-transmembrane domains, we note that Jordan et al., 2005 claim ongoing trends related to the origin of the genetic code. Specifically, they claim that the amino acids incorporated earliest into the genetic code (Trifonov, 2000) tend to be lost, while more recently acquired amino acids tend to be gained. Jordan et al., 2005 used a parsimony-based method of ancestral sequence reconstruction to detect trends in the amino acid composition of existing protein-coding sequences over time. The method has, however, been much criticized (Goldstein and Pollock, 2006; Hurst et al., 2006; McDonald, 2006).

Our amino acid phylostratigraphy slopes do not correlate with the flux values of Jordan et al., 2005 (Figure 3—figure supplement 2). However, we find that among ancient domains, the very oldest domains are enriched in the amino acids that were hypothesized to have been recruited into the genome first, both in nontransmembrane domains (Spearman’s ρ = −0.62, p = 0.004), and in transmembrane domains (Spearman’s ρ = −0.55, p = 0.01). This suggests that amino acid availabilities at different stages of life could have affected the composition of proteins born at those stages, in an effect that persists to this day. Specifically, it suggests that amino acids recruited later into the genetic code increased in abundance only slowly, creating long-term trends in amino acid availability among different ancient origin dates.

The order of recruitment ranking is taken from Trifonov, 2000, who estimated a possible consensus chronology from the expectations of 40 amino acid ranking criteria, including such diverse factors as results from experiments into primordial conditions, amino acid complexity, and thermostability. However, 3 of the 40 criteria were based on the amino acid compositions of protein assemblages. Therefore, to ensure non-circularity, we recalculated the consensus order, excluding these criteria and calculating the mean rank, as in Trifonov, 2000 prior to their proposed second step of smoothing by a filtering procedure. This recalculation had very little effect on the order of amino acids, and our correlation in ancient nontransmembrane pfams remains significant (Figure 4A, Spearman’s ρ = −0.58, p = 0.008), although the correlation for ancient transmembrane pfams is no longer on its own robust to multiple testing (Figure 4B, Spearman’s ρ = −0.47, p = 0.04).

Figure 4 with 2 supplements see all
Download asset Open asset
Ancient amino acid phylostratigraphy slopes reflect the order of recruitment of amino acids into the genetic code.

Phylostratigraphy slopes for non-transmembrane (A) and transmembrane (B) pfams are in units of percentage points of composition per billion years, with lines indicating the standard errors on the slopes. Phylostratigraphy slopes for ancient pfams are calculated over all lineages, and include all pfams over 2101 MY. Consensus order is modified from Trifonov, 2000 (see description in Results) and is given as rank data, such that the first amino acid to be recruited is given a rank of 1. Regression slopes are shown as dashed lines. Late-recruited amino acids are rare in the most ancient non-transmembrane pfams (A) (Spearman’s ρ = −0.58, p = 0.008) and transmembrane pfams (B) (Spearman’s ρ = −0.47, p = 0.04). Points are color-coded by their metabolic costliness to produce (measured as the number of high energy phosphate bonds required for synthesis, plus the energy lost due to the precursors used in the synthesis), as estimated for aerobic conditions in yeast (Raiford et al., 2008), such that the most costly amino acids are blue, and the least costly are yellow. (C) The correlation (Spearman’s ρ) between amino acid phylostratigraphy slope and order of recruitment over different subsets of our data set. X-axis labels indicate the minimum age of pfams included, in billion years. Lines are the standard errors of the Spearman’s ρ values, calculated using the Fisher transformation (Fisher, 1915).

The order of recruitment ranking is highly speculative, and correlations with it could be driven by a single correlated amino acid physiochemical property. Specifically, amino acids with very simple structures (A, G, and V) have steeper positive phylostratigraphy slopes than any other amino acids, and this alone may drive our observed correlation between phylostratigraphy slopes and order of recruitment. We note that order of recruitment is not significantly correlated with phylostratigraphy slopes in more recent lineages (Figure 4—figure supplement 1). This is in agreement with the broader hypothesis of Jordan et al., 2005 that the order of recruitment shaped early amino acid availabilities and hence protein compositions, but is hard to reconcile with alternative explanations in which amino acid simplicity is directly causally responsible.

Our interpretation of these results is that amino acids that were introduced late into the genetic code remained relatively rare for a prolonged period of time, well after the genetic code was complete. To determine for how long, we ask when the correlation between phylostratigraphy slope and order of recruitment is strongest (Figure 4C). The strongest correlation is with slope over the three oldest phylostrata in our data set. These span the pfams found in last universal common ancestor (LUCA), through to pfams present in all extant eukaryotic but no prokaryotic lineages. The strength and significance of the correlation between order of recruitment and phylostratigraphy slope decrease if younger phylostrata are included, suggesting that from the LECA onward, the order of recruitment into the genetic code no longer shaped amino acid availability.

We do not believe our order of recruitment results to be a byproduct of a correlation with amino acid costliness. This is because there is no more than marginal significance for correlations between phylostratigraphy slope and the metabolic cost of amino acid production in yeast (Figure 4—figure supplement 2; amino acids are color-coded by costliness in Figure 4; Raiford et al., 2008).

Trends in clustering

Finally, we consider whether there are any trends in amino acid order. Specifically, a temporal trend in the value of a ‘clustering’ metric has previously been reported for mouse gene families (Foy et al., 2019). This clustering metric calculates the degree to which hydrophobic amino acids tend to lie close together along the primary sequence, normalized to ensure that values do not depend on length or amino acid frequencies (Irbäck et al., 1996) (see Materials and methods for a full technical description). Not only is there no causal link by which amino acid composition affects clustering, there is also no significant correlation (Pearson’s p = 0.06 across non-transmembrane pfams) We confirm, using our considerably larger data set, that young pfams (Figure 5A, slope = −0.037, R2 = 0.028, p = 1 × 10−50), and young proteins (Figure 5—figure supplement 1A; slope = −0.031, R2 = 0.033, p = 5 × 10−132) have more clustering of their hydrophobic amino acids. In Figure 5—figure supplement 1B, we show that our improved methods for assigning gene age result in a steeper slope, that is, a larger effect size, for mouse genes than that reported by Foy et al., 2019: −0.056 instead of −0.045. We note that while the proportion of variance in clustering that is explained by age is small, this is driven at least in part by large denominator variance due by domain function or random variation.

Figure 5 with 2 supplements see all
Download asset Open asset
Young domains have more clustered hydrophobic amino acids (A), and the trend in clustering with age is consistent across time and taxonomic groups (B).

n (number of pfams) = 8002, 3100, 2456, and 183 for all, ancient, animal, and plant groups, respectively. Clustering has an expected value of 1 for randomly distributed amino acids. Results are shown for non-transmembrane pfams. (A) Each data point consists of the average across all homologous instances of a pfam, across all species in which the pfam occurs. Phylostratigraphy assigns these to age classes, dated using timetree. To visualize our results, every age class is represented by a weighted box plot, where the width of the plot indicates the number of pfams in that age class. The median is shown in red, with the boxes representing upper and lower quartiles (the 75th and 25th percentile), and the whiskers indicating 9 and 91 quantiles. For age classes with only a single pfam, values are presented as red dots. For clarity of presentation our plots do not show outliers, although these are included in our linear regression models. The blue line is the linear regression slope calculated over all lineages. Slope = −0.039, R2 = 0.030, p = 3 × 10−54. Slope represents the average decrease in clustering score per billion years. (B) Phylostratigraphy slopes (in units of change in clustering per billion years) for pfams are calculated over different phylostrata subsets of the data (Figure 5—figure supplement 2), and shown with their 95% confidence intervals. The left-most point corresponds to the regression slope in (A).

Unlike the trend in ISD, the trend in clustering with age is remarkably consistent across age categories (Figure 5—figure supplement 2, Figure 5B: slope = −0.067, p = 0.0005 for animal-specific domains and slope = −0.025, p = 0.008 for ancient domains). While the confidence intervals increase for subsets that contain less data (such as the set of young plant pfams, for which we have little power), they always overlap the clustering slope calculated over all lineages and all phylostrata. This suggests that trends in clustering may be universal. If new sequences were born with similar clustering values, this would imply that the dispersion of hydrophobic amino acids change at an approximately constant rate over the whole of evolutionary time, with the observed trend resulting from the fact that the change has had more time to take place in older sequences. Under this interpretation, the data is also compatible with a slower contemporary rate of change of clustering in ancient pfams.

Discussion

We have shown that in our data set of eukaryote genomes, there is a universal trend for young genes and pfam domains to have clustered hydrophobic amino acid residues, while old domains and genes have more evenly interspersed hydrophobic amino acids. In contrast, trends in amino acid content depend on which lineage is examined. Only animal trends are driven by young protein-coding sequences having high ISD. Ancient trends, from LUCA to the LECA, are driven by the order of recruitment of amino acids into the genetic code. We have less power to decipher plant trends, although we note a particularly strong excess of cysteine, especially relative to the generally low levels at which abundant cellular cysteine is incorporated into plant proteins as a whole.

The effect sizes of sequence age on amino acid frequencies are small. The question is whether they are biologically important. Kosinski et al., 2020 estimated the marginal effects of substituting one amino acid on the fitness associated with expressing a random peptide in E. coli. Encouragingly, they found that those marginal effects correlate with our animal phylostratigraphy slopes, suggesting that our results are indeed driven by the forces that facilitate gene birth. For example, a single switch of one proline in a 50-mer changes genotype frequency by 3% over the course of one cycle of experimental evolution. We can think of that switch as a two percentage point change in proline frequency, that is, of the order of change we see per billion years. While this fitness effect is small, it nevertheless has the potential to be biologically relevant.

Our analysis relies on homology detection of pfam domains. It is important to consider whether systematic errors in homology detection could bias our results. Our improved hmmer-based homology detection methodology should reduce the false negative rate, and thus the scope for homology detection bias. But both our ISD and our clustering trends are stronger rather than weaker than those previously reported using blastp-based methods (Foy et al., 2019; Wilson et al., 2017). Homology detection bias also cannot explain why trends in ISD or amino acid composition are lineage-specific, nor the absence of correlation with amino acids’ evolutionary changeabilities. Additionally, old sequences are expected to be most affected by homology detection bias (Moyers and Zhang, 2017), but it is more recent animal domains that drive the ISD result. Overall, we find substantial evidence contradicting the suggestion that homology detection bias drives trends.

The rejection of homology detection bias does not imply that distant homology is always successfully detected, but rather that errors in age assignment might be primarily random rather than systematic. Evidence is accumulating that de novo gene birth is not, as once thought, rare (or even non-existent) after the origin of life. A recent analysis exploited synteny to find that most taxonomically restricted genes are not the product of divergence, but have instead emerged de novo (Vakirlis et al., 2020). Another recent study excluded genes with insufficient power to detect ancient homology, and similarly found that evidence for de novo origin persisted in a substantial number of cases (Weisman et al., 2020). To get correct trends given random failures in homology detection, all we need is for apparent ages to be correlated with true ages.

Sequences can also be assigned as older than they really are, when they are incorrectly classified as present in distantly related species (false positives), due either to overly permissive homology detection cutoffs, or to contamination (see, for example, Longo et al., 2011; Merchant et al., 2014). In addition, some domains may be present in species due to horizontal gene transfer, which is problematic for the purpose of our analysis as it means that our annotation based on presence in a species does not reflect the evolutionary history of that domain (Liebeskind et al., 2016; Ravenhall et al., 2015). We used a number of filters to remove false positives and likely products of HGT from our analysis (see Materials and methods). False positives are not expected to be subject to systematic error with respect to ISD or clustering (Moyers and Zhang, 2018). In agreement with this, we found that these steps reduced the amount of ‘noise’ in our data set, making the trends we observed stronger with each improved iteration of our quality control.

Our results are consistent with the hypothesis that selection acts to reduce the aggregation propensity of proteins. Polypeptide chains are intrinsically prone to forming aggregates that are toxic to cells (Chiti and Dobson, 2017). Young animal genes avoid aggregation by being on average more disordered and hydrophilic (Wilson et al., 2017), with their few hydrophobic residues tending to cluster together along the protein chain (Foy et al., 2019). Kosinski et al., 2020 directly showed that high ISD tends to make random de novo sequences less harmful when expressed in E. coli. However, in order to be retained, sequences must also consistently serve a function within their host cells; this often requires them to fold, which is promoted by higher hydrophobicity. It is possible that pfam domains eventually reach a stable level of disorder, a balance between selection for folding, selection against misfolding, and mutational pressure. This may be why we do not see trends in ISD in ancient domains.

Frustration between the two aims of selection might be resolved by sequences evolving more sophisticated strategies for folding while avoiding aggregation. One such strategy appears to be the ordering of amino acid residues (Bertram and Masel, 2020; Foy et al., 2019), with older, more evolved sequences having more evenly dispersed hydrophobic amino acid residues. Remarkably, we find that the trend in clustering is common to all lineages (animal, plant, and pre-LECA) included in our analysis, with no evidence for saturation at a stable level of clustering. Hydrophobic amino acids are more evenly dispersed than would be expected from chance in old proteins, ruling out a trend toward a steady state set by mutation (Foy et al., 2019), The trend lowers aggregation propensity compared to the expected value based on their amino acid compositions (Foy et al., 2019), which is compatible with the trend being due to selection having had more time to decrease the clustering of hydrophobic amino acids. If so, then even over vast amounts of evolutionary time, selection continues to reduce levels of clustering in protein domains rather than reaching a steady state. This may be because decreasing clustering is difficult. For a sequence to reduce levels of hydrophobic amino acid clustering while maintaining a similar amino acid composition, many substitutions are needed on a rugged fitness landscape (Bertram and Masel, 2020).

Disordered regions can contain ‘sticky’ amino acids at protein–protein interaction surfaces (Levy et al., 2012), whose hydrophobicity within a hydrophilic region will decrease clustering. More abundant proteins have fewer sticky residues within their disordered regions (Dubreuil et al., 2019). Older domains generally have higher protein abundance (Carvunis et al., 2012), which should cause lower stickiness and higher clustering, not the lower clustering seen in our results. Our clustering trend is therefore not driven by trends in protein ‘stickiness’, which is a possible contributor to aggregation propensity.

Protein length is another hidden factor that might contribute to the trends we see. It is a priori likely that younger proteins are shorter (Van Oss and Carvunis, 2019) with a correspondingly higher surface to volume ratio, and thus a proportionately smaller hydrophobic core and higher ISD score. However, this cannot explain our results, because while younger animal domains do tend to be shorter, shorter animal domains are actually less disordered (R2 = 0.0025, p = 0.007). The reasons for this are unclear; length and disorder are, as expected from surface to volume considerations, negatively correlated in ancient domains and in recent plant domains.

Differences between taxa are surprising, with young plant domains having different trends in amino acid frequencies from young animal domains, and not sharing their preference for high ISD. Plant domains as a whole are more ordered, with ancient domains in plants also having lower ISD in plants than the same domains in animals (although the difference in ISD is most pronounced in lineage-specific domains). We speculate that plants are simply better at coping with proteostasis than other taxonomic groups, using mechanisms other than high ISD to protect themselves from protein misfolding and aggregation. Plant proteomes are rich in predicted amyloidogenic regions (Antonets and Nizhnikov, 2017a), and yet there is a paucity of evidence for aggregate formation in plants under normal conditions (Antonets and Nizhnikov, 2017b). This is perhaps due in part to the presence of molecules such as polyphenols, which inhibit aggregation (Velander et al., 2017). The effectiveness of plant molecules at inhibiting aggregation has also been related to plant longevity (Mohammad-Beigi et al., 2019). Young, ordered plant proteins may therefore be under weaker selection against aggregation than the same protein would be in an animal. Our data set does not have enough plant genomes to give us power to detect whether selection acts on clustering in plants.

The degree of cysteine enrichment in young plant domains is striking, standing out despite the low power of our plant analyses. We hypothesize that young plant domains might be cysteine-rich because cysteine is relatively available in the cell. Producing cysteine is the final step in the assimilation of sulfur by plants. Additionally, cysteine production enables the detoxification of reactive oxygen species produced during photosynthesis and is essential to chloroplast function (Bermúdez et al., 2010). In the cytosol, levels of cysteine appear to play a role in maintaining redox capacity and in pathogen resistance (Romero et al., 2014). Despite the abundance of cysteine in the cytosol, its incorporation into protein sequences, while inexpensive in terms of raw materials, might be deleterious in terms of risk of oxidative damage and unstable disulfide bonds. This might explain why cysteine is progressively lost after high abundance at the time of birth. After cysteine, the two amino acids that are most enriched in young plant pfams are glutamic acid and aspartic acid, E and D. These are the most abundant amino acids in plants as a whole (Kumar et al., 2017) and may therefore also be enriched in young plant domains due to their high availability.

Genomic GC content might determine the amino acid composition of young genes, causing high %GC taxa to spawn more hydrophilic and thus disordered young genes (Ángyán et al., 2012). Our animal vs. plant data have 50% and 47% GC, respectively (Welch’s t-test, p = 0.0001), which could account for several percentage points of difference in ISD (Ángyán et al., 2012), but not for the much larger difference in young domain ISD seen in Figure 1A. We note that GC content evolves quickly between species and is variable on a finer taxonomic scale than is considered in this study. For example, monocot and dicot plants have a greater average difference in GC content than do animals and plants (Li and Du, 2014; Romiguier et al., 2010; Šmarda et al., 2014). Differences in GC content can also not explain the relative excess of cysteine in young plant domains. From inspection of the genetic code, we expect cysteine to easily accommodate any genomic %GC between 1/3 and 2/3, and this is confirmed by its prevalence as a function of species %GC (Li et al., 2015).

Our data are consistent with a hypothesis of selection acting in the usual way, that is, that proteins born de novo go on to discover, through mutation, more sophisticated strategies for folding while avoiding aggregation and/or misfolding. A second, not mutually exclusive possibility is that the trends we observe are due to differential retention of pfam domains and genes, rather than descent with modification. Under this hypothesis, newborn sequences are diverse, and sequences with favored properties are lost less often during subsequent long-term evolution. The relative contribution of these two levels of selection could be investigated if quality control could be made stringent enough to obtain quantitatively reliable loss rates and diversification rates for pfam domains.

In conclusion, trends in the evolution of amino acid composition show surprising differences between taxonomic groups and over different spans of time. Amino acid availability at time of de novo birth has a strong effect on domain composition. Even after billions of years of evolution, ancient domains remain influenced by which amino acids were more abundant as the genetic code was being formed. De novo birth and subsequent trends in plants may also be shaped by amino acid availabilities. In young animal domains, amino acid composition trends may be due to selection for high ISD. However, the inferred ancestors of eukaryote protein sequences show a universal trend toward reducing their levels of hydrophobic clustering with age, no matter their amino acid composition, presumably in an attempt to find a balance between folding and misfolding. The only claim in evolutionary biology that is close to comparable in timescale is the increase in body size in some taxa over the 3.5 billion-year history of life (Cope, 1885; Heim et al., 2015; Payne et al., 2009).

Materials and methods

Data compilation

Request a detailed protocol

We compiled a data set of whole-genome sequences with a sequencing status of ‘complete’ from the GOLD database (Mukherjee et al., 2019, accessed August 7, 2018), resulting in a list of 1138 unique species. We then downloaded species from the Ensembl Biomart interface databases (Kinsella et al., 2011; Zerbino et al., 2018, fungi V40, plant V40, metazoan V40, protest v40, and main V93, accessed between July 31, 2018 and September 12, 2018). Of the resulting 306 species, 19 were not found in the GOLD species list and thus were excluded. We next searched the NCBI RefSeq repository (O'Leary et al., 2016, accessed between September 27, 2018 and October 26, 2018) for the remaining GOLD species, excluding archaeal, bacterial, and viral genomes due to lack of phylogenetic resolution. Of the resulting 948 species, 344 had been annotated, were in the GOLD species list, and were not also in Ensembl Biomart.

We further required all species in our data set to be present in TimeTree (Hedges et al., 2006), eliminating 163. Duplicates were removed (seven species), as were species using an alternative coding alphabet (two species), or species for which there were genome annotation quality issues (such that we were able to identify unusually many more pfams from blast searches of intergenic regions than were present in the currently annotated genome, and/or encountered technical problems such as many annotated genes being rich in stop codons) (four species), resulting in a final data set of coding sequences of 455 unique species.

Protist genomes were retained, but used only as outgroups for dating (for a review, see Eme et al., 2014), while full protein properties were calculated for a core set of 435 animal, plant, and fungal species. Protein-coding sequences of the 455 species in the curated list were then downloaded from the previously listed Ensembl repositories from their respective FTP sites, or from NCBI’s FTP RefSeq repository. For a full set of species used in the analysis, see Supplementary file 2. Following genome acquisition, mitochondrial and chloroplast genes were removed from the data set.

Finally, genes often have multiple annotated transcripts/proteins, resulting in non-independent datapoints. We chose a single transcript to represent each gene, specifically the closest homolog to the most closely related sister species in our data set, identified using reciprocal Blastp (Altschul et al., 1990). If a gene failed to produce any hits below an e-value cutoff of 10−3, the longest gene transcript was chosen instead.

Domain annotation

Request a detailed protocol

Ensembl provided Pfam annotations for all protein transcripts, which are based on InterProScan with default parameters. These annotations were manually downloaded from the BioMart web interface and paired with their corresponding protein using Ensembl’s unique transcript identifiers. To make our data set internally consistent, we replicated Ensembl’s annotation methodology by processing all NCBI sequences with InterProScan (Jones et al., 2014, v.5.20–69.0 downloaded June 20, 2018) using default parameters. All sequences without an associated Pfam were discarded from our analyses.

Domain filtering

Request a detailed protocol

Because a well-resolved tree is crucial for correct age assignment, we focus on eukaryotic genes. Contamination is a major problem in gene annotation (Kryukov and Imanishi, 2016; Merchant et al., 2014; Salzberg, 2017), so we filtered our data set for likely contaminants, in addition to horizontally transferred genes. While horizontally transferred genes are not contaminants, their unusual form of inheritance means is likely to produce incorrect estimates of gene ages. The same applies to the unique evolutionary history of organelles.

Out of 17,929 available pfams, we first excluded 4913 pfams that appeared in our eukaryotic data set despite being annotated as occurring in prokaryotes, but not eukaryotes (as annotated in the tree file ‘pfamA_species_tree.txt’, accessed on 17 May, 2019 from https://pfam.xfam.org/).

Next we excluded pfams via keyword search. Of the 5076 pfams annotated as exclusive to eukaryotes, 218 pfams that contained the terms ‘mitochondria’ or ‘chloroplast’ in the interproscan abstract were excluded, leaving 4858 pfams. Of these eukaryote-specific pfams, we further excluded 561 pfams for which there was no interpro abstract available, leaving 4297 pfams out of the 4858.

Of the 7773 pfams annotated as occurring both in eukaryotes and in bacteria, archaea, or viruses, we excluded 2064 that contained the strings ‘viral’, ‘virus’, ‘bacter*‘, ‘capsid’, ‘bacillus’, ‘Pilus’, ‘Pilin’, ‘mitochondria’, or ‘chloroplast’ unless they also contained one of the terms ‘eukary*‘, ‘vertebrat’, ‘fung*‘, ‘metazoa’, ‘plant’, ‘mammal’, ‘insect’, ‘yeast’, ‘human’, ‘all organisms’, ‘human’, ‘antibod’, and ‘immune’, leaving 5709 pfams out of the 7773. Of these ancient pfams shared with prokaryotes, we further excluded 1112 pfams without an available interpro abstract, leaving 4597 pfams.

The 167 pfams that did not have species annotation information available in the pfam database (https://pfam.xfam.org/) were removed from the data set, resulting in an overall total of 8894 pfams.

The phylogenetic distribution of some remaining pfams appeared indistinguishable from chance contamination. To formalize this criterion, we compared the number of losses inferred by Dollo parsimony to that occurring under the hypothesis that all hits are random contamination, performing this test for all pfams reported in less than half the total number of species. For each possible number of species containing the pfam, we simulated a distribution of the inferred number of losses when that number of species was chosen at random. We then used the mean and variance for the resulting distribution to calculate z-scores for the number of inferred losses from the actual data on each pfam, and rejected the pfam if z-scores > −2, a total of 688 pfams. This resulted in a final set of 8209 pfams.

We also remove all genes that contain a pfam excluded as described above.

Age assignment

Pfams

Request a detailed protocol

Pfam domains were assigned a date of evolutionary origin using TimeTree, using the ‘build a timetree’ option. This option returns the mean divergence times calculated over all published molecular time estimates. We assumed each pfam originated halfway between the node of the most recent common ancestor of all species in which the pfam originated, and the node prior to that. To add additional temporal resolution to the oldest domains, we identified any Pfams in our data set likely to have been present in the first eukaryotic common ancestor, defined as those pfams present both in non-protist eukaryotes and in the protist group ‘Excavata’ (species included in our data set are: Naegleria gruberi, Giardia intestinalis, Giardia lamblia, Spironucleus salmonicida, Angomonas deanei, Bodo saltans, Leishmania braziliensis, Leishmania donovani, Leishmania infantum, Leishmania major, Leishmania mexicana, Leishmania panamensis, Leptomonas pyrrhocoris, Leptomonas seymouri, Perkinsela sp., Strigomonas culicis, Trypanosoma brucei, Trypanosoma cruzi, Trypanosoma rangeli, Trichomonas vaginalis), which is considered a possible outgroup to the other eukaryotes, and with an estimated date of 2230 MYA (Hedges et al., 2001; Hedges et al., 2006; Parfrey et al., 2011).

Beyond TimeTree, pfams inferred to have been present in the LUCA (Weiss et al., 2016) were dated as 4090 MY old. For those pfams that are found in eukaryotes and bacteria, but not archaea, and also those pfams that are found in eukaryotes and archaea, but not bacteria, that is, for those that emerged after LUCA but prior to the emergence of eukaryotes, we assigned them the age of 3145 MY, as the halfway point between the emergence of eukaryotes and LUCA. Dates of these ancient events are imprecise. Recent work has also highlighted the difficulty in identifying domains present in LUCA (Berkemer and McGlynn, 2020).

Full genes

Request a detailed protocol

We dated each gene according to the age of its oldest Pfam. Pfams represent highly conserved sequences, so it is unlikely that genes would be categorized as older based on non-Pfam sequences. Note that while different parts of a gene have different ages, here we classify a protein’s age based on the age of its oldest pfam. We excluded genes without any annotated pfams from this analysis, which helps to remove sequences that are not truly genic, but may result in an underestimate of disorder in our data set as a whole if it is the youngest, most disordered genes that are least likely to have annotated domains.

Homology groups

Request a detailed protocol

Homologous sequences share a common evolutionary origin, and therefore treating them as independent datapoints is a form of pseudoreplication (Wilson et al., 2017). To correct for this, an average taken across each Pfam was treated as a single datapoint in our domain analyses.

For genes, this was achieved by grouping according to their oldest Pfam. For genes with multiple Pfams that are all equally the oldest, a cluster analysis was performed. We started by considering the frequency of the concurrence of two Pfams, A and B, P(AB). If P(AB| A or B) ≥ 50% given all protein-coding genes in our data set, then a link was established. Following all pairwise comparisons, Pfams were grouped together using single-link clustering, and each group was assigned a unique gene homology ID. An average across all genes sharing a gene homology ID was then treated as a single datapoint in our whole gene analyses.

The homology group files used in our analyses are available at https://doi.org/10.6084/m9.figshare.12037281.

Metrics

Transmembrane status

Request a detailed protocol

Genes and pfams were assigned as either transmembrane or non-transmembrane using the program tmhmm (Krogh et al., 2001), which predicts the number and position of transmembrane helices within a protein sequence. A gene was designated as transmembrane if it contained over 18 amino acids that are predicted to be in a transmembrane helix (Krogh et al., 2001). A pfam was designated as transmembrane if it overlapped with a predicted transmembrane helix by a minimum of 50% of the pfam length, or if a transmembrane helix overlapped with a pfam by a minimum of 50% of the helix length; 50% was chosen as an arbitrary albeit reasonable cutoff for the designation of transmembrane status.

Disorder

Request a detailed protocol

Disorder predictions were made for each sequence using IUPred 2 (Dosztányi et al., 2005a; Mészáros et al., 2018), which assigns a score between 0 and 1 to each amino acid. Each protein’s ISD score was calculated by averaging over the values of all amino acids. To determine the disorder of each Pfam, we averaged over only the relevant subset. In order to obtain results with interpretable units, we show untransformed results; however, results are qualitatively similar if data is transformed (using a Box-Cox transform with the optimal value of lambda to achieve normality) prior to analysis. We also recalculated IUPred scores after first removing cysteine residues from protein sequences.

Amino acid composition

Request a detailed protocol

The fractional amino acid composition was found by counting the occurrences of each of the 20 amino acids in each protein or pfam and dividing by the length of the sequence. In order to obtain results with easily interpretable units, we show results for the untransformed proportions. Our qualitative results do not change if amino acid fractions are arcsine transformed (the most appropriate transform for proportions [Sokal and Rohlf, 1995]) prior to analysis (results not shown).

Clustering

Request a detailed protocol

To determine clustering scores, we followed Foy et al., 2019 and compared the variance of the hydrophobicity among blocks of s = 6 consecutive amino acids to the mean hydrophobicity, to produce a normalized index of dispersion. If the length L of a protein was not a multiple of 6, we took the average of all possible p = (L modulo 6) frames, truncating the ends appropriately. The six most hydrophobic amino acids Leu, Ile, Val, Phe, Met, and Trp were scored as +1, and the remaining standard amino acids to −1, as in Irbäck et al., 1996 and Foy et al., 2019. Non-standard amino acid residue abbreviations were also scored as +1 or −1, interpreted as follows: B as corresponding to D or N and hence −1; J as corresponding to either I or L and hence +1; Z as corresponding to either E or Q and hence −1; U as selenocysteine and O as pyrrolysine, both scored as −1.

For a sequence of length N = ⌊L/s⌋ in frame 1≤ f≤p, the scores for each block k = 1,..., N/s of six amino acids were summed to σk,f, with the full truncated sequence summing in total to

Mf= ∑k=1Nsσk,f.

The normalized index of dispersion for a particular frame was then calculated as

ψf= sN∑k=1N/s1K (σk,f-sMf/N)2

with the normalization factor

K=sN2-Mf2N2-N 1- sN.

The total normalized index of dispersion was then calculated by taking the average over all possible frames

ψ = 1p∑f=1pψf.

Clustering scores were again analyzed as raw scores for ease of interpretation of slope units; however, results are qualitatively similar if the data are Box-Cox transformed prior to analysis, as for disorder.

Statistical analyses

Request a detailed protocol

All statistical work in this research was performed using R. Our major results throughout the manuscript are based on ‘phylostratigraphy slopes’, which come from simple linear models of the effect of age on the particular protein metric of interest (i.e. ISD, percent amino acid composition, or clustering). Linear models include all homology groups. The function tidy from the broom package was used to convert model outputs to data frames for analysis and plotting, and to calculate 95% confidence intervals for regression slopes (Robinson, 2014). We used the base function lm to perform basic linear models. The lm function was also used to perform weighted least squares regression on variances. Plots were generated in R, using packages ggplot2 and gridExtra. The parametric boxplots were drawn using the quantile function in base R.

Data availability

Request a detailed protocol

All scripts used in this work can be accessed at: https://github.com/MaselLab/ProteinEvolution copy archived at swh:1:rev:1c3b5dcc3dbce35acc1bbbfaafa29fc6398ecee8. Our raw data and homology files used in our analyses are available at https://doi.org/10.6084/m9.figshare.12037281.

Data availability

All scripts used in this work can be accessed at: https://github.com/MaselLab/ProteinEvolution (copy archived at https://archive.softwareheritage.org/swh:1:rev:1c3b5dcc3dbce35acc1bbbfaafa29fc6398ecee8/). Our raw data, and homology files used in our analyses, are available at https://doi.org/10.6084/m9.figshare.12037281.

The following data sets were generated
    1. James JE
    2. Willis SM
    3. Nelson PG
    4. Weibel C
    5. Kosinski LJ
    6. Masel J
    (2020) figshare
    Data from: Universal and taxon-specific trends in protein sequences as a function of age.
    https://doi.org/10.6084/m9.figshare.12037281

References

    1. Holm L
    2. Sander C
    (1994) Parser for protein folding units
    Proteins: Structure, Function, and Genetics 19:256–268.
    https://doi.org/10.1002/prot.340190309
  98. Book
    1. Sokal R
    2. Rohlf FJ
    (1995)
    Biometry
    W. H Freeman.

Decision letter

  1. Christian R Landry
    Reviewing Editor; Université Laval, Canada
  2. Diethard Tautz
    Senior Editor; Max-Planck Institute for Evolutionary Biology, Germany
  3. Taraneh Zarin
    Reviewer; University of Toronto, Canada

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

Acceptance summary:

The authors present a comprehensive analysis of the composition of protein domains across millions of years of evolution to examine whether evolution is directing proteins towards specific regions of the available sequence space. One key finding is a reduction with time of hydrophobic clustering in proteins, consistent with evolution to reduce aggregation propensity as proteins age. The results presented show that although natural selection can change some protein features rapidly, other changes may be more difficult to achieve and can thus continue to be improved over long-term evolution.

Decision letter after peer review:

[Editors’ note: the authors submitted for reconsideration following the decision after peer review. What follows is the decision letter after the first round of review.]

Thank you for submitting your work entitled "Universal and taxon-specific trends in protein sequences as a function of age" for consideration by eLife. Your article has been reviewed by a Senior Editor, a Reviewing Editor, and two reviewers. The following individuals involved in review of your submission have agreed to reveal their identity: Taraneh Zarin (Reviewer #2).

Our decision has been reached after consultation between the reviewers. Based on these discussions and the individual reviews below, we regret to inform you that your work will not be considered further for publication in eLife.

Although the questions addressed in the manuscript are of potential interest to a broad community, the organization of the manuscript and some of the method descriptions make it difficult to know what are the specific questions being asked (Introduction) and how the data is compiled and analyzed (results). Some of the results appear to be significant but with very small effect sizes, which may require a more careful examination and interpretation. Since the study is focusing on protein domains, it would also be useful to define in the introduction what kind of protein domains are being analyzed. For many people, a protein domain has a tertiary structure, which means that by definition it is ordered. It is therefore not sure what a disorder predictor is predicting for such domains.

Reviewer #1:

Prior phylostratigraphic studies have reported correlations between various properties and gene age, where gene age is estimated from the last common ancestor of the species that harbor homologs of that gene. Failure to identify distant homologs is a source of error in gene age estimates. The prevalence of such errors, whether the error introduces systematic bias, and the impact of errors on phylostratigraphic inference is currently a subject of debate. This manuscript seeks to address these issues through the use of Pfam Hidden Markov models, instead of blastp searches, for homolog detection. The manuscript presents results that are relevant to two lines of inquiry:

– Does the use of Pfam HMMs reduce error in gene age estimates via greater sensitivity? Do the age estimates obtained with Pfam HMMs exhibit evidence of systematic bias that could lead to erroneous conclusions?

– What do gene age estimates obtained using Pfam HMMs tell us about the forces that govern the evolution of amino acid composition in proteins of different ages?

Both questions are important. However, there are substantial difficulties with the analysis and interpretation, as presented here, that must be addressed before is possible to assess the quality of the evidence and how well it supports the conclusions with respect to either question.

I) Lack of a formal hypothesis testing framework: In its current form, this manuscript does not address either of the above questions with sufficient formality and rigor to make a convincing case that the evidence supports a particular set of conclusions. To address this, the introduction should (1) summarize prior work, discuss open questions and unresolved controversies; (2) specifically state which hypotheses will be tested in the current study, with a discussion of the testable predictions that flow from these hypotheses; (3) describe how these predictions will be tested; and (4) what steps will be taken to eliminate confounding factors and rule out alternate hypotheses.

For the first question (mitigation of error due to more accurate homology identification), steps 1 and 2 are handled reasonably well. However, testable predictions are mentioned haphazardly in the results and confounding factors are not adequately addressed. The second question (what forces act on the amino acid composition encoded by genes as they age) is not formally stated as a target of inquiry, but introduced almost in passing in the middle of the manuscript.

II) Prior art: This manuscript states that "A key innovation is to use Pfam domains rather than BLASTp as our unit of homology". However, Pfam HMM models have been brought to bear on the question of homology detection and gene age estimates in a prior work that is not cited or discussed in the current manuscript: Jain et al., 2019. The authors will wish to read that article to determine to what extent it duplicates, contradicts or is complementary to the results in the current manuscript. In addition, they may find that the Jain et al. article offers useful ideas for methodological refinements and/or casts new light on how the results of the current manuscript should be interpreted. In particular, Jain et al., deal explicitly with lineage- and family-specific differences in evolutionary rates, topics that are not explored in depth in the current manuscript.

III) Methodology: The methodology in the current manuscript is not sufficiently well described to allow a reviewer to fully assess the results or another scientist to reproduce them. Technical terms are used without definition. Descriptions of how various quantities are actually calculated are lacking. Attention is given to statistical significance (p-values), but not to effect sizes. The presentation relies heavily on summary statistics and derived data, in ways that may obscure trends in the underlying raw data. Insufficient information is given about the data used, and it is not clear what data and supplementary information will be made available to the reader, other than the eight supplementary figures. Specific examples of these problems are given in the detailed comments below, but this is not an exhaustive list.

IV) Demonstrating that the evidence supports the conclusions: The manuscript contains a number of strongly word assertions without sufficient demonstration that these conclusions are in fact supported by the evidence. The logic underlying the assertions is not spelled out and alternate explanations are not ruled out or even discussed. The manuscript tends to conflate the description of an observation with the interpretation of what that observation means.

– In one example of these issues, the manuscript concludes that the use of Pfam HMMs resulted in improved detection accuracy. While this is highly plausible, the evidence for this conclusion is not rigorously demonstrated. The main evidence presented is that phylostratigraphy slopes for mouse genes are steeper when Pfam HMMs are used for homology detection, compared to blastp. However, the manuscript does not demonstrate that steeper slopes are incontrovertible evidence of more accurate homology detection. A comparison of the gene sets obtained using the two methods would provide a more convincing and direct assessment. Is there a decrease in false negatives without an increase in false positives? The properties of sequences in the two sets should also be examined for evidence of systematic bias.

The discussion of this issue exemplifies the tendency to confuse cause and effect. In the Results section, the statement "Our improved methodology increased the steepness of the relationship between gene ISD and gene age…" does not make a clear distinction between the observation (steeper slopes) and the inference (that steeper slopes are due to the new methodology and that the new methodology is an improvement). The inference must be demonstrated, not simply stated.

– In a second example, the argument is made that the observed relationships between amino acid composition and age are unaffected by systematic error in homology detection because different trends are observed in plant and animal protein domains and "homology detection bias is expected to create similar patterns for all taxa". No argument is presented to support the statement that "homology detection bias is expected to create similar patterns for all taxa." While this might hold under some conditions, it is not clear to me why this should be true in taxonomic lineages that differ substantially in GC content or evolutionary rates.

Even assuming that the "similar patterns for all taxa" prediction is valid, only one taxonomic comparison, plants versus animals, is offered as evidence. If subsets of animal species are compared, for example, are dissimilar patterns also observed? Further, could the different patterns in plants and animals be due to issues with the underlying data? The Materials and methods section describes a substantial effort to obtain high quality genomes. (This is one of the strongest sections in the Materials and methods.) Despite, and perhaps because of, this effort, the number of animal genomes exceeds the number of plant genomes by almost a factor of four. Domain discovery, modeling, and annotation in plant genomes lags substantially behind animal genomes, in part because of the relative dearth of proteomic data in plants. The fact that the slope for recent plant domains is positive and not statistically significant (p=0.1), as well as the large variance in plant slopes shown in Figure 1B, all suggest that there may be problems with the plant data that could be responsible for different trends.

The manuscript contains other, similar problems with interpretation of evidence, in addition to the two examples given above.

Reviewer #2:

James et al., comprehensively outline the challenges in phylostratigraphy and homology detection while carefully applying these methods to detect age-dependent trends in protein sequences. Their method improves signal for previously reported trends such as decreased hydrophobicity and increased hydrophobic clustering in young protein sequences. Interestingly, the authors find increased intrinsic structural disorder (ISD) in young animal domains, but not young plant domains. Their method also allows them to gain insights into changes in amino acid frequency with gene age. These results would be of general interest to eLife readers. The following could help clarify the scope of these results:

1) The authors should elaborate on the diversity of the taxa that are included. Although there are an impressive 435 species included in the analysis, there are only 5 non-plant and non-animal species, all of which are fungi. This is understandably due to quality control for data that are included in the analysis, but it's not clear if it's fair to use the term "universal" for the trends observed throughout the paper if they are based (mostly) on plants and animals.

2) The authors should comment on whether or not there are systematic differences in the plant vs. animal species that are included, and how that could affect the results of the study. For example, is the GC content of the included plant genomes a concern? There seem to be some reported differences in GC content of different plant species (monocots vs. dicots) [Kawabe and Miyashita, 2003; Li and Du, 2014, Šmarda et al., 2014] - are these species broadly sampled in this study? If not, this should be clearly stated.

[Editors’ note: further revisions were suggested prior to acceptance, as described below.]

Thank you for sending your article entitled "Universal and taxon-specific trends in protein sequences as a function of age" for peer review at eLife. Your article is being evaluated by Diethard Tautz as the Senior Editor, a Reviewing Editor, and two reviewers.

As you will see, the first reviewer was overall satisfied with the revisions, although some issues remained. The second reviewer raised important points, many of which are still valid for the latest version of the manuscript and concern some of the main conclusions.

Reviewer #2:

Having gone through the appeal letter and the current/revised version of the manuscript, I feel that the authors have addressed the reviewers' concerns, where possible and within the scope of the paper. The overall goals of the study are now more clearly defined, and the elaboration on the methods, data, and alternative explanations that could underlie the observations is helpful.

Reviewer #3:

Trends in summary statistics of proteins between plants and animals based on grouping age of origin show that bias in detecting homologous proteins is accounted for and not the strongest effect in results. Authors convincingly show (1) there are minor differences in degree of clustering of hydrophobic residues between proteins with ancient origins and more recent animal proteins; (2) animal and plant proteins with more recent origins have different degrees of predicted disorder; (3) Amino acid frequency trends are surprisingly consistent between ancient transmembrane and non-transmembrane domains, less surprisingly ancients assessed in plants and animal. Other claims are not convincing.

1) A main claim of the manuscript, that "events during the earliest stages of life continue to have an impact on the composition of ancient sequences" is not supported. The research, as described in the "general assessment" above, is focused on amino acid composition and other protein summary statistics compared among age of origin classes. Even if most of the analysis is accepted as correct, this data can only be interpreted as evidence for change if you assume the starting composition of the proteins in these age classes is the same.

2) No evidence was provided that differences in summary statistics were driven by selection, as opposed to slow relaxation to equilibrium. This should not be implied in the Discussion section. Further, there is no evidence that anything has reached any kind of "optimum". These unfounded claims, for me, detract from the overall analysis.

3) The difference in plant and animal recent disorder prediction measures is a good point to refute some previous analyses, but it contradicts the idea that disorder is necessarily a selection-drive property of new proteins.

4) The slope in excess clustering of hydrophobic amino acids is framed as being consistent, but this is an unconventional way of putting it given that the plant variation is so high, and the mean is in the opposite direction. Further, the ancient proteins are in the same direction, but barely so (just into apparent significance based on direction).

5) The effect size of the change in excess clustering among these groups is not big compared to the magnitude of the excess clustering overall. It was not clear to me that this is enough to matter, particularly over short time periods such as a few hundred million years. It would be good to demonstrate that this average magnitude effect would have a large impact on disorder or stickiness.

6) Because the clustering effect size is so small, I was not convinced that it would be difficult to alter the amino acid composition to this degree if required by selection. That is, my sense is that protein composition is moderately malleable to allow a steady state to be achieved fairly quickly (eg millions of years) such that the effect of origin age would not be detected on this time scale.

7) While the correlation in composition between transmembrane and non-transmembrane AA usage frequencies is impressive, three amino acids, G, A, and V, are strong outliers. These are the simplest three amino acids, and while they were possibly recruited early into the genetic code, the cited inference is fairly speculative, and was itself a prediction based largely on the amino acid simplicity. The "order of recruitment" correlation with slope (Figure 4) consists of 17 amino acids with little correlation, and these early three that have a large positive slope. This does not convince that age of recruitment continues to affect frequencies. There are many possible amino acid property vectors, and clarity on what has driven these differences is lacking.

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

Author response

[Editors’ note: The authors appealed the original decision. What follows is the authors’ response to the first round of review.]

Although the questions addressed in the manuscript are of potential interest to a broad community, the organization of the manuscript and some of the method descriptions make it difficult to know what are the specific questions being asked (Introduction)

We would like to appeal this decision. Reviewer 2 was very positive, while reviewer 1 misunderstood what the manuscript was about, as we detail in this document. We acknowledge weaknesses in our writing of the Introduction, which led reviewer 1 to the incorrect idea that our focus was on establishing the technical superiority of certain methodology, rather than broader biological findings that we study here using methodologies that have elsewhere been shown to be improvements. We have fixed this flaw in the Introduction. However, we believe the reviewer could have resolved their confusion about the central focus of the manuscript if they had referred back to the Title and Abstract, both of which we continue to find quite clear, and entirely incompatible with reviewer 1’s characterization of what the two questions of our paper are.

Guided by a focus on the wrong two questions, reviewer 1 provided detailed comments about peripheral findings (e.g. a great deal of focus on comparisons to previous whole-gene work in mouse found in Figure 1—figure supplement 2 and Figure 5—figure supplement 1) as though they were central, and made no specific comments or mentions of what were actually our central findings.

and how the data is compiled and analyzed (results). Some of the results appear to be significant but with very small effect sizes, which may require a more careful examination and interpretation.

We now more thoroughly discuss the magnitude of the effect sizes in the manuscript and their potential to be biologically relevant, drawing on the results of Kosinski et al., (2020) in order to convert effect sizes into fitness terms. We also note here that all the effects observed in this study are due to factors experienced during the de novo birth of protein sequences, rather than functional specificity or chance variation, and as such we do not expect large signal : noise ratios. We expect function to largely determine protein properties, a point we now make explicitly in the manuscript.

Since the study is focusing on protein domains, it would also be useful to define in the introduction what kind of protein domains are being analyzed. For many people, a protein domain has a tertiary structure, which means that by definition it is ordered. It is therefore not sure what a disorder predictor is predicting for such domains.

This was previously done in the second half of a paragraph in the Introduction, which gives our operational definition of domains as things present in the pfam database. Our use of pfams follows standard practice for high-throughput analysis of domains. We have split that paragraph into two in order to better highlight the definition and added the two sentences: “The pfam database has developed over time, and now curates sets of clearly homologous sequences that are no longer required to correspond to compact folds. We continue, however, to use the term “domain” in conjunction with pfams.”

Reviewer #1:

Prior phylostratigraphic studies have reported correlations between various properties and gene age, where gene age is estimated from the last common ancestor of the species that harbor homologs of that gene. Failure to identify distant homologs is a source of error in gene age estimates. The prevalence of such errors, whether the error introduces systematic bias, and the impact of errors on phylostratigraphic inference is currently a subject of debate. This manuscript seeks to address these issues through the use of Pfam Hidden Markov models, instead of blastp searches, for homolog detection.

The switch from blastp to the hmms used by pfam is indeed one of the ways our manuscript is an advance on previous work. Another is the extensive quality filters that we introduce, e.g. to remove contaminants. A third, required for our major results, is the dramatic increase in the number of focal species relative to past phylostratigraphy studies.

The manuscript presents results that are relevant to two lines of inquiry:

– Does the use of Pfam HMMs reduce error in gene age estimates via greater sensitivity?

This is not a question that our manuscript addresses, and as such was stated as fact, accompanied by citations from the Eddy lab. We take it as well-established prior knowledge that hmmer, when used with a multiple sequence alignment, has dramatically higher sensitivity than blastp performed in a pairwise manner, with no loss of specificity. We have now expanded upon this point in the introduction and have added additional references which compare homology detection methods, to emphasize the demonstrable improvement of HMM-based methods.

It is logically entailed that improving homology detection reduces errors in age estimates. We also cite the simulations of Moyers and Zhang as evidence that improvement in age estimates resulting from greater sensitivity will, as naively expected, reduce artefactual trends rather than strengthening them.

Upon careful re-reading of our submitted manuscript, we cannot understand how it could have been read as having this question as a central line of inquiry, but we nevertheless hope that these edits make this clearer.

Do the age estimates obtained with Pfam HMMs exhibit evidence of systematic bias that could lead to erroneous conclusions?

Something like this is indeed part of our study, albeit merely as a necessary validation step along the way to scientific insight rather than a primary scientific focus. However, we note that effect sizes matter here. Systematic bias much smaller in size than the actual trends would not be a problem, so our question is not, as phrased by reviewer 1, a question of the existence of bias. The question for our purposes is instead whether the trends we see are due to bias.

We have added a citation of related recent work from our lab (Kosinski et al., 2020) that validates our amino acid trends against an external experimental dataset on the fitness effects of peptides, further ruling out an important role for homology detection bias, and supporting our interpretation of phylostratigraphy trends as reflective of the process of gene birth.

– What do gene age estimates obtained using Pfam HMMs tell us about the forces that govern the evolution of amino acid composition in proteins of different ages?

This is not an accurate representation of the related question we address. We ask how protein sequences born at different times are different. Differences we find might be explained by the same forces governing their subsequent evolution but accumulating over different amounts of time, since they were all born with similar properties. Alternatively, they might be due to different conditions of life at the time of birth – as we conclude with respect to our oldest age groups.

We cannot (with the Materials and methods of this paper, this is something we are working on) and do not attempt to identify differences in subsequent evolutionary forces in proteins of different ages. We are uncomfortable in any case with the term “forces” in evolution, which in this case include regression to the mean – this could be seen as a force of mutation bias, but more causally important is the bias in the initial state. Reviewer 1 seems to have misunderstood the central questions and findings of the manuscript.

Both questions are important.

Beyond the details of these questions above, we agree that both the question of whether trends are real biological signals vs. artefacts, and the question of what evolutionary cause they have if real, are important. The reviewer has omitted any mention of the primary question we address (as clearly reflected in both our Title and our Abstract), which is whether patterns are the same vs. different in the 3 groups we address – animals, plants, and ancient protein sequences that predate the divergence of these two groups.

However, there are substantial difficulties with the analysis and interpretation, as presented here, that must be addressed before is possible to assess the quality of the evidence and how well it supports the conclusions with respect to either question.

I) Lack of a formal hypothesis testing framework: In its current form, this manuscript does not address either of the above questions with sufficient formality and rigor to make a convincing case that the evidence supports a particular set of conclusions. To address this, the introduction should (1) summarize prior work, discuss open questions and unresolved controversies; (2) specifically state which hypotheses will be tested in the current study, with a discussion of the testable predictions that flow from these hypotheses; (3) describe how these predictions will be tested; and (4) what steps will be taken to eliminate confounding factors and rule out alternate hypotheses.

This is partly a matter of personal style. We have two main scientific questions rather than one question with a short enough list of hypotheses to be easily listed. Our two questions are not what the reviewer highlighted above, but are as follows: Question #1: what non-artefactual patterns are there in protein properties as a function of age? Question #2: are these patterns universal or taxon-specific? To use the framework preferred by this reviewer, and moreover to do so honestly, we would need to list as a separate hypothesis every conceivable pattern that our methods could uncover. This would not in our view improve clarity, given the large number of possible patterns considered in a combinatorial manner across our 3 taxonomic groups. There are just too many different possible outcomes to coherently list them all.

We agree that the Introduction needed a clearer focus on our two main questions. We have made substantial edits to fix this, in particular to what are now the fifth and last paragraphs, and believe it to now be much clearer.

For the first question (mitigation of error due to more accurate homology identification), steps 1 and 2 are handled reasonably well. However, testable predictions are mentioned haphazardly in the results and confounding factors are not adequately addressed. The second question (what forces act on the amino acid composition encoded by genes as they age) is not formally stated as a target of inquiry, but introduced almost in passing in the middle of the manuscript.

It would have been helpful for the reviewer to tell us where they thought it was introduced in passing in the middle of the manuscript, so we could clarify that this was not a question that the current manuscript is able to address. Phylostratigraphy trends are informative about the forces shaping gene birth, not about the forces changing genes as they age.

We made a number of changes to the Introduction to try to make this clearer earlier. The manuscript pivots at a certain point from confirming previously reported trends in ISD and clustering to discovering new trends in amino acid composition. This was previously introduced extremely briefly in the last paragraph of the Introduction. By splitting that paragraph in two, we have given this more emphasis, to better prepare the reader. We also added material to the fifth paragraph of the Introduction, to introduce broad hypotheses concerning patterns for the phylostratigraphic trends of the 20 amino acids.

II) Prior art: This manuscript states that "A key innovation is to use Pfam domains rather than BLASTp as our unit of homology". However, Pfam HMM models have been brought to bear on the question of homology detection and gene age estimates in a prior work that is not cited or discussed in the current manuscript: Jain et al., 2019. The authors will wish to read that article to determine to what extent it duplicates, contradicts or is complementary to the results in the current manuscript.

We now reference this paper in our manuscript. It primarily focuses on the question of orthology detection which, while also important, is not the main question in our manuscript, which deals with trends in structural (rather than functional, as in Jain et al.,) properties with age. In Jain et al., proteins are taken as the unit of homology, however the authors incorporate pfam information in order to parameterize their evolutionary models. This is considerably different to the way in which pfams were used in the present analysis as the unit of homology.

In addition, they may find that the Jain et al. article offers useful ideas for methodological refinements and/or casts new light on how the results of the current manuscript should be interpreted.

We note that the proposed traceability index in Jain et al. is based on the results of BlastP orthology searches, and so while relevant to older phylostratigraphy results it is not applicable to our manuscript, which use multiple alignment based HMM searches, which are able to more accurately detect even distant homology. Incorporating the Jain et al. methodology would be computationally expensive considering our much larger dataset, and is unlikely to result in an improvement in gene age estimation.

In particular, Jain et al., deal explicitly with lineage- and family-specific differences in evolutionary rates, topics that are not explored in depth in the current manuscript.

We do not address evolutionary rates in our manuscript because we do not study how proteins evolve over time, but only the factors influencing the birth process. See also our more detailed response below to a more detailed suggestion from reviewer 1 as to why we might want to look at evolutionary rates.

III) Methodology: The methodology in the current manuscript is not sufficiently well described to allow a reviewer to fully assess the results or another scientist to reproduce them. Technical terms are used without definition. Descriptions of how various quantities are actually calculated are lacking. Attention is given to statistical significance (p-values), but not to effect sizes.

We now define the abbreviation ISD upon first use – we were unable to determine which other technical terms the reviewer is referring to. We now also fully describe how phylostratigraphy slopes are calculated in both the Results and Materials and methods.

Effect sizes are represented in our paper by regression slopes. These were always visible in the Figures, and we have now increased the degree to which we discuss them in the text. Some of our arguments are based on the direct comparisons of effect sizes.

We consider slope (i.e. signal) to be the effect size rather than R2 (i.e. signal: noise ratio). We have also however also increased our discussion of R2 values. For these, visual representation takes the form of box plots.

The presentation relies heavily on summary statistics and derived data, in ways that may obscure trends in the underlying raw data.

We are not sure what is meant. The only processing we do of raw data prior to fitting the clearly described linear models is to take the average over all instances of a pfam such that a pfam constitutes only a single datapoint to the analysis. To treat them as separate datapoints would introduce a massive amount of pseudoreplication. To use a random effect term, as we did in Wilson et al., 2017 and Foy et al., 2019, is not computationally viable for the much larger number of species analyzed here.

Some confusion here may come from our box plots, which are shown only to enable the visualization of our raw data, and do not enter the statistical pipeline.

Our amino acid slopes do fall into the category of summary statistics / derived data. Using one summary statistic per amino acid is necessary in order to analyze their relationships with other sets of 20 numbers corresponding to other amino acid properties. We believe this illuminates rather than obscures. We do retain confidence intervals / weights for these downstream analyses of amino acid summary statistics.

Insufficient information is given about the data used, and it is not clear what data and supplementary information will be made available to the reader, other than the eight supplementary figures. Specific examples of these problems are given in the detailed comments below, but this is not an exhaustive list.

Extensive raw data tables were available at figshare, as reported in the manuscript.

IV) Demonstrating that the evidence supports the conclusions: The manuscript contains a number of strongly word assertions without sufficient demonstration that these conclusions are in fact supported by the evidence. The logic underlying the assertions is not spelled out and alternate explanations are not ruled out or even discussed. The manuscript tends to conflate the description of an observation with the interpretation of what that observation means.

We respond to more specific criticisms below, however we hope that our improved Introduction somewhat addresses these concerns, which we suspect stem from misunderstanding as to what we believe is known from prior literature vs. what questions the current study addresses.

– In one example of these issues, the manuscript concludes that the use of Pfam HMMs resulted in improved detection accuracy. While this is highly plausible, the evidence for this conclusion is not rigorously demonstrated.

This is explicitly assumed, not concluded, as the revised Introduction now makes clearer, citing published evidence in support.

The main evidence presented is that phylostratigraphy slopes for mouse genes are steeper when Pfam HMMs are used for homology detection, compared to blastp. However, the manuscript does not demonstrate that steeper slopes are incontrovertible evidence of more accurate homology detection.

These results in Figure 1—figure supplement 2 and Figure 5—figure supplement 1 are evidence that the trends are not artefactual. The logic is that when methods are improved, the effect size of the trend (i.e. the slope) gets larger.

A comparison of the gene sets obtained using the two methods would provide a more convincing and direct assessment. Is there a decrease in false negatives without an increase in false positives?

There is no ground truth available for the age of different genes. Without it, there is no assessing the absolute rates of false negatives and false positives. What the reviewer is asking for is logically impossible.

The properties of sequences in the two sets should also be examined for evidence of systematic bias.

Without there being a ground truth available for gene ages, we do not believe that comparisons of gene data sets could provide any evidence of bias.

We also note that the comparison of old vs. new is possible only for genes not pfams. This gene comparison was included in supplementary figures, rather than the main text, because it is not a major focus of our analysis, which is focused on pfams. It is simply one of multiple supporting pieces of evidence that the trends we look at are unlikely to be mere artifacts.

Our sets of genes and pfams are available in the Supplementary data on figshare, for anyone who wants to do their own analysis.

The discussion of this issue exemplifies the tendency to confuse cause and effect. In the Results section, the statement "Our improved methodology increased the steepness of the relationship between gene ISD and gene age…" does not make a clear distinction between the observation (steeper slopes) and the inference (that steeper slopes are due to the new methodology and that the new methodology is an improvement). The inference must be demonstrated, not simply stated.

As discussed above, the improved nature of the methodology is not an inference but a premise. That steeper slopes are obtained with the new methodology is an observation, which supports the inference that the slopes are “real” i.e. have a biological and not artefactual cause.

– In a second example, the argument is made that the observed relationships between amino acid composition and age are unaffected by systematic error in homology detection because different trends are observed in plant and animal protein domains and "homology detection bias is expected to create similar patterns for all taxa". No argument is presented to support the statement that "homology detection bias is expected to create similar patterns for all taxa." While this might hold under some conditions, it is not clear to me why this should be true in taxonomic lineages that differ substantially in GC content or evolutionary rates.

Homology detection bias is defined as what happens when some types of sequences are systematically more prone to undetectable homology, leading to systematic error in which sequences are assigned too young an age. It has the potential to affect our trends when the protein property of interest is correlated with the ability to detect homology.

If one taxonomic group has a faster evolutionary rate, then homology detection becomes harder, and so there is a greater opportunity for homology detection bias within that group. But this does not change whether or not the property of interest is substantially correlated with the ability to detect homology.

Similarly, a %GC that has extreme divergence away from 50% might make homology detection harder. This seems in any case moot, because the taxonomic lineages in question are “plants”, “animals” and “pre-LECA”, which do not systematically differ from one another in this way. We have now added a section to the Discussion to address this point.

We do discuss at length the possibility that more changeable amino acids make it harder to detect homology, which has more potential to be a confounding factor in our work. We note that it is standard practice to use the same BLOSUM substitution matrix, regardless of which taxonomic group is being studied. This shows that the consensus in the field is that the factors that would change amino acid compositions which are most prone to homology detection bias do not vary substantially among taxa.

Even assuming that the "similar patterns for all taxa" prediction is valid, only one taxonomic comparison, plants versus animals, is offered as evidence.

Our third group is inferred pre-LECA ancestors, so in fact there are 3 pairwise comparisons across 3 taxonomic groups, with none of them showing similarity in amino acid trends.

If subsets of animal species are compared, for example, are dissimilar patterns also observed? Further, could the different patterns in plants and animals be due to issues with the underlying data? The Materials and methods section describes a substantial effort to obtain high quality genomes. (This is one of the strongest sections in the Materials and methods.) Despite, and perhaps because of, this effort, the number of animal genomes exceeds the number of plant genomes by almost a factor of four. Domain discovery, modeling, and annotation in plant genomes lags substantially behind animal genomes, in part because of a the relative dearth of proteomic data in plants. The fact that the slope for recent plant domains is positive and not statistically significant (p=0.1), as well as the large variance in plant slopes shown in Figure 1B, all suggest that there may be problems with the plant data that could be responsible for different trends.

We agree with the reviewer, and our manuscript already stressed, that the plant analyses (including but not limited to the ISD slope mentioned above) are underpowered. The same would be true if e.g. we were to compare vertebrates to insects – there is too little data on the latter. We note that the ancient trends, which the reviewer does not comment on, do have sufficient power, and that they can be compared to the animal trends. Note that Figure 2A shows that the ancient trends can be pulled out just as easily from the plant data as from the animal data, suggesting that the problem is with the number of plant-specific protein domains, not the plant genomes per se.

The manuscript contains other, similar problems with interpretation of evidence, in addition to the two examples given above.

While a more detailed list would of course be welcome, we think our responses to the examples above illustrate that the problem seems on a systematic basis to be one of style and communication, in which the purpose of the study was not clear to the reviewer, rather than our errors in the interpretation of evidence (for a purpose that the reviewer thinks we had or should have had, but which we did not). Our edits to the Introduction should have fixed the communication issue, clarifying the purpose of our study.

Reviewer #2:

James et al., comprehensively outline the challenges in phylostratigraphy and homology detection while carefully applying these methods to detect age-dependent trends in protein sequences. Their method improves signal for previously reported trends such as decreased hydrophobicity and increased hydrophobic clustering in young protein sequences. Interestingly, the authors find increased intrinsic structural disorder (ISD) in young animal domains, but not young plant domains. Their method also allows them to gain insights into changes in amino acid frequency with gene age. These results would be of general interest to eLife readers. The following could help clarify the scope of these results:

1) The authors should elaborate on the diversity of the taxa that are included. Although there are an impressive 435 species included in the analysis, there are only 5 non-plant and non-animal species, all of which are fungi. This is understandably due to quality control for data that are included in the analysis, but it's not clear if it's fair to use the term "universal" for the trends observed throughout the paper if they are based (mostly) on plants and animals.

We have added a supplementary figure outlining the phylogeny of species used in this analysis. We note that the ancient pfam domains included in our analysis predate the split of animals/fungi and plants, with some estimated to have been present in LUCA. Our claim to universality is based on being present in 3 lineages: animal, plant, and pre-LECA. We have edited the Abstract and other places to make clearer the fact that all reported trends, including those that predate eukaryotes, are assessed in eukaryotes.

2) The authors should comment on whether or not there are systematic differences in the plant vs. animal species that are included, and how that could affect the results of the study. For example, is the GC content of the included plant genomes a concern? There seem to be some reported differences in GC content of different plant species (monocots vs. dicots) [Kawabe and Miyashita, 2003; Li and Du, 2014, Šmarda et al., 2014] - are these species broadly sampled in this study? If not, this should be clearly stated.

We now discuss possible differences in genomic GC content between the two groups. Thank you for the references, which we include.

[Editors’ note: what follows is the authors’ response to the second round of review.]

Reviewer #3:

Trends in summary statistics of proteins between plants and animals based on grouping age of origin show that bias in detecting homologous proteins is accounted for and not the strongest effect in results. Authors convincingly show (1) there are minor differences in degree of clustering of hydrophobic residues between proteins with ancient origins and more recent animal proteins; (2) animal and plant proteins with more recent origins have different degrees of predicted disorder; (3) Amino acid frequency trends are surprisingly consistent between ancient transmembrane and non-transmembrane domains, less surprisingly ancients assessed in plants and animal. Other claims are not convincing.

We thank the reviewer for their comments and we believe that we have addressed their concerns. Our interpretations of results have been softened throughout the revised manuscript, with particular attention to what the role of selection might or might not have been. Detailed responses are given below.

1) A main claim of the manuscript, that "events during the earliest stages of life continue to have an impact on the composition of ancient sequences" is not supported. The research, as described in the "general assessment" above, is focused on amino acid composition and other protein summary statistics compared among age of origin classes. Even if most of the analysis is accepted as correct, this data can only be interpreted as evidence for change if you assume the starting composition of the proteins in these age classes is the same.

The authors seem to agree that the results "are not evidence of composition change" but "due to sequences born at different times". In the early version I looked at, I interpreted the main claim about "impact from early events" to be arguing for composition change. I think we are in agreement about this now.

We agree with the reviewer that our results are not evidence for composition change; our interpretation (“main claim” above) is that proteins born at different times have different amino acid compositions, due to the differing availability of amino acids at different stages of life. We have made this interpretation clearer in the revised manuscript, both in the abstract and in the Results section.

2) No evidence was provided that differences in summary statistics were driven by selection, as opposed to slow relaxation to equilibrium. This should not be implied in the Discussion section. Further, there is no evidence that anything has reached any kind of "optimum". These unfounded claims, for me, detract from the overall analysis.

The Discussion still states that “Our results are consistent with the hypothesis that selection acts to reduce the aggregation propensity of proteins.” This weak claim remains true, and the remainder of this paragraph cites other work (not the current manuscript) in support. We have been more careful in our discussion of selection, removed mentions of “optima”, and added explicit discussion of the possibility of relaxation to (mutational) equilibrium.

3) The difference in plant and animal recent disorder prediction measures is a good point to refute some previous analyses, but it contradicts the idea that disorder is necessarily a selection-drive property of new proteins.

We do not believe that the non-universality of trends is evidence against selection. Neither do we believe that if we had found trends to be universal, that would on its own have been evidence for selection. We see universality and selection as separate questions. We discuss the possibility that de novo proteins in plants may be under different selection pressures, due to the presence of aggregation inhibiting molecules such as polyphenols in plant cells, and we now include a mention of this in the Results section so as to remove any apparent contradiction.

4) The slope in excess clustering of hydrophobic amino acids is framed as being consistent, but this is an unconventional way of putting it given that the plant variation is so high, and the mean is in the opposite direction. Further, the ancient proteins are in the same direction, but barely so (just into apparent significance based on direction).

This section has been altered to be more intuitive and stressing that we have little power in plant data. However, it is clear that the overall trend in clustering overlaps the trend for ancient domains, which is not the case for the trends in ISD. While the error bars are larger, the trend remains highly significant for ancient domains (see Figure 5—figure supplement 2), not barely so as the reviewer states. We have included in the Results text the slopes and p-values for each of the two significant taxonomic groups in order to assure the reader of the significance of this result, in addition to continuing to display slopes and associated confidence intervals in Figure 5B. We also added values for the number of data points (n) to the Figure 5 legend.

5) The effect size of the change in excess clustering among these groups is not big compared to the magnitude of the excess clustering overall. It was not clear to me that this is enough to matter, particularly over short time periods such as a few hundred million years.

In the previously revised version of the manuscript that the reviewer did not see we highlighted and interpreted the small effect size of our result in the second paragraph of the Discussion.

It would be good to demonstrate that this average magnitude effect would have a large impact on disorder or stickiness.

Metrics of disorder or stickiness are dominated by amino acid composition, which clustering excludes by design. We have added confirmation of the lack of a significant correlation between clustering and disorder. Our interest in clustering come precisely because it may affect a variety of protein properties in a manner that is independent from amino acid frequencies. We have improved our discussion of the clustering metric in the Results section to avoid confusion.

6) Because the clustering effect size is so small, I was not convinced that it would be difficult to alter the amino acid composition to this degree if required by selection.

Similarly, in point #6, the confusion may have to do with communication about ongoing compositional change. In addressing this point though, it would be good to clarify what is meant by "composition". I understood it to be quite general, such that composition change would include the substitutions required to affect clustering. When the authors say clustering is "not affected by amino acid composition", I think they are using it to mean something narrower, equivalent to amino acid frequencies. I don't think we are disagreeing here, but miscommunicating, possibly due to my lack of precision in word choice for which I apologize. I want to see some citation or evidence, though, if the authors are going to argue that it requires "many", but I think we're roughly in agreement here. The authors' point out that there are differences in the proteins that lead to this clustering difference (no matter the number of amino acid substitution required) stands.

We agree this may be a miscommunication. Yes, by “amino acid composition” we mean amino acid frequencies, and we have replaced this term throughout the manuscript. The distinction between amino acid frequencies and order is now clearer. Our clustering metric is designed to be independent from the former, in order to allow focused study of the latter.

Changing clustering requires many substitutions distributed across the entire protein-coding sequence. We have included further discussion of the clustering metric in the revised manuscript. Largely, we refer the reader/reviewer to Bertram and Masel, 2020, a manuscript devoted in large part to this point.

That is, my sense is that protein composition is moderately malleable to allow a steady state to be achieved fairly quickly (eg millions of years) such that the effect of origin age would not be detected on this time scale.

As we do detect systematic trends in both clustering and amino acid composition with protein age, even among ancient domains, we are unsure how to interpret the reviewer’s comment. This is a major reason why we think our findings are so interesting – because they do not confirm the intuition of the reviewer and others. Empirically, we found that clustering has not reached a steady state. The fact that our finding is so unintuitive in being found at all is what makes it remarkable, despite its small effect size. We discuss several interpretative hypotheses as to what could account for this result.

Even if there were a steady state, it would not be one by mutation but by mutation-selection balance: hydrophobic amino acids are more evenly dispersed than would be expected from chance in older proteins. We now fully discuss this point in the seventh paragraph of the Discussion.

7) While the correlation in composition between transmembrane and non-transmembrane AA usage frequencies is impressive, three amino acids, G, A, and V, are strong outliers.

We assume the reviewer is referring to these as outliers in terms of the magnitude of trends in frequencies. In the Figure 3 analysis, we use Spearman’s ranked correlation coefficient, which disregards these magnitudes. In addition, we note these amino acids fall very neatly along the dashed x=y line shown in Figure 3A.

These are the simplest three amino acids, and while they were possibly recruited early into the genetic code, the cited inference is fairly speculative, and was itself a prediction based largely on the amino acid simplicity.

The inferred order of recruitment is indeed speculative by its nature, which we have highlighted in the revised manuscript (Results,). However, simplicity does seem a reasonable predictor of order of recruitment into biological use.

The "order of recruitment" correlation with slope (Figure 4) consists of 17 amino acids with little correlation, and these early three that have a large positive slope. This does not convince that age of recruitment continues to affect frequencies. There are many possible amino acid property vectors, and clarity on what has driven these differences is lacking.

On point 7, I am less certain but think we are in agreement. I did not mean to imply the effect should be thrown out, but I think the authors get the point that with three amino acids (GAV) almost wholly responsible for the effect, inferences about which physicochemical properties are responsible for driving the effect are highly speculative, even if reasonable. It sounds like the authors intend to highlight the possible alternatives, so I think this aspect is likely to be fine.

In the analyses shown in Figure 4 we again use Spearman’s, and so the magnitude of the slope will not affect the analysis. It is not surprising that a Spearman’s correlation among 20 points can be destroyed by selectively removing the 3 points that a posteriori most contribute to the correlation, especially when those are the 3 points with the strongest signal in a parametric analysis. This fact is captured within the Figure 4a Spearman’s rho = -0.58, p = 0.008.

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

Article and author information

Author details

  1. Jennifer E James

    Department of Ecology and Evolutionary Biology, University of Arizona, Tucson, United States
    Contribution
    Data curation, Software, Formal analysis, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing
    For correspondence
    jejames@email.arizona.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0518-6783

  2. Sara M Willis

    Department of Ecology and Evolutionary Biology, University of Arizona, Tucson, United States
    Contribution
    Data curation, Software, Formal analysis, Investigation, Visualization, Methodology, Writing - original draft
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1605-6426

  3. Paul G Nelson

    Department of Ecology and Evolutionary Biology, University of Arizona, Tucson, United States
    Contribution
    Data curation, Formal analysis, Methodology, Writing - original draft
    Competing interests
    No competing interests declared

  4. Catherine Weibel

    1. Department of Physics, University of Arizona, Tucson, United States
    2. Department of Mathematics, University of Arizona, Tucson, United States
    Contribution
    Data curation, Software, Investigation, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1837-5209

  5. Luke J Kosinski

    Department of Molecular and Cellular Biology, University of Arizona, Tucson, United States
    Contribution
    Resources, Formal analysis, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8146-5955

  6. Joanna Masel

    Department of Ecology and Evolutionary Biology, University of Arizona, Tucson, United States
    Contribution
    Conceptualization, Supervision, Funding acquisition, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    masel@email.arizona.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7398-2127

Funding

John Templeton Foundation (60814)

  • Joanna Masel

National Institutes of Health (GM-104040)

  • Joanna Masel

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

Acknowledgements

This work was supported by the John Templeton Foundation (60814) and the National Institutes of Health (GM-104040). The authors would also like to thank Heather Meyer for helpful discussions.

Senior Editor

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

Reviewing Editor

  1. Christian R Landry, Université Laval, Canada

Reviewer

  1. Taraneh Zarin, University of Toronto, Canada

Publication history

  1. Received: March 28, 2020
  2. Accepted: January 5, 2021
  3. Accepted Manuscript published: January 8, 2021 (version 1)
  4. Version of Record published: January 21, 2021 (version 2)

Copyright

© 2021, James 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

  • 1,533
    Page views
  • 171
    Downloads
  • 9
    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)

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

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

  1. Jennifer E James
  2. Sara M Willis
  3. Paul G Nelson
  4. Catherine Weibel
  5. Luke J Kosinski
  6. Joanna Masel
(2021)
Universal and taxon-specific trends in protein sequences as a function of age
eLife 10:e57347.
https://doi.org/10.7554/eLife.57347

Categories and tags

Further reading

    1. Evolutionary Biology

    Dating the origin and spread of specialization on human hosts in Aedes aegypti mosquitoes

    Noah H Rose, Athanase Badolo ... Carolyn S McBride
    Research Article Updated

    The globally invasive mosquito subspecies Aedes aegypti aegypti is an effective vector of human arboviruses, in part because it specializes in biting humans and breeding in human habitats. Recent work suggests that specialization first arose as an adaptation to long, hot dry seasons in the West African Sahel, where Ae. aegypti relies on human-stored water for breeding. Here, we use whole-genome cross-coalescent analysis to date the emergence of human-specialist populationsand thus further probe the climate hypothesis. Importantly, we take advantage of the known migration of specialists out of Africa during the Atlantic Slave Trade to calibrate the coalescent clock and thus obtain a more precise estimate of the older evolutionary event than would otherwise be possible. We find that human-specialist mosquitoes diverged rapidly from ecological generalists approximately 5000 years ago, at the end of the African Humid Period—a time when the Sahara dried and water stored by humans became a uniquely stable, aquatic niche in the Sahel. We also use population genomic analyses to date a previously observed influx of human-specialist alleles into major West African cities. The characteristic length of tracts of human-specialist ancestry present on a generalist genetic background in Kumasi and Ouagadougou suggests the change in behavior occurred during rapid urbanization over the last 20–40 years. Taken together, we show that the timing and ecological context of two previously observed shifts towards human biting in Ae. aegypti differ; climate was likely the original driver, but urbanization has become increasingly important in recent decades.

    1. Genetics and Genomics
    2. Evolutionary Biology

    Regulatory and coding sequences of TRNP1 co-evolve with brain size and cortical folding in mammals

    Zane Kliesmete, Lucas Esteban Wange ... Wolfgang Enard
    Research Article

    Brain size and cortical folding have increased and decreased recurrently during mammalian evolution. Identifying genetic elements whose sequence or functional properties co-evolve with these traits can provide unique information on evolutionary and developmental mechanisms. A good candidate for such a comparative approach is TRNP1, as it controls proliferation of neural progenitors in mice and ferrets. Here, we investigate the contribution of both regulatory and coding sequences of TRNP1 to brain size and cortical folding in over 30 mammals. We find that the rate of TRNP1 protein evolution (ω) significantly correlates with brain size, slightly less with cortical folding and much less with body size. This brain correlation is stronger than for >95% of random control proteins. This co-evolution is likely affecting TRNP1 activity, as we find that TRNP1 from species with larger brains and more cortical folding induce higher proliferation rates in neural stem cells. Furthermore, we compare the activity of putative cis-regulatory elements (CREs) of TRNP1 in a massively parallel reporter assay and identify one CRE that likely co-evolves with cortical folding in Old World monkeys and apes. Our analyses indicate that coding and regulatory changes that increased TRNP1 activity were positively selected either as a cause or a consequence of increases in brain size and cortical folding. They also provide an example how phylogenetic approaches can inform biological mechanisms, especially when combined with molecular phenotypes across several species.

    1. Developmental Biology
    2. Evolutionary Biology

    A novel gene REPTOR2 activates the autophagic degradation of wing disc in pea aphid

    Erliang Yuan, Huijuan Guo ... Yucheng Sun
    Research Article

    Wing dimorphism in insects is an evolutionarily adaptive trait to maximize insect fitness under various environments, by which the population could be balanced between dispersing and reproduction. Most studies concern the regulatory mechanisms underlying the stimulation of wing morph in aphids, but relatively little research addresses the molecular basis of wing loss. Here, we found that, while developing normally in winged-destined pea aphids, the wing disc in wingless-destined aphids degenerated 30-hr postbirth and that this degeneration was due to autophagy rather than apoptosis. Activation of autophagy in first instar nymphs reduced the proportion of winged aphids, and suppression of autophagy increased the proportion. REPTOR2, associated with TOR signaling pathway, was identified by RNA-seq as a differentially expressed gene between the two morphs with higher expression in the thorax of wingless-destined aphids. Further genetic analysis indicated that REPTOR2 could be a novel gene derived from a gene duplication event that occurred exclusively in pea aphids on autosome A1 but translocated to the sex chromosome. Knockdown of REPTOR2 reduced autophagy in the wing disc and increased the proportion of winged aphids. In agreement with REPTOR’s canonical negative regulatory role of TOR on autophagy, winged-destined aphids had higher TOR expression in the wing disc. Suppression of TOR activated autophagy of the wing disc and decreased the proportion of winged aphids, and vice versa. Co-suppression of TOR and REPTOR2 showed that dsREPTOR2 could mask the positive effect of dsTOR on autophagy, suggesting that REPTOR2 acted as a key regulator downstream of TOR in the signaling pathway. These results revealed that the TOR signaling pathway suppressed autophagic degradation of the wing disc in pea aphids by negatively regulating the expression of REPTOR2.