A gene-level test for directional selection on gene expression

Abstract

Most variants identified in human genome-wide association studies and scans for selection are noncoding. Interpretation of their effects and the way in which they contribute to phenotypic variation and adaptation in human populations is therefore limited by our understanding of gene regulation and the difficulty of confidently linking noncoding variants to genes. To overcome this, we developed a gene-wise test for population-specific selection based on combinations of regulatory variants. Specifically, we use the QX statistic to test for polygenic selection on cis-regulatory variants based on whether the variance across populations in the predicted expression of a particular gene is higher than expected under neutrality. We then applied this approach to human data, testing for selection on 17,388 protein-coding genes in 26 populations from the Thousand Genomes Project. We identified 45 genes with significant evidence (FDR<0.1) for selection, including FADS1, KHK, SULT1A2, ITGAM, and several genes in the HLA region. We further confirm that these signals correspond to plausible population-level differences in predicted expression. While the small number of significant genes (0.2%) is consistent with most cis-regulatory variation evolving under genetic drift or stabilizing selection, it remains possible that there are effects not captured in this study. Our gene-level QX score is independent of standard genomic tests for selection, and may therefore be useful in combination with traditional selection scans to specifically identify selection on regulatory variation. Overall, our results demonstrate the utility of combining population-level genomic data with functional data to understand the evolution of gene expression.

Keywords: evolution; gene regulation; human evolution; quantitative genetics; selection.

Fig. 1.

We adapted the Qx statistic to test for selection on regulatory variants. a) Spearman $\rho$ between observed and predicted expression in 1 kG for 7,251 JTI models trained in GTEx LCLs, b) and that $\rho$ plotted vs. the in-sample training $R^2$. c) Schematic of Qx calculation as applied to JTI models. The $Q_X$ score is based on the JTI effect sizes and allele frequencies across populations for the set of regulatory variants for each gene. The F matrix contains frequencies across the same populations for variants that were the same frequency in the JTI study population, but were not associated with expression of that gene, thereby modeling the expected covariance for these variants. $Q_X$ is higher when the regulatory variants for a gene exceed those expected patterns.

Fig. 2.

Using a gamma-corrected P-value, we identify 45 genes with evidence of selection. a) Power curves for each P-value method, based on simulations (Methods). We calculated the power for the gamma-corrected version of the $Q_X$ test, as well as for 2 variations of the effect-permuted test. In the first, we drew the effect sizes from the simulation that modeled the corresponding selection strength, and for the other from the neutral effect model. FOS, fitness optimum shift. b) The $Q_X$ score can be decomposed into its $F_{ST}$-like component and its LD-like component. Significant ($\mathrm{FDR}<0.1$) genes in 1 kG for the gamma-controlled and effect-permuted P-values are highlighted in blue and red, respectively, while the red line indicates where $F_{ST}=LD$. The $Q_X$ score for each gene is obtained by adding the two components together. The Spearman rank correlation between the components is 0.15. c) Manhattan plot and d) QQ-plot for gamma-corrected $Q_X$  P-values for 17,388 genes. The horizontal line in (c) at −log₁₀(p) = 3.95 corresponds to $\mathrm{FDR}<0.1$.

Fig. 3.

The $Q_X$ statistic is not correlated with other selection statistics. Pairwise heatmap of Spearman rank correlations between $Q_X$  P-values and various selection-related scores.

Fig. 4.

Different combinations of variant effects can drive predicted differences. a) Median predicted expression in each 1 kG population for the top gene in each peak of the gamma-corrected P-values. For display purposes, for each gene is standardized across populations. b) for FADS1 there is one primary haplotype (tagged by rs174549), while for c) ACO2 there are 3 variants driving the upregulation in PEL. Cells are colored by the product of JTI effect size times effect allele frequency in each population. Values in b) and c) are mean-centered for each variant.