The regulatory genome constrains protein sequence evolution: implications for the search for disease-associated genes

Data

We utilized evolutionary conservation scores from BioMart (Smedley et al., 2015) for human-chimp and human-mouse comparisons. In our analyses, we collected synonymous (dS) and nonsynonymous (dN) substitution rates for 23,816 genes from these comparisons. Proteins vary by two to three orders of magnitude in their rate of evolution (dN/dS).

We utilized GTEx v6p release expression data for 53 tissues involving 8,555 RNA-Seq samples. The mean RPKM and median RPKM were calculated for each gene in each tissue. The variance and ratio of variance to mean were also calculated for each gene in each tissue to estimate inter-individual transcriptional variability (see sample size for each tissue in Table S1). For each gene, we identified the tissue in which the gene shows maximum expression variance. A list of the 53 tissues with information on sample size is in Table S1. In addition to the multi-region sampling of the brain, “frontal cortex” and “cerebellar hemisphere” (obtained after receipt by the brain bank) were sampled by the GTEx Consortium in duplicate (“cortex” and “cerebellum,” respectively, as the first replicate obtained at initial tissue collection). Availability of lymphoblastoid cell lines (LCLs, derived from blood) and cultured primary fibroblasts (from skin) provided an opportunity to compare the expression variance and network constraints (see below) in the cell types with those in the tissues of origin.

For all correlation analyses presented, we used log2-transformed RPKM values. We calculated the Spearman’s ρ and the corresponding p-value with evolutionary rate for all predictors using the R statistical software. To fit a nonlinear model (e.g., of the effect of expression on evolutionary rate), we assumed a generalized additive model using cubic regression spline as “smooth function.”

To test for the robustness of our results to technical and experimental confounding, we calculated the Probabilistic Estimation of Expression Residuals (PEER) (Stegle et al., 2012) factors (which are based on Factor Analysis) and used the residual to test the correlation of expression level and expression variability with evolutionary rate. We compared results across all tissues (Spearman correlation) between pre-PEER and post-PEER analysis.

To quantify tissue expression breadth, we applied the tau (τ) statistic to the 53 GTEx tissues, with the exception that the frontal cortex was the only brain region used so that the highly correlated tissue samples from the brain would not bias the estimate (Yanai et al., 2005; Kryuchkova-Mostacci & Robinson-Rechavi, 2016). ${\displaystyle \tau =\sum _{i=1}^n\frac{1-\widehat{x_i}}{n-1}}$ where ${\displaystyle \widehat{x_i}=x_i∕\underset{1\le j\le n}{max}\left\{x_j\right\}}$ and x_i provides the expression value for the gene in tissue i, and n is the number of tissues. We note that τ yields a score between zero and one, with zero indicating the same expression across tissues and a value of one indicating strong tissue-specific expression. Genes with similarly high levels or similarly low levels across tissues would have high “expression breadth,” which here thus refers to the extent of similarity across tissues. In the actual data, τ attains a minimum of 0.204 and median of 0.82 for protein-coding genes, underscoring the high tissue specificity of a large set of protein-coding genes. Thus, τ also reflects inter-tissue variability (to be distinguished from inter-individual variability in each tissue) in the expression of a gene.

We utilized protein network data in Homo sapiens obtained from STRING v10 (Szklarczyk et al., 2015) to investigate the correlation of the number of interactions with evolutionary rate and expression.

We examined several subsets of genes that define a spectrum of phenotypic effects. We obtained and curated a compendium of Mendelian disease genes from the Online Mendelian Inheritance in Man (OMIM). Essential genes were obtained from a study (Georgi, Voight & Bućan, 2013) that identified human orthologs of mouse genes with known lethal phenotype from the Mouse Genome Database (Blake et al., 2011). LOF-tolerant genes, which can be inactivated without obvious clinical effect, were collected from a comprehensive survey of LOF variants in protein-coding genes (MacArthur et al., 2012). We also investigated immune response genes and olfactory genes, as annotated by the Gene Ontology Consortium (Gene Ontology Consortium, 2015).

Tissue-anchored model of evolutionary rate

We identified the tissue (MaxTissue) in which a gene attains its maximum inter-individual expression variance (MaxVariance)—and, thus, perhaps the tissue in which the gene exhibits its full range of functional activity—and tested the MaxTissue’s association with evolutionary rate. We should emphasize that the MaxTissue is not necessarily the only tissue in which a gene functions, and indeed a gene may have multiple functions in different tissues. A rejection of the null hypothesis would indicate support for the hypothesis (i.e., the “tissue-anchored model”) that when genes are classified into their MaxTissues, the variance ${\sigma }_{across}^2$ in evolutionary rate across tissues is greater than the variance ${\sigma }_{within}^2$ of the tissues.

The MaxTissue’s association with evolutionary rate was tested using the non-parametric Kruskal–Wallis test. To determine whether rejection of the null hypothesis is driven by a single outlier tissue or a class of tissues, we calculated the median evolutionary rate and a metric, MaxTissue-δ, for each MaxTissue, defined as ${\displaystyle 3.14\left(v_{75}-v_{25}\right)∕\sqrt{n}}$ where v₇₅ and v₂₅ are the 75th and 25th percentile for evolutionary rate respectively and n is the number of genes. We also performed empirical analysis by permutation-based tissue assignment (N = 1000), preserving the gene count in the observed MaxTissue configuration, for the entire protein-coding set with available human-mouse data on rates of nonsynonymous and synonymous substitutions. For each permutation of the entire set (RandomTissue), the H-statistic from the Kruskal–Wallis test was calculated. The empirical p-value was the proportion of the total number of permutations in which the RandomTissue H-statistic matched or exceeded the observed value from the MaxTissue assignment.

We can define a novel evolutionary signature, the “evolutionary rate of a tissue,” as the median evolutionary rate of the genes for which the tissue is a MaxTissue. MaxTissue- δ is thus a kind of “confidence interval” for this evolutionary signature. Under the tissue-anchored model, there is a statistically significant difference between tissues for this evolutionary signature. (We note that this signature is to be distinguished from the median evolutionary rate of all expressed genes in a tissue).

The tissue-anchored model starts from the observation that the correlation between within-tissue features (e.g., expression level) and evolutionary rate significantly varies by tissue, raising the centrality of a cross-tissue analysis. The model would suggest that variation between developmental programs in which transcription and other within-tissue features exert their effect may be contributing to variation in evolutionary rate. In each MaxTissue, we ranked the genes according to expression variance from high to low. We then identified the significant Gene Ontology biological processes (Benjamini–Hochberg FDR < 0.05) for the top genes (N = 100) in a MaxTissue. The presence of significant, non-overlapping biological processes between MaxTissues with significant difference in median evolutionary rate (Kruskal–Wallis test) would suggest the importance of developmental processes for constraining evolutionary rate.

Tissue-anchored model: constraint of developmental programs on evolutionary rate

To further investigate the physiological and developmental mechanisms underlying the tissue-anchored model of evolutionary rate, we sought to determine to what extent the MaxTissue for a gene could predict the primary affected tissue of developmental disorders associated with the gene. We utilized a unique UK-wide collaborative resource (generated by over 180 clinicians across 24 regional genetics services), the Deciphering Developmental Disorders (DDD) study (Martin et al., 2018) which has conducted genome-wide genotyping and whole-exome sequencing of children with developmental disorders and their parents. This resource provides a list of genes associated with developmental disorders as well as affected organs. For each gene, we considered the significance of the overlap of MaxTissue with the primary affected organ (using a chi-square test of the null hypothesis of the independence of the MaxTissue and affected tissue).

Testing for independent effects given noisy omics data and estimation of spurious correlation

Partial rank correlation analysis of omics data between two variables, D (a determinant) and K (evolutionary rate), while controlling for a third variable X may generate spurious results (Liu, 1988; Drummond, Raval & Wilke, 2006). Consider the case in which a noisy version X′ is distributed with mean equal to X and variance equal to ${\sigma }_{X^{{}^{\prime }}}^2$: ${\displaystyle X^{{}^{\prime }}\sim \left(X,{\sigma }_{X^{{}^{\prime }}}^2\right).}$

We assume that D and K have the following distributions: $D\sim \left(X,{\sigma }_D^2\right)$ and $K\sim \left(X,{\sigma }_K^2\right)$. Then the partial rank correlation r_DK|X′ between D and K given X′ simplifies to the following expression: ${\displaystyle r_{DK|X^{{}^{\prime }}}=\frac{{\sigma }_{X^{{}^{\prime }}}^2}{f\left({\sigma }_D^2,{\sigma }_K^2,{\sigma }_{X^{{}^{\prime }}}^2\right)}.}$

Here the denominator is some function of all the variances. Importantly, the numerator is nonzero in the presence of noise in X even under the null hypothesis (i.e.,  r_DK|X = 0), leading to a spurious correlation. To determine whether the determinant D contributes to K after controlling for X (given that conventional partial rank correlation analysis generates spurious correlations when applied to noisy biological data), we performed permutation analysis, in which D and X′were shuffled together (thus preserving their correlation) n times (here n = 1,000) and the partial rank correlation between D and K was assessed within each such permutation null set. This generates an empirical p-value for the significance of the observed nonzero correlation coefficient (as the proportion of the permutation null sets that match or exceed the observed correlation coefficient) as well as quantifies the magnitude of the spurious correlation. The permutation null distribution M₀ for r_DK|X′ can be used to estimate an adjusted partial rank correlation coefficient: ${\displaystyle \widehat{{\vartheta }_{DK|X}}=\sqrt{\widehat{r_{DK|X^{{}^{\prime }}}^2}-{\left(\widehat{E\left(M_0\right)}\right)}^2}}$ where $\widehat{r_{DK|X^{{}^{\prime }}}}$ is the observed partial rank correlation coefficient and E(.) is the expectation operator. We call this approach Empirical Partial Rank Correlation Analysis (EPRCA), which facilitates a test for independent effects on evolutionary rate. For example, since expression level and tissue breadth were found to be correlated, we tested their independent effect on evolutionary rate using this approach. The approach, in addition, provides an estimate of the extent of the spurious correlation (from the mean and standard deviation of the permutation null distribution).

Correlation between expression and evolutionary rate given co-expression

The correlation in gene expression (for genes that belong to the same co-expression network) may bias our estimate of the correlation between gene expression and evolutionary rate. We therefore fit a model using Generalized Least Squares (GLS) to account for the non-independence of genes in a co-expression network (shown here, for convenience, for the nonsynonymous substitution rate, but our approach extends, without loss of generality, to dN/dS): ${\displaystyle G=\pi dN+\epsilon }$ ${\displaystyle var\left(\epsilon \right)={\sigma }^2\Delta }$ where dN is the vector of nonsynonymous substitution rates for a set of genes, G is a vector of gene expression, π is the effect size, Δ is the (known) gene expression covariance matrix, and σ² is the unknown (absolute) scale. The GLS effect size estimate solves the following minimization problem: ${\displaystyle \hat{\pi }=\underset{\pi }{\text{argmin}}{\left(G-\pi dN\right)}^T{\Delta }^{-1}\left(G-\pi dN\right)}$ which implies: ${\displaystyle \hat{\pi }={\left(d{N}^T{\Delta }^{-1}dN\right)}^{-1}d{N}^T{\Delta }^{-1}G}$ ${\displaystyle var\left(\hat{\pi }\right)={\left(d{N}^T{\Delta }^{-1}dN\right)}^{-1}{\sigma }^2.}$

We can view the gene expression traits as mapping to a phylogenetic tree such that the covariance matrix Δ captures the covariance between each pair of tips in the tree. The GLS model implements regression that accounts for the phylogeny.

Null phylogenies and null networks

To assess the significance of the correlation between amino acid substitution rate and expression while adjusting for branch assignment (defined in Supplementary Information), we shuffled the evolutionary rate estimates for genes while preserving gene branch assignment within the vertebrate phylogeny. Furthermore, we implemented a node degree-preserving permutation method to generate null gene networks, which were used to control for potential bias that may affect the correlation between evolutionary rate and expression. For branch assignment or node degree, we calculated the empirical p-value as the proportion of phylogeny- or degree- preserving permutations (out of 10,000), respectively, for which the permutation correlation test statistic matched or exceeded the observed statistic in the actual data.

Confounding due to tissue diversity sampling

Despite the comprehensiveness of the GTEx resource, the collection of tissues examined here is still only a partial subset of all human tissues and each tissue is composed of many cell types. We therefore evaluated the robustness to tissue sampling of the correlation between τ and evolutionary rate by calculating the correlation on each of 100 random subsets of 10 and 30 tissues (selected from the 44 tissues).

Cis heritability of gene expression

We estimated the proportion of gene expression variance captured by local genetic variation and quantified its contribution to evolutionary rate. We considered the following linear mixed model: ${\displaystyle Y=W+T+\epsilon }$ ${\displaystyle var\left(Y\right)=A{\sigma }^2+B{\sigma }_T^2+I{\sigma }_{\epsilon }^2.}$ where Y is the residual gene expression phenotype n-dimensional vector after adjusting for hidden factors (with n equal to the number of samples in the reference transcriptome dataset), W is the polygenic cis contribution to gene expression for SNPs within 1 Mb of the gene, A is the genetic relatedness matrix estimated from the local polymorphism data (Yang et al., 2010), T is the polygenic trans contribution from the remaining common variants (MAF > 0.05) in the genome, and B is the genetic relatedness matrix estimated from these trans-variants on the other chromosomes. The variance of the polygenic cis burden W is Aσ² while that of the polygenic trans burden is $B{\sigma }_T^2$; the remaining variance attributable to environmental regulation is $I{\sigma }_{\epsilon }^2$. We estimated these variances using restricted maximum likelihood, as implemented in GCTA (Yang et al., 2011), allowing us to quantify the SNP-based cis heritability of gene expression as $h_Y^2={\sigma }^2∕\left({\sigma }^2+{\sigma }_T^2+{\sigma }_{\epsilon }^2\right)$. We used DGN whole blood samples (N = 922) for maximal power.

Gradient boosted regression for modeling of evolutionary rate

Selection of informative features among many predictors and incorporation of possibly nonlinear effects into a functional form for evolutionary rate may provide insights into potential causal factors and their relative contributions. The task is a variable selection and model choice problem. Although prediction is not our primary aim in this context, a modeling approach that is robust to overfitting and to redundancy (or multicollinearity) is desired. We therefore modeled dN/dS, based on the human-mouse comparison, using gradient boosted regression trees. The approach combines, in an iterative fashion, otherwise “weak” models or classifiers into a “strong” learner. Using a loss function $L\left(y,F\left(x\right)\right)$ on the evolutionary rate y and the model $F\left(x\right)$ built on a vector x of features, the approach seeks to incrementally boost the prediction: ${\displaystyle M_0\left(x\right)=argmi{n}_{\alpha }\sum _{i=1}^{n}L\left(y_i,\alpha \right)}$ ${\displaystyle M_p\left(x\right)=M_{p-1}\left(x\right)+argmi{n}_h\sum _{i=1}^{n}L\left(y_i,M_{p-1}\left(x_i\right)+h\left(x_i\right)\right)}$

Here x_i is a vector of features for y_i, n is the number of training set observations, α is a constant, p is the iteration index, M_p is the p-th model, and h is a base learner (tree) fitted to improve on the model M_p−1. (The residual $y_i-F\left(x_i\right)$ can be interpreted as a negative gradient.)

The models were fitted using the “gbm” R package with a squared error loss function. We utilized a range of values for model parameters. We fitted N_t (100, 500, 1000, 10,000) trees with interaction depth of D_i(1, 4, 5), indicating a model with up to D_i-way interactions (which limits the number of nodes in a tree), and a shrinkage parameter of 0.01 (which shrinks the contribution of each tree). We also performed 10-fold cross validation to evaluate the generalizability of the model. Each combination of model parameters provides a “feature importance score”, generating a predictor importance spectrum that shows the relative importance of each input feature in predicting evolutionary rate. The feature importance score provides support for how MaxTissue is defined in the proposed tissue-anchored model. In addition, using 100 and 500 iterations, N_t = 1000, D_i = 4, and shrinkage parameter of 0.01, we trained on 70% of the data to predict the complement set and to obtain a distribution of out-of-sample adjusted R².

Gene expression imputation using PrediXcan

PrediXcan is an approach for estimating the genetically determined component of gene expression (Gamazon et al., 2015a) using only the germline genetic profile; this estimate is then tested for association with disease risk. Since the disease trait is not likely to modify the germline genetic profile, an observed association comes with a proposed causal direction. Thus, even when directly measured expression is available, its genetic component, though likely to be correlated with the total expression level, provides additional mechanism-relevant information not influenced by disease. Gene expression imputation models were generated using the GTEx reference transcriptome panel in 44 tissues (each with at least 70 samples), as previously described (Gamazon et al., 2015a). Each gene expression model (consisting of selected SNPs within 1 MB of the given gene and their additive regulatory effect size on expression) comes with a measure of imputation performance, namely, the 10-fold cross-validation R² (the square of the correlation between imputed and observed expression) estimated within each tissue (from GTEx or DGN). We note that the best eQTL for a gene also provides an example of an imputation model, though possibly suboptimal in predictive performance. Thus, the practice of testing such an eQTL for its association with a trait is a specific instance, if suboptimal, of the PrediXcan framework. In contrast, our study utilized multi-SNP imputation models for improved performance. The imputation R² also provides an estimate of the aggregate effect of local genetic variation on gene expression and can be evaluated, in each tissue, for its correlation with evolutionary rate.

Methodological implications

We assessed several methodological implications on the search for disease-associated genes using both PrediXcan and conventional GWAS. We explored three scenarios.

1. Extent of local genetic control or cis heritability of gene expression: Cis heritability provides an upper bound on how well local genetic variation may be utilized to impute gene expression. We calculated the Spearman correlation between the imputation R² (derived from 10-fold cross validation for each tissue within the PrediXcan framework (Gamazon et al., 2015a) and evolutionary rate to evaluate to what extent the genetically determined component of the expression of conserved genes, relative to the complement set, may be reliably imputed. We estimated the reduction in statistical power to detect disease associations for conserved genes (dN/dS < 0.01) relative to fast-evolving genes (dN/dS > 1) in transcriptome-wide association studies, assuming a modestly sized GWAS (N = 1000) and using the empirical distribution of R² for the gene sets in a typical tissue (skeletal muscle, chosen for its sample size). The significance was assessed using the Mann–Whitney U test.

2. Relevance of cross-tissue imputation models: We assessed the utility of cross-tissue (versus single-tissue) predictors by evaluating the impact of expression tissue specificity (τ) on gene expression imputation. Genes with τ close to zero have equal expression levels across all tissues, and in this case, tissue-dependent imputation models are likely to suffer from degraded imputation performance; on the other hand, tissue-specific expression profiles would require tissue-dependent imputation models. We tested for enrichment of broadly expressed genes (defined at multiple thresholds, τ < 0.30, 0.40, and0.50) among conserved genes using random sampling from the set of protein-coding genes without replacement (N = 1000). For protein-coding genes, τ has a distribution with a minimum of 0.20 and median of 0.82, indicating high tissue specificity.

3. Stratified analysis of GWAS data using eQTL information: Disease-associated loci identified by GWAS have been shown to be enriched for cis eQTLs (Nicolae et al., 2010; GTEx Consortium, 2015). Indeed, incorporation of eQTL information into GWAS analysis has been shown to improve (quite substantially, for some traits and tissues) the false discovery rate. We calculated the Spearman correlation between the absolute magnitude of the largest cis genetic effect (within 1 Mb of target gene) and evolutionary rate. A strong positive correlation would suggest, for instance, that conserved genes may be less likely to be detected as eGenes. This observation would have important implications on the task of attribution of gene mechanisms to GWAS loci. We calculated the same correlation for the set of trait-associated SNPs (p < 5 × 10⁻⁸) found in the NHGRI GWAS catalog that are in linkage disequilibrium (r² > 0.80) with a best eQTL using the eQTL’s effect size.

Functional enrichment analysis

We performed functional enrichment analyses using the DAVID Bioinformatics Database (https://david.ncifcrf.gov/) and assuming the human genome as background. We analyzed genes with the lowest 10% of τ (i.e., high expression breadth) as well as genes that pass a more stringent threshold (i.e., the 100 most widely expressed genes) to find significantly enriched functional annotations and known pathways for tightly regulated genes. At the other end of the distribution, the 100 genes with the most tissue-specific expression profile were evaluated to identify enriched functional annotations. Similarly, we analyzed genes with 1 to 10 interactions to characterize the genes that map to the oldest branch of the phylogeny with unexpectedly low connectivity. We also tested the 100 genes at the other end of the node degree distribution to characterize the genes that occupy central network positions.

Expression level and variance as tissue dependent correlates of rate of protein sequence evolution

We utilized RNA-seq data derived from 53 human tissues to estimate the mean and variance of expression patterns for each gene (the number of samples for each tissue can be found in Table S1). Throughout the paper, “expression level” refers to the mean expression level. To assess evolutionary conservation, for each orthologous gene we estimated dN/dS rates by comparing human sequences separately with chimpanzee and mouse sequences. We present results from both human-mouse and human-chimp divergence comparisons (Fig. 1), as the human-chimp comparison may be more relevant to the evolution of genes involved in human diseases (discussed later) while the human-mouse comparison is less sensitive to small numbers of mutations.

The rate of nonsynonymous substitution (dN) based on the human-chimp comparison shows a significantly negative correlation (Spearman’s ρ ranges from −0.32 to −0.17, p < 1.41 × 10⁻¹¹⁰) with mean expression level in all 53 tissues. Whenever we present a range for a test statistic, such as Spearman’s ρ, for all tissues, the p-value presented is always the highest or least significant. The same observation holds (Spearman’s ρ ranges from −0.32 to −0.15, p < 1.35 × 10⁻⁷⁷) for evolutionary rate, dN/dS. Stronger correlations with similar trends across tissues are obtained from the human-mouse comparison (Spearman’s ρ ranges from −0.48 to −0.24, p < 1.37 × 10⁻²³⁵ for dN and expression level, and Spearman’s ρ ranges from −0.42 to −0.19, p < 1.35 × 10⁻¹⁴⁸ for dN/dS and expression level).

The uniform analytic pipeline applied to the breadth of tissues (from a large number of the same individuals) enables us to investigate differences between tissues in the correlation between expression in a tissue and support for either purifying selection or positive selection. For each gene in a given tissue, we calculated the mean (to quantify level) and variance or ratio of variance to mean (to quantify inter-individual variability) of gene expression. Notably, the observed expression correlation with sequence divergence varies considerably across tissues and displays tissue clustering (Fig. 1A; human chimp divergence). The brain regions (as previously noted Duret & Mouchiroud, 2000) and, interestingly, cervical spinal cord show the largest absolute magnitude of correlation of expression with evolutionary rate while testis, whole blood, and liver are at the other end of the distribution in both divergence comparisons. The clustering for the remaining tissues enables comparisons of the contribution of physiological systems, including reproductive, immune, and gastrointestinal systems, to sequence divergence (Table S2). Notably, although neural tissues tend to cluster together, gastrointestinal systems and reproductive tissues do not. After the neural tissues, the alimentary canal –in particular, the muscularis mucosae of the esophagus, gastro-esophageal junction, and sigmoid colon –has among the largest absolute effects of expression in primary cells on nonsynonymous substitution rate in both divergence datasets.

We examined the tissue distribution of the correlation of expression with silent substitution rate. Selection on synonymous sites has been previously noted in other taxa (Rocha & Danchin, 2004; Larracuente et al., 2008; Zhou et al., 2016). A negative correlation is observed for dS in most tissues (Table S2) though at lower magnitude than for dN despite the strong correlation between dN and dS (Spearman’s ρ = 0.415, p < 1 × 10⁻³⁰⁰) based on human-chimp divergence. However, this relationship may merely be due to selective constraints on amino acid changes. Previous studies have indeed attributed the correlation between dN and dS to neighboring effects (Duret & Mouchiroud, 2000), and suggested that neighboring sites may influence a site’s mutation rate (Aggarwala & Voight, 2016). Interestingly, the tissue patterns of the correlation of gene expression with synonymous divergence differ from that with non-synonymous divergence. For example, the highest absolute correlation with dS, based on human-chimp divergence, is in fibroblasts, lymphoblastoid cell lines (LCLs), and skeletal muscle (and not in any of the brain regions); in contrast, the cortex now shows the lowest absolute correlation (Table S2). The tissue dependence of the variation of silent substitution rate with expression becomes more significant when using human-mouse divergence, with aorta, sigmoid colon, and tibial nerve (for example) displaying significantly higher absolute effects than the brain regions (Fig. 1A). Collectively, these results support the notion that silent sites may be to some degree under selective constraints.

We assessed the contribution of inter-individual transcriptional variability within each tissue, using both the variance and the ratio of variance to mean (given the correlation between expression mean and variance), to evolutionary rate. Expression variance explains similar levels of variation in evolutionary rate as expression mean in all tissues (Fig. 1B), as expected from the correlation between expression mean and variance. Highly significant results on the effect of expression variance in each tissue on nonsynonymous substitution rate and evolutionary rate are obtained with human-mouse divergence (Spearman’s ρ ranges from −0.46 to −0.18, p < 2.20 × 10 ⁻¹²⁶ for dN and expression variance, and Spearman’s ρ ranges from −0.42 to −0.15, p < 5.20 × 10⁻⁹⁰ for dN/dS and expression variance). Using the ratio of variance to mean to quantify variability, we continue to observe a significant effect on evolutionary rate in each tissue (Spearman’s ρ ranging from −0.39 to −0.08, p < 3.49 × 10⁻²⁶).

To evaluate the robustness of our conclusions, we tested the extent to which hidden and unmeasured confounders in gene expression measurements, including any population structure, may bias our analyses. We calculated the Probabilistic Estimation of Expression Residuals (PEER) factors (Stegle et al., 2012) and utilized the residual (thereby reducing the impact of technical artefacts, population structure, and other hidden confounders that may be present in the measurement of gene expression) for downstream analysis. We continued to observe a significant correlation between variance (using residual expression after PEER adjustment) and evolutionary rate in each tissue (Fig. 1C), with results between the pre-PEER and the PEER-adjusted analysis being significantly correlated (Spearman’s ρ = 0.997).

The central nervous system (CNS) tissues, including the brain and spinal cord, consistently show the largest absolute effects of within-tissue expression variability on evolutionary rate using either variance or the ratio of variance to mean (Table S2 and Fig. 1A), but there is also significant variation in observed effects across the CNS tissues. Based on the human-mouse comparison, highly variable genes (defined as the top 1000 in expression variance) in the cortex are under significantly greater purifying selection than highly variable genes in the subcortical regions (Wilcoxon rank sum p = 1.4 ×10⁻⁴). On the other hand, highly variable genes in the cervical spinal cord evolve more rapidly than highly variable genes in the brain regions (Wilcoxon rank sum range from p = 5.17 ×10⁻⁹ for anterior cingulate cortex BA24 to p = 0.007 for substantia nigra). Whole blood and LCLs tend to show lower contribution of expression variance to variation in evolutionary rate than solid tissues, suggesting that evolutionary analysis in these relatively accessible tissues may not accurately reflect the transcriptome-wide contribution to molecular evolution. Although cell type heterogeneity and Epstein-Barr virus transformation may globally alter gene expression profiles in whole blood and LCLs respectively, we note that spleen, a key component of the immune system, also displays among the lowest absolute effects of expression variance on evolutionary rate (Spearman’s ρ=−0.217, p = 3.45 ×10⁻¹⁸⁷ compared to Spearman’s ρ=−0.419 for anterior cingulate cortex). Of interest to pharmacogenomics studies, highly variable genes in the liver, a key tissue for drug metabolism and transport (Chhibber et al., 2017), tend to show relatively accelerated evolution (Wilcoxon rank sum p = 2.33 ×10⁻¹⁵ for the comparison of muscle and liver), consistent with the hypothesis that pharmacogenes have evolved as a defense mechanism against the accumulation of harmful xenobiotics.

A gene’s MaxTissue is significantly associated with rate of protein sequence evolution

We find that the tissue (“MaxTissue”) in which a gene attains its maximum inter-individual variance (“MaxVariance”) –and perhaps also the tissue in which the gene shows the full range of its functional activity –is significantly associated with evolutionary rate (human-mouse comparison, Kruskal–Wallis test p = 5.55 ×10⁻²⁸⁴). (We note that tissue-specific (τ > 0.95) essential genes demonstrate significantly higher expression variance than tissue-specific non-essential genes (e.g., Mann–Whitney U p = 1.96 ×10⁻⁶⁰ in whole blood), supporting the potential value of considering the tissues of high or maximal variance.) The significantly greater variation in evolutionary rate across MaxTissues (than within) suggests a tissue-anchored model (see Methods) in which the developmental program or physiological process in which transcription takes place may constrain protein sequence evolution. To illustrate, among 3 selected MaxTissues (cortex, liver, and testis) with significant variance (p = 5.55 ×10⁻²⁸⁴) in median dN/dS of the genes in a MaxTissue, the 100 genes with the highest expression variance from each MaxTissue are significantly enriched (Benjamini–Hochberg FDR < 0.10) for non-overlapping sets of functional annotations (see Methods and Fig. 2A). In cortex, dendrite, cell junction, and neuron projections are among the functional categories implicated (p = 4.6 ×10⁻³); in liver: blood microparticle, complement and coagulation cascades, and metabolism of xenobiotics by cytochrome P450 (p = 3.7 ×10⁻¹⁰); in testis: spermatogenesis, chromosome condensation, and multicellular organism development (p = 5.74 ×10⁻⁴). The non-overlapping functional categories and the significant difference (p = 5.55 ×10⁻²⁸⁴) in median evolutionary rate between MaxTissues suggest that variation in the conservation of certain MaxTissue-specific biological and developmental processes significantly contributes to the variation in the rate of protein evolution. Furthermore, permutation analysis that preserves the gene count in the observed MaxTissue configuration finds no random assignment of gene to tissue (RandomTissue; N = 1000 datasets) of the entire protein-coding set that generates a (Kruskal–Wallis) statistic as extreme (see Methods). As expected, the genes with the highest expression variance in RandomTissues show no enriched functional categories. The association is not driven by a single tissue or class of tissues. Using the MaxTissue- δ metric (see Methods) to estimate the variability in evolutionary rate and to rank the MaxTissues, the cortex (0.075), esophagus (0.061), and vagina (0.061) show among the highest values while fibroblasts (0.013), skeletal muscle (0.014), and whole blood (0.015) among the lowest. The other brain regions are spread throughout (e.g., anterior cingulate cortex (0.016), hippocampus (0.017), and caudate basal ganglia (0.051)). Removing the cortex as MaxTissue, the observed association remains significant (p = 1.04 ×10⁻²⁸⁴).

A gene’s MaxTissue predicts primary affected tissue for developmental disorders associated with the gene

To provide additional insights into how developmental programs within the MaxTissues may constrain protein sequence evolution, as proposed by the tissue-anchored model, we tested the extent to which the MaxTissue for a gene could predict the primary affected tissue for developmental disorders associated with the gene using a large-scale, independently curated resource (see Methods). We found that MaxTissue significantly overlapped (p = 3.5 ×10⁻⁴¹) with the primary affected tissue (Fisher’s exact test odds ratio = 4.62), suggesting a potential mechanism for the tissue-anchored model since developmental disorders are likely to affect fitness and therefore the evolutionary rate of the protein.

In addition to MaxVariance, we tested an alternative cross-tissue feature for its effect on evolutionary rate. We used the tau (τ) statistic to capture the level of similarity of expression across tissues (i.e., expression breadth) in GTEx (Yanai et al., 2005; Kryuchkova-Mostacci & Robinson-Rechavi, 2016) (see Methods). In the distribution of τ by gene type (Fig. S1), protein-coding genes display the widest variation; in contrast, lincRNAs show substantial tissue specificity in expression compared to protein-coding genes (Fig. S2). Importantly, using tissue downsampling (see Methods), we find that less comprehensive tissue catalogs than presented here can have substantial variability in the estimated tissue specificity of a gene expression trait (Fig. 2B), with significantly higher variance in τ among genes with higher τ (Spearman’s ρ = 0.359, p < 1.0 × 10⁻³⁰⁰).

We see a significant positive correlation between dN and τ (Spearman’s ρ = 0.32, p = 1 ×10⁻³⁰⁰) and dN/dS and τ (Spearman’s ρ = 0.25, p = 7.5 ×10⁻²⁴³) from the human-mouse comparison, consistent with the notion that intensity of selection on nonsynonymous sites is strongly determined by how broadly expressed the gene is (Duret & Mouchiroud, 2000).

We conjectured that the difference in tissue specificity between non-coding and protein-coding genes may indicate functionally relevant differences between regulatory genes and their targets and provide insights into how expression breadth may constrain evolutionary rate. Indeed, functional enrichment analysis of the lowest 10% of τ scores (i.e., genes with the broadest expression profile) implicates functional annotations relating to transcription (Benjamini–Hochberg adjusted p = 6.19 ×10⁻³⁶) and mRNA processing (Benjamini–Hochberg adjusted p = 1.57 ×10⁻³⁴) (Huang, Sherman & Lempicki, 2009). More stringently, the 100 genes with the greatest expression breadth are enriched in Ubl conjugation pathway (n = 13, Benjamini–Hochberg adjusted p = 3.7 ×10⁻³), a highly conserved eukaryotic gene regulatory mechanism that frequently promotes protein-protein interactions (Hochstrasser, 2009). Thus, a possible explanation for the importance of expression breadth in constraining evolutionary rate is the conservation of genes involved in multi-tissue gene regulation. In contrast, the 100 genes with the most tissue-specific expression profile are enriched for secreted proteins (n = 31, Benjamini–Hochberg adjusted p = 8.5 ×10⁻⁹) with gene products, for instance, located in extracellular region or involved in signaling.

Comparing the two cross-tissue features, we find that τ explains only 0.88% (p = 6.0 ×10⁻⁴⁰) of the variability in MaxVariance. Using EPRCA, we find that evolutionary rate and MaxVariance are significantly correlated after controlling for expression breadth (Spearman’s ρ=−0.151, permutation p < 0.001). These results suggest that expression breadth and MaxVariance, as predictors of between-tissue effects, may independently constrain evolutionary rate.

Joint feature analysis of rate of protein sequence evolution (dN/dS)

Besides MaxTissue and MaxVariance (above), we identified additional correlates of dN/dS. Heritability provides a measure of the potential of a trait to respond to selection (Roff, 2000). Gene expression level as a quantitative trait has a heritable component. For many genes, expression level can be predicted, to a degree measured by heritability and in a tissue-dependent manner, based on genetic polymorphism data (Gamazon et al., 2015a). We find a significant relationship between the cis (common variant) heritability of gene expression and dN/dS (Supplementary Information and Fig. 2C). Conserved genes tend to have lower cis heritability than other genes (Mann–Whitney U p = 1.2 ×10⁻¹⁷). Although gene age and PPI node degree have been previously investigated (and although the method for gene age estimation may generate methodological artifacts (Moyers & Zhang, 2017), we report their relationship with the unique expression data analyzed here (controlling for potential confounders) and find these variables to be significant associated with dN/dS (Supplementary Information).

The features we have examined are mutually correlated (Fig. 3A), posing a challenge to the search for causal factors. A comparison of the generalized additive models based on the human-mouse comparison using AIC (which is asymptotically equivalent to leave-one-out cross-validation for ordinary linear regression (Stone, 1977)), is shown in Table S4. Notably, tissue expression breadth is the model with the best univariate fit.

We then implemented an approach that performs variable selection and model choice (given the number of correlated predictors) and incorporates potentially nonlinear effects (given the univariate observation on the presence of such effects). We utilized gradient boosted regression trees to model evolutionary rate (see Methods), based on human-mouse divergence, jointly analyzing all the features. Notably, we find that MaxVariance, among all features, has the highest feature importance score using a range of model parameters for the number of trees and the number of interactions (see Methods) (Fig. 3B), providing further support to the proper centrality of this feature in our tissue-anchored model of evolutionary rate. In particular, MaxVariance has a stronger influence on predicting evolutionary rate than maximum expression level. We find that the combined model explains an average (and median) of 62% of the variability in evolutionary rate (Fig. 4A) based on out-of-sample prediction (see Methods), with MaxVariance consistently having a higher relative importance than maximum expression level (Fig. 4B). The relative importance of MaxVariance in the gradient boosted model is consistent with the single-model finding that a significant association with evolutionary rate characterizes the MaxTissue configuration and is not observed for a random tissue reconfiguration (permutation p < 0.001; see Methods).

We also considered specific gene sets to determine to what extent the observed relationships between evolutionary rate and expression features apply (Supplementary Information). Mendelian disease genes and essential genes display substantially lower MaxVariance effect, in terms of variability explained, on evolutionary rate than randomly drawn sets (N = 1,000) of proteins, indicating that MaxVariance and gene function may interact to influence evolutionary rate (Supplementary Information and Fig. 4C). As in the full set of protein-coding genes, when restricted to Mendelian disease genes for all model assumptions examined, gradient boosted modeling analyses indicate that MaxVariance continues to have a higher relative importance than maximum expression level.

Methodological implications on genomic studies of disease

Our collective findings have important methodological implications for the search for disease-associated genes and the study of gene function. We explore three specific applications (using the human-chimp comparison for illustration).

The genetically determined component of gene expression in a tissue (see Methods) is typically tested for association with disease in transcriptome-wide association studies (e.g., using PrediXcan (Gamazon et al., 2015a)). Our study would suggest that imputation of genetically determined expression for conserved genes using genetic variants within the cis region of a gene (1 Mb) should be less accurate than for fast-evolving genes. Genes with less accurate expression imputation would suffer from lower statistical power in an association analysis. We note that even if directly measured gene expression data are available, estimation and inference on the genetically determined expression should yield additional information on underlying disease mechanisms.

Our hypothesis is supported by the significant positive correlation across genes between the out-of-sample imputation R² for gene expression (see Methods) and dN/dS (Spearman’s ρ = 0.09, p = 8.6 ×10⁻¹⁴ for adipose subcutaneous; Spearman’s ρ = 0.057, p = 2.6 ×10⁻⁶ for skeletal muscle; and Table S6 for the remaining tissues). For most tissues, imputation performance is significantly lower for conserved genes than for fast-evolving genes (Fig. 5A). The differential performance has important implications for mapping disease-associated genes. Indeed, in a tissue (skeletal muscle, chosen for its sample size) with significant differential imputation performance (Bonferroni-adjusted p < 0.05 for number of tissues tested) and a modestly sized genome-wide association study (N = 1000), we estimate that conserved genes would have significantly lower statistical power than fast-evolving genes (median power 0.16 versus 0.41, Mann–Whitney U test p = 2.2 ×10⁻⁶; see Methods). Furthermore, a smaller proportion of conserved genes (0.32 versus 0.41) would have at least 80% power (Fig. 5B). Finally, the significant correlation with out-of-sample imputation R ² confirms the significant positive correlation between evolutionary rate and cis heritability, for which estimation, using linear mixed models, depends on a rather strong assumption of polygenicity of gene expression.

A second application involves the utility of multi-tissue imputation models. We find that genes involved in multi-tissue gene regulation are notable for the relative strength of purifying selection acting on them, contributing to the strong correlation between expression breadth and evolutionary rate. Our study would suggest that imputation of gene expression for conserved genes could be substantially improved using multi-tissue (versus tissue-specific) genetic predictors. Consider the hypothetical case where a gene is broadly expressed and, indeed, has equal expression across all tissues (i.e., τ = 0). Clearly, an imputation model that varies by tissue would be inconsistent with the gene’s true expression profile and, thus, have suboptimal imputation performance. On the other hand, a gene with highly tissue-specific expression profile (e.g., τ = 1) would require tissue-dependent models for optimal imputation.

Consistent with this hypothesis, among the most highly conserved genes (human-chimp dN/dS < 0.01, N = 2,750), a significant enrichment (p < 0.001) for broadly expressed genes can be found (see Methods). In general, 4.41% of the variance in evolutionary rate is explained by expression breadth. These genes that are subject to the strongest purifying selection show a significant 15% decrease in tissue specificity (τ) (Wilcoxon rank sum p = 9.57 ×10⁻⁷⁸) relative to the remaining protein-coding genes.

For a final application, we investigated whether our findings would also yield methodological insights into traditional single-variant analyses of complex traits. Since trait-associated SNPs identified by GWAS have been shown to be enriched for cis eQTLs, we asked where the target genes of these regulatory variants lie in the overall distribution of evolutionary rate. We calculated the Spearman correlation between the best cis eQTL effect (within 1 Mb of target gene) in absolute value and evolutionary rate. We find a significant positive correlation (Spearman’s ρ = 0.204 between absolute cis effect size and dN/dS, p = 1.3 ×10⁻¹³³), suggesting that conserved genes may be less likely to be detected as eGenes. Furthermore, using reproducible trait-associated variants (p < 5 × 10⁻⁸) curated in the NHGRI GWAS catalog that are in linkage disequilibrium (r ² > 0.80) with a best eQTL, this correlation becomes even greater (Spearman’s ρ = 0.36, p = 6 × 10⁻³²), suggesting that conserved genes may be under-represented as target genes of variants identified by genome-wide association studies. This finding raises the possibility that attribution of gene mechanisms to GWAS loci using only cis eQTLs has underestimated the influence of this class of genes.