Disentangling selection on genetically correlated polygenic traits via whole-genome genealogies

Abstract

We present a full-likelihood method to infer polygenic adaptation from DNA sequence variation and GWAS summary statistics to quantify recent transient directional selection acting on a complex trait. Through simulations of polygenic trait architecture evolution and GWASs, we show the method substantially improves power over current methods. We examine the robustness of the method under stratification, uncertainty and bias in marginal effects, uncertainty in the causal SNPs, allelic heterogeneity, negative selection, and low GWAS sample size. The method can quantify selection acting on correlated traits, controlling for pleiotropy even among traits with strong genetic correlation (|r_g|=80%) while retaining high power to attribute selection to the causal trait. When the causal trait is excluded from analysis, selection is attributed to its closest proxy. We discuss limitations of the method, cautioning against strongly causal interpretations of the results, and the possibility of undetectable gene-by-environment (GxE) interactions. We apply the method to 56 human polygenic traits, revealing signals of directional selection on pigmentation, life history, glycated hemoglobin (HbA1c), and other traits. We also conduct joint testing of 137 pairs of genetically correlated traits, revealing widespread correlated response acting on these traits (2.6-fold enrichment, p = 1.5 × 10^-7). Signs of selection on some traits previously reported as adaptive (e.g., educational attainment and hair color) are largely attributable to correlated response (p = 2.9 × 10^-6 and 1.7 × 10^-4, respectively). Lastly, our joint test shows antagonistic selection has increased type 2 diabetes risk and decrease HbA1c (p = 1.5 × 10^-5).

Keywords: GWAS; adaptation; ancestral recombination graph; polygenic selection; stratification; thrifty gene.

Figure 1

PALM power, calibration, and robustness to uncorrected stratification and ascertainment (A) Left: Power/false positive rate (FPR) of PALM and tSDS. Error bars denote 95% Bonferroni-corrected confidence intervals. Right: PALM selection gradient estimates $\left(\hat{\omega }\right)$. Error bars denote 25th–75th percentiles (thick) and 5th–95th percentiles (thin) of estimates; see Table 1 for more details of $\hat{\omega }$ moments and error. Markers and colors in (A) also apply to (B) and (D). (B) FPR of PALM and tSDS applied to neutral simulations with uncorrected population stratification, simulated with 1000 Genomes data. We used baseline values of ${\sigma }_S=0.1,N_{TSI}/N_{GBR}=1\%$, $M={10}^3$,$h^2=50\%$ by using SNPs ascertained at $\text{p}<5\times {10}^{-8}$. Error bars denote 95% Bonferroni-corrected confidence intervals. (C) Comparison of PALM with true versus Relate-inferred trees, causal versus GWAS-ascertained tag SNPs, and true marginal SNP effects (solid) versus GWAS-estimated SNP effects (hatched). Error bars denote 95% Bonferroni-corrected confidence intervals. (D) Varying polygenicity (M) of the polygenic trait. Error bars denote 95% Bonferroni-corrected confidence intervals. Baseline parameters for all simulations except (C) were our constant-size model with $M={10}^3,$with Scz under positive selection and testing Scz for selection. In (A) and (B), we use Relate-inferred trees and estimated SNP effects at the causal SNPs; in (D), we use Relate-inferred trees and estimated effects at tag SNPs. In all panels, we use a 5% nominal FPR (dashed horizontal line) and simulated ${10}^3$replicates. In (D), light/saturated colors signify neutral/selected simulations.

Figure 2

Joint testing for polygenic adaptation controls for pleiotropy (A) We simulated a cluster of four traits (I–IV) modeled after real human heritability and genetic correlation estimates for schizophrenia (I), bipolar disorder (II), major depression (III), and anorexia (IV), with selection to increase trait I in the last 50 generations. (B and C) We ran marginal and joint tests for selection on these four traits. While marginal selection tests were well powered, they were strongly biased by even fairly low genetic correlations. Conducting a joint test fully controls for genetic correlations while retaining high power to detect and isolate selection on trait I. Simulations (1,000 replicates) were done under our constant effective population size model with $\varrho =60\%,M=\mathrm{1,000},$with trait I under positive selection.

Figure 3

Simulations of joint testing power and calibration (A–F) Differing the degree of pleiotropy $\varrho$ (A), the trait truly under selection (B), the polygenicity M of the traits (C), antagonistic selection on two traits with positive genetic correlation (D), pairwise tests for selection (trait I under selection) (E), and pairwise tests for correlated response (trait I under selection) (F). In (A)–(D), red/pink/blue bars indicate estimates of selection for traits under positive selection/neutrality/negative selection. In (E) and (F), heatmap is colored by positive rate (also text in boxes; standard errors in parentheses). Dashed horizontal lines indicate 5% nominal significance level, and black lines indicate 95% Bonferroni-corrected confidence intervals. Baseline parameters for all simulations (1,000 replicates under each scenario) were our constant-size model with $\varrho =60\%,M=\mathrm{1,000},$with trait I under positive selection. In (A) and (B) and (D), joint tests are performed on trait I/trait III and trait I/trait II, respectively. (E) Diagonal elements correspond to marginal test for selection.

Figure 4

Schematic of “intractable” GxE and induced heritability/genetic correlation Consider Planet A in which there is a homogeneous environment where a policy states that only the fairest (i.e., lightest-skinned) teenagers are admitted to college. If skin color is heritable, then so will be years of education (EduYears) (with some minor attenuation due to truncation with respect to skin color). Likewise, since genetic variants that modulate skin color will also modulate the likelihood of being admitted to college, there will be a genetic correlation between the two traits. However, on Planet B, such a policy might not exist, and accordingly, the heritability of EduYears will be attenuated (as well as its genetic correlation with skin color, which will go to 0%). Unless a study extends into Planet B, GxE effects are indistinguishable from genetic effects.

Figure 5

Estimates of the selection gradient on 56 human traits The selection gradient $\left(\hat{\omega }\right)$ was estimated using 1000 Genomes Great British (GBR) individuals and summary statistics from various GWASs (see Table S4 for full results), with standard errors (${\widehat{se}}_{\omega }$) estimated via block bootstrap ($Z=\hat{\omega }/{\widehat{se}}_{\omega }$). Starred traits indicate significance at FDR = 0.05.

Figure 6

Correlated response in real traits (A) Expanded view of antagonistic selection on glycated hemoglobin (HbA1c) versus type 2 diabetes (T2D). We estimate individual SNP selection coefficients by taking the maximum-likelihood estimate $\hat{s}$ for each SNP. We plot this value against the joint SNP effect estimates for HbA1c and T2D. Colored lines represent isocontours of $s\left(\beta \right)={\beta }_{HbA1c}{\hat{\omega }}_{HbA1c}+{\beta }_{T2D}{\hat{\omega }}_{T2D}$, the estimate of the Lande transformation from SNP effects to selection coefficients, where $\hat{\omega }$ is inferred jointly for the two traits (Table 3). The purple-green color gradient illustrates expected selection coefficients under $\hat{\omega }$ (background) versus individual SNP selection coefficient estimates (rings). Ring diameter is proportional to SNP selection log-likelihood ratio. (B) Enrichment of correlated response in analysis of genetically correlated traits. Enrichment in the tails of the distribution of our test statistic for correlated response $R(\text{p}=1.5\times {10}^{-7}$, binomial test), which had 2.6-fold enrichment at the nominal 5% level. We assessed $n=2\times 137=274$estimates of correlated response on 137 trait pairs with Bonferroni-significant ${\text{p}}_{r_g}<0.005/\frac{56}{2}$ and $|r_g>0.20|$. Red area indicates pointwise 95% CI of the survival curve.