An approach to identify gene-environment interactions and reveal new biological insight in complex traits

Abstract

There is a long-standing debate about the magnitude of the contribution of gene-environment interactions to phenotypic variations of complex traits owing to the low statistical power and few reported interactions to date. To address this issue, the Gene-Lifestyle Interactions Working Group within the Cohorts for Heart and Aging Research in Genetic Epidemiology Consortium has been spearheading efforts to investigate G × E in large and diverse samples through meta-analysis. Here, we present a powerful new approach to screen for interactions across the genome, an approach that shares substantial similarity to the Mendelian randomization framework. We identify and confirm 5 loci (6 independent signals) interacted with either cigarette smoking or alcohol consumption for serum lipids, and empirically demonstrate that interaction and mediation are the major contributors to genetic effect size heterogeneity across populations. The estimated lower bound of the interaction and environmentally mediated heritability is significant (P < 0.02) for low-density lipoprotein cholesterol and triglycerides in Cross-Population data. Our study improves the understanding of the genetic architecture and environmental contributions to complex traits.

Introduction

Current genome-wide association studies (GWAS) focus on detecting genetic variants that lead to different phenotype means across genotype groups^1,2, and have identified a large number of genetic loci that, in some cases, explain large proportions of the trait’s SNP-heritability^3,4,5. While it is commonly agreed that complex traits are influenced by genetic and environmental factors and their interactions^6,7,8,9, there is a long-standing disagreement about the magnitude of the $G\times E$ contribution to heritability because of different theoretical models and assumptions^10,11. As pointed out in ref. ¹², arbitrarily defined parameterizations of genetic effects with non-additive gene actions may explain the same degree of genetic variation as the currently prevailing additive model. Thus, while using additive genetic models such as polygenic risk scores to predict individual quantitative or qualitative phenotypes has become standard⁵, these models may not be fully informative in understanding genetic architecture.

Interactions are often studied secondarily in comparison to additive variance, whose advantage is to explain most of the correlations among relatives and fit natural selection model well^10,13. Theoretical studies have demonstrated that a significant portion of variance can be explained by an additive model even when the genetic contribution to a phenotype is purely through $G\times E$¹⁴. This limitation is one of the key factors explaining the low power of approaches modeling interactions conditional on additive variance. As a result, studies focusing on detecting $G\times E$ at the genome-wide level are seldom considered as primary analyses. By contrast, the joint evidence of main genetic and $G\times E$ effects, in addition to the $G\times E$ alone, is tested in the Gene-Lifestyle Interactions (GLI) Working Group within the Cohorts for Heart and Aging Research in Genetic Epidemiology Consortium (CHARGE)^9,15, where only a modest number of genetic loci have been identified through testing for $G\times E$ alone^16,17,18,19. The joint test limits our ability to determine to what degree the currently identified loci reflect evidence for $G\times E$ contribution, making it difficult to understand the precise interplay between genetic and environmental factors.

Concurrently, Mendelian randomization (MR) has been developed and widely applied to study causal relationships between risk exposures and outcomes in the post-GWAS era^20,21. Although an MR approach has been used to explain $G\times E$²², the underlying connection between testing pleiotropic variants in the MR framework and the detection of $G\times E$ is currently unclear. Here, we conceptually connect $G\times E$ with MR framework, illuminate their similarities and demonstrate that the test of horizontal pleiotropy in MR^23,24 can be used for detecting $G\times E$. Based on this principle, one can identify $G\times E$ using existing available GWAS and GWIS summary statistics. We applied this approach to the summary statistics from the Global Lipids Genetics Consortium study (GLGC, n = 1.65 M)³ and the summary statistics in the interaction analysis with cigarette smoking and alcohol drinking in the CHARGE GLI working group¹⁷, with replication using direct interaction tests performed in the UK Biobank (UKBB) data. Although the UKBB data accounted for about one third of sample in the GLGC consortium, theoretical work suggests that such replication is statistically independent (Supplementary Note)

Results

Testing $G\times E$ and mediation based on Mendelian randomization (MR)

Traditionally a genome-wide interaction study (GWIS) with an environmental exposure on a quantitative trait Y is modeled through a linear regression:

$$Y={\beta }_{0}+{\beta }_{1}G+{\beta }_{2}E+{\beta }_{3}G\times E+\epsilon,$$

(1)

where ${\beta }_{1},{\beta }_{2}$ and ${\beta }_{3}$ correspond to the ‘main’ effect of $G$ (in the presence of $E$), the main effect of $E$ and the interaction effect of $G\times E$, respectively, and $\epsilon$ is a random noise. Here $G$, $E$, and $G\times E$ represent a genotype value, environmental factor, and their interaction respectively. For simplicity, we do not include any covariates, but it will not affect the general conclusion. The interaction effect is evaluated by the direct test statistic ${T}_{{direct}}={\hat{\beta }}_{3}^{2}/{{{{\mathrm{var}}}}}\left({\hat{\beta }}_{3}\right)$, where ${\hat{\beta }}_{3}$ refers the estimate from the regression model (1). Theoretical work indicates that the test statistics for the main effect ${\beta }_{1}=0$ and the interaction effect ${\beta }_{3}=0$ are correlated, with the correlation coefficient equal to $-{\mu }_{E}/\scriptstyle\sqrt{{\mu }_{E}^{2}+{\sigma }_{E}^{2}}$, where ${\mu }_{E}$ and ${\sigma }_{E}^{2}$ are the mean and variance of the environmental factor in the data¹⁴. However, the power of the direct test is usually low because of the collinearity between $G$ and $G\times E$ which induces a covariance between the estimates of ${\beta }_{1}$ and ${\beta }_{3}$. This covariance produces uncertainty (i.e. larger standard error) which by itself reduces power for testing either ${\beta }_{1}$ or ${\beta }_{3}$, even if the underlying true model includes $G\times E$ alone (i.e., ${\beta }_{1}=0$ and ${\beta }_{3}\ne 0$)^10,14.

In practice, a GWAS is routinely conducted first when studying the genetic contribution to a trait, which is typically done through a linear regression model without including environmental factors, i.e.,

$$Y={\alpha }_{0}+\alpha G+\epsilon,$$

(2)

where we refer to $\alpha$ as the ‘marginal’ effect from a GWAS (in the absence of $E$) to differentiate from the main effect ${\beta }_{1}$ in model (1). We show that $\alpha -{\beta }_{1}=\frac{\rho {\sigma }_{E1}}{{\sigma }_{G1}}{\beta }_{2}+({\mu }_{E1}+\frac{\rho {\sigma }_{E1}}{{\sigma }_{G1}}){\beta }_{3}$, where $\rho$ is the mediation contribution of $G$ through $E$, ${\mu }_{E1}$, ${\sigma }_{E1}$, and ${\sigma }_{G1}$ represent the environmental mean, standard deviation, and genotype standard deviation in GWAS data, respectively, suggesting that testing the hypothesis H₀: α–β₁ = 0 for the difference between the marginal and main effects is equivalent to testing for the combined effect of $G\times E$ and mediation, and further reduces to testing for the $G\times E$ when $G$ and ${E}$ are independent (i.e., $\rho=0,$ Supplementary Note). This hypothesis can be tested by the statistic $\,{T}_{{diff}}$$={(\hat{\alpha }-{\hat{\beta }}_{1})}^{2}/{{{{\mathrm{var}}}}}(\hat{\alpha }-{\hat{\beta }}_{1})$, where $\hat{\alpha }$, ${\hat{\beta }}_{1}$, and their corresponding standard errors are estimated from the GWAS and GWIS analyses, respectively. In fact, ${T}_{{diff}}$ and ${T}_{{direct}}$ are equivalent when GWAS and GWIS are performed in the same data. We verified this property using real data analysis in the GLI studies¹⁷, from which the summary statistics of the marginal, main, and interaction effects are available and the marginal effect was obtained after adjusting for $E$. We observed that the correlation between the statistics of the ${T}_{{diff}}$ and the direct test is 0.98 for LDL and current smoking (Supplementary Fig. 5). However, GWAS is often performed in a much larger sample than the GWIS because of data availability. The environmental exposure may have different distributions in cohorts for conducting GWAS and GWIS (i.e., different mean and variance). Furthermore, models (1) and (2) are likely to be performed by two different groups of investigators, which will bring variation across studies in trait definitions, trait measurement procedures, quality control procedures, and covariates. Moreover, the summary statistics are obtained through meta-analyses in both GWAS and GWIS analyses, which can bring additional variation and confounding factors, including population stratification and cryptic relatedness, leading to a potentially invalid comparison between the marginal and main effects. In fact, it has been reported that the confounding of population stratification is not sufficiently corrected in large GWAS^25,26. Therefore, directly using ${T}_{{diff}}$ to screen the genome can be biased even for testing the combined contribution of interaction and mediation.

To overcome this bias, we note that the marginal effect estimate $\hat{\alpha }$ and the main effect estimate ${\hat{\beta }}_{1}$ have a linear relationship,

$$\hat{\alpha }=\theta {\hat{\beta }}_{1}+\frac{\rho {\sigma }_{E1}}{{\sigma }_{G1}}{\hat{\beta }}_{2}+\left({\mu }_{E1}+\frac{\rho {\sigma }_{E1}}{{\sigma }_{G1}}\right){\hat{\beta }}_{3},$$

(3)

where $\theta$ reflects the contribution of main effect to marginal effect, which converges to 1 when GWAS and GWIS are conducted using homogeneous measurements of phenotypes and environments (“Methods”). The genetic variants with no $G\times E$ and no mediation will fall on the regression line but the variants with $G\times E$ or mediation will depart from this line. We do not expect this pattern to be systematically impacted by the variation across studies. Therefore, we search the genetic variants that depart from this regression line to test the combined effect of $G\times E$ and mediation, providing $\theta$ can be correctly estimated. This idea is conceptually the same as the MR framework when we introduce a pseudo exposure $\widetilde{X}$, representing a polygenic score of the trait (Fig. 1). We do not need to construct this pseudo exposure in our analysis because we work directly on the summary statistics under the MR framework. We then estimate the causal effect $\theta$ of the pseudo exposure $\widetilde{X}$ on trait $Y$ in the MR framework and the $G\times E$ effect or mediation through $E$ is tested in the same way as testing for horizontally pleiotropic variants²³. In doing so, we first select a set of independent variants associated with trait $Y$ and perform the inverse variance weighted analysis to estimate θ, ${{{{{\rm{denoting}}}}}}$ ${{{{{\rm{as}}}}}}$ $\hat{\theta }$. Second, we test the $G\times E$ or mediation of a genetic variant through $E$ by the statistic ${T}_{{MR\_GxE}}=\frac{{\left(\hat{\alpha }-\hat{\theta }{\hat{\beta }}_{1}\right)}^{2}}{{{{{\mathrm{var}}}}}\left(\hat{\alpha }-\hat{\theta }{\hat{\beta }}_{1}\right)} \sim {\chi }_{1}^{2}$. This test can be performed by the iterative Mendelian randomization and pleiotropy (IMRP) approach^23,27. The statistic ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$ is an asymptotically unbiased test for testing the combined effect of $G\times E$ and mediation through $E$ (Supplementary Note).

Fig. 1: Illumination of Mendelian randomization and G × E.

A Left panel: the path diagram of the MR, where U refers to all confounders. Genetic variants (G) contributing to outcome Y through mediation of exposure X are often selected as the valid genetic instrumental variables (black paths). Genetic variants contributing to Y through both black and red paths independently are horizontal pleiotropic variants. Genetic variants contributing to Y through confounders (U) are invalid instrumental variables and need be blocked (x). Right panel: a scatter plot of effect sizes of genetic instrumental variants for an exposure and an outcome. Each + corresponds to the 95% confidence intervals of the exposure effect size (horizontal line segment) and the outcome effect size (vertical line segment). The horizontal pleiotropic variants (red +) depart from the regression line and can be separated from the variants with no pleiotropic effect (blue +). B Left panel: the $G\times E$ framework, with the goal of testing $G\times E$. Instead of an explicit exposure, we create a pseudo exposure $\widetilde{X}$, which can be viewed as a polygenic score for trait Y based on marginal effect sizes. However, our analysis does not require estimating this pseudo exposure. The genetic variants associated with the pseudo exposure $\widetilde{X}$ but not through either the environment E or $G\times E$ are valid instrumental variables. The genetic variants interacting with E can be viewed the same as horizontally pleiotropic variants in the MR framework. Genetic variants associated with Y via mediation through E can contribute to both the pseudo exposure $\widetilde{X}$ and Y, and thus have similar effects as $G\times E$ and cannot be distinguished from $G\times E$. Thus, testing the combined effect of interaction and mediation is conceptually equivalent with testing the horizontally pleiotropic effect in the MR framework. Right panel: a scatter plot of genetic variants for GWIS main effects and GWAS marginal effects. Each + corresponds to the 95% confidence intervals of the GWIS main effect size (horizontal line segment) and the GWAS marginal effect size (vertical line segment). Like the horizontal pleiotropic variants in the MR framework, $G\times E$ variants (red +) depart from the regression line and can be separated from variants with no $G\times E$ assuming no mediation.

Two-step procedure for testing $G\times E$

Note that ${T}_{{MR\_GxE}}$ likely tests for the combined effect of $G\times E$ and mediation unless $G$ and ${E}$ are independent (i.e., $\rho=0$). To test for $G\times E$, we propose a two-step procedure by using ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$ to screen the whole genome and then performing ${T}_{{direct}}$ on the variants surviving the ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$ screen. We set the significance level at 5 × 10⁻⁸ for the first step (${T}_{{MR\_GxE}}$ test), and the significance level at 0.05/X for the step 2 ${T}_{{direct}}$ test, where X is the number of independent significant variants in the first step/test. This two-step procedure can increase power at the screening step when there is interaction and mediation and increases power at the direct testing step by substantially reducing the multiple comparison burden. ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$ and ${T}_{{direct}}$ are not independent (Supplementary Note), and therefore, the variants detected by the two-step procedure could still reflect the contribution of mediation and $G\times E$, and it is necessary for further replication by performing ${T}_{{direct}}$ in independent data. To mitigate this problem, we can exclude the variants identified through GWAS of $E$, which could represent large mediation effect.

Type I error rate and power of ${T}_{{MR\_GxE}}$ and the two-step procedure

We first performed a series of simulations to investigate the type-I error rate and power of ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$ in the absence of mediation. In simulations we observed that $E$($\hat{\theta }$) is close to 1 and the estimate $\hat{\theta }$ converges to 1 when sample size increases, which is expected by theoretical prediction (Fig. 2A and Supplementary Fig. 6a). The direct estimate of the interaction effect ${\beta }_{3}$ as well as of $\left(\hat{\alpha }-{\hat{\beta }}_{1}\hat{\theta }\right)/{\mu }_{E}$ is also unbiased (Fig. 2B, Supplementary Fig. 6b), although the standard error of $\left(\hat{\alpha }-{\hat{\beta }}_{1}\hat{\theta }\right)/{\mu }_{E}$ is affected by the environmental means in GWAS and GWIS. When no mediation is present, the type-I error rates for both ${T}_{{MR\_GxE}}$ and the direct test are well controlled (Fig. 2C and Supplementary Fig. 6c)). The power of ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$ depends on multiple parameters, including ${\mu }_{E}$ and allele frequency in GWAS and GWIS and is less powerful than ${T}_{{Direct}}$ when the environmental mean in GWAS is lower (Fig. 2D and Supplementary Fig. 6d). Additional simulations for the estimates of $\hat{\theta }$, interaction effect $\left(\hat{\alpha }-{\hat{\beta }}_{1}\hat{\theta }\right)/{\mu }_{E}$, type-I error rate and power are presented in Supplementary Figs. 7–9.

Fig. 2: Simulation performance of T_{MR_GXE} and the two-step procedure.

A–D No medication was present. The simulation details were described in “Methods”. A Box plots of $\hat{\theta }$ in simulations under different environments in GWAs data. The top and bottom edges of the box plots represent the 25th and 75th percentiles of $\hat{\theta }$, and the horizontal middle line represents the 50th percentile. The vertical bars extend from the 25th (or 75th) percentile of $\hat{\theta }$ to the minimum (or maximum) value of simulated data. $E\left(\hat{\theta }\right)$ is close to 1 as expected. B Box plots of the direct estimate of ${\beta }_{3}$ in GWIS (top panel) and by $\left(\hat{\alpha }-{\hat{\beta }}_{1}\hat{\theta }\right)/{\mu }_{e}$ through MR- $G\times E$ analysis (bottom panel). The box plots are interpreted the same as in (A) accordingly. Both the estimates of ${\beta }_{3}$ and that by $\left(\hat{\alpha }-{\hat{\beta }}_{1}\hat{\theta }\right)/{\mu }_{e}$ are unbiased. Here s = −1 refers to the scenario when the main effect and interaction effect have opposite effect directions; s = 0 refers to no main effect; and s = 1 refers to the scenario when the main effect and interaction effect have the same effect direction. C Type I error rate comparison between ${T}_{{MR\_GxE}}$ and the direct test for different main and interaction effect directions. Both ${T}_{{MR\_GxE}}$ and the direct test maintain the type I error rate well. D Power comparison between ${T}_{{MR\_GxE}}$ and the direct test for different main and interaction effect directions. E, F 20 variants were tested when mediation was present or not. The simulation details were described in ”Methods”. E Type I error comparison for ${T}_{{Direct}}$, ${T}_{{MR\_GxE}}$ and two-step procedure. The dash lines represent the 95% confidence interval. F Power comparison for ${T}_{{Direct}}$, ${T}_{{MR\_GxE}}$ and two-step procedure.

We next investigated the performance of ${T}_{{Direct}}$, ${T}_{{MR\_GxE}}$ and the two-step procedure when mediation is present and multiple variants are tested. We generated 20 independent variants with one variant having mediation, interaction, or both. All three tests have well controlled type I error rates when mediation is absent (Fig. 2E and Supplementary Fig. 10A). When mediation is present, the type-I error rate was still well controlled, although inflation can be observed for the two-step test and ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$ when $E$ contributes to 5% of the outcome variation and the samples between GWAS and GWIS are completely overlapped (Supplementary Fig. 10B, C). This inflation was caused by mediation and quickly disappeared when the overlapping rate between GWAS and GWIS subjects was reduced. The statistical power of ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$ and the two-step procedure for testing $G\times E$ was much more improved than ${T}_{{Direct}}$ when mediation was present (Fig. 2F and Supplementary Fig. 10D–F).

Identifying gene-smoking and gene-alcohol drinking interactions to serum lipids

We applied the two-step procedure to search for genetic variants interacting with cigarette smoking and alcohol drinking for serum lipids, using the summary statistics of high-density lipoprotein cholesterol (HDL-C), low-density lipoprotein cholesterol (LDL-C), and triglycerides (TG) from the GLGC (n = 1.65 M) and the CHARGE GLI (n = 134 K). To mitigate the effects of mediation through cigarette smoking or alcohol drinking, we excluded all loci with P-value < 5 × 10⁻⁷ reported in the early GWAS of cigarette smoking status or alcohol drinking²⁸, which represent relatively large effect sizes of variants on cigarette smoking and alcohol drinking. We observed that $\hat{\theta }$ ranged from 0.92−1.33, 0.95−1.62, 0.83−1.25, 0.87−1.37, and 0.95–1.28 for European, African, Asian, Hispanic, and Cross-population data, respectively (Supplementary Data S1). The departure of $\hat{\theta }$ from 1 suggests that the phenotype treatments, analysis protocols, and corrections for population structure were not identical between the GLGC and CHARGE GLI consortiums. For example, CHARGE GLI performed a natural logarithmic transformation to the lipid measurements, whereas GLGC further performed an inverse normal transformation. The number of principal components (PCs) for correcting populations was also different between GLGC and CHARGE GLI. Despite these discrepancies, we did not observe an inflation for ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$, with the genomic control $\lambda$ values being close to 1 (range 0.93–1.05, Supplementary Data S2).

Using ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$ to screen the genome, we observed 15 genome-wide significant loci consisting of 17 independent signals (P < 5 × 10⁻⁸, ${r}^{2} < 0.1$), including 4 and 5 loci for LDL-C, 7 and 5 loci for HDL-C, and 5 and 6 loci for TG, interacting with cigarette smoking and alcohol drinking or mediating through them, respectively (Fig. 3A–C, G–I, Supplementary Data S3a). All but 3 loci have been reported to be associated with either cigarette smoking or alcohol drinking in the recent largest GWAS study with over 3 million samples²⁹, suggesting the contribution of both $G\times E$ and mediation. Since we already excluded the cigarette smoking and alcohol drinking variants identified from a relatively smaller study²⁸, these detected variants should represent modest mediation effects. Locus-specific plots of all significantly associated loci are presented in Supplementary Fig. 11, which suggest that multiple protein-coding genes are present in these loci. Strikingly, all the loci have previously been mapped to lipids traits except RPL5P26 on chromosome 10. The $G\times E$ or mediation loci are clearly departing from most of the lipids-associated variants (Fig. 3D–F, J–L). The population-specific ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$ results are presented in Supplementary Fig. 12, which are also consistent with the Cross-population results, although the main contribution comes from the European population.

Fig. 3: Manhattan plots, marginal and main effect size comparisons.

The circle Manhattan plots of gene × alcohol drinking interactions for A LDL-C; B HDL-C; and C TG, respectively. The genome-wide significant loci are presented in red dots. The marginal and main effect sizes corresponding to alcohol drinking for D LDL-C; E HDL-C, and F TG, respectively. The colored circles represent the genome-wide significant loci and gray circles represent insignificant loci by ${T}_{{MR\_GXE}}$ test. The circle Manhattan plots of gene × cigarette smoking interactions for G LDL-C; H HDL-C; and I TG, respectively. The marginal and main effect sizes corresponding to cigarette smoking for J LDL-C; K HDL-C, and L TG, respectively.

By applying the two-step procedure, we observed that 8 of the 17 independent signals are significant when using the direct test ${T}_{{Direct}}$ after correcting for the 17 tests and 4 environmental factors (Table 1, P < 7.35 × 10⁻⁴). In comparison to the direct test in GWIS, the two-step procedure identified more $G\times E$ signals for each of the three lipid traits and four environmental factors (Supplementary Data S3b). This provides additional support for the enhanced statistical power of the two-step procedure. The tissue enrichment analysis using the GWAS-based pathway analysis tools, MAGMA³⁰ and FUMA³¹, suggest that these loci are enriched in liver, hippocampus, small intestine, and stomach tissues (Supplementary Fig. 13). Multiple loci were colocalized with expression quantitative trait loci (eQTLs) in the corresponding liver, lung, and blood tissues in the genotype-tissue expression database (GTEx)³² (Supplementary Fig. 14).

Table 1 Interaction loci screened by ${T}_{{MR\_GxE}}$ and followed by the direct test ${T}_{{Direct}}$ in GLI (two-step test) and replicated by the direct test ${T}_{{Direct}}$ in UK Biobank

Independent replication

We next attempted to replicate the evidence for these 8 independent signals in the UKBB. Although the UKBB data accounted for about one third of samples in the GLGC consortium, the direct test statistic ${T}_{{Direct}}$ calculated in the UKBB is independent of ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$, so are the ${T}_{{Direct}}$ test statistics calculated in UKBB and CHARGE GLI, thus qualifying this as an independent replication (Supplementary Note). Six of the 8 signals are replicated in the UKBB after adjusting for 32 tests (P < 1.56 × 10⁻³), and 5 of them were genome-wide significant by the ${T}_{{Direct}}$ test in combined CHARGE GLI and UK Biobank data (Table 1). All 8 independent signals have the same interaction direction in CHARGE GLI and UKBB except LPL, which is not significant in UKBB (Supplementary Data S3a). The CETP and SMARCA4 loci are the only two loci with no reported mediation evidence through either cigarette smoking or alcohol drinking.

We now aim to understand if the interaction evidence observed in this study has an alternative explanation³³ because of linkage disequilibrium (LD) with a variant which has causal effect on cigarette smoking or alcohol drinking. To examine this, we searched if there exists a variant(s) at each of the loci in Table 1 explaining the observed interaction evidence in the UKBB. However, we did not observe such variants (Supplementary Fig. 15), suggesting that the interaction evidence presented in Table 1 is genuine. In total, we identified 5 loci consisting of 6 independent signals that have evidence of interaction with either cigarette smoking or alcohol drinking.

$G\times E$ interaction and mediation to SNP heritability

Since $\hat{\alpha }-{\hat{\beta }}_{1}\hat{\theta }$ refers to the combined interaction and mediation contribution to the marginal effect, we use $\hat{\alpha }-{\hat{\beta }}_{1}\hat{\theta }$ to estimate the heritability contributed by the interaction and mediation through the LD score (LDSC) regression³⁴. Note that this heritability is a lower bound of the phenotype variance contributed by the $G\times E$ and mediation through E and is a part of the heritability estimated through the marginal effect, which is often referred to as the SNP-heritability in GWAS. In both Cross-Population (Fig. 4A) and European population (Fig. 4B), we observed significant interaction and mediation heritability (P < 0.03) with ever cigarette smoking for LDL-C, and alcohol consumption or cigarette smoking for TG, suggesting that the heritability estimates based on marginal effects also include significant contributions from $G\times E$ and mediation through the corresponding environment factors (Supplementary Data S4).

Fig. 4: The estimated heritability of HDL-C, LDL-C, and TG using LDSC regression.

A Cross-Population. B European population. X-axis represents heritability in percentage. Y-axis represents the corresponding heritability estimated in percentage (marginal.effect: marginal effect heritability; current.drinking: gene and current drinking interaction effect heritability; regular.drinking: gene and regular drinking interaction effect heritability; current.smoking: gene and current smoking interaction effect heritability; ever.smoking: gene and ever smoking interaction effect heritability). Marginal effect heritability refers to the heritability estimated through the marginal effect $\hat{\alpha }$, and interaction effect heritability refers to the heritability estimated through $\hat{\alpha }-\hat{\theta }{\hat{\beta }}_{1}$. The percentage number displayed on the right side of each bar represents the estimated heritability, and the corresponding 95% confidence interval shown as horizontal error bars. For the marginal effect analysis, the sample size is 1.65 M and 1.32 M for cross-population and European population analysis, respectively. For the interaction effect analysis, the sample size is 134 K and 80 K for cross-population and European population analysis, respectively.

$G\times E$ interaction and mediation to heterogeneity of genetic effect sizes across populations

As noted in Eq. (3), the marginal effect estimate of a genetic variant in GWAS consists of the $G\times E$ and mediation contribution when the $G\times E$ and mediation occur. Because of the environment heterogeneity across populations, we expected that the marginal effect sizes of the variants will be less correlated across populations for the variants with than without $G\times E$ interaction or mediation. We calculated the marginal effect size correlations between European, African, Hispanics, and Eastern Asian for these variants reported in Graham et al³ after excluding the variants in Supplementary Data S3a where their $G\times E$ interactions or mediations were observed in this study. Similarly, we calculated the marginal effect size correlations for the variants in Supplementary Data S3a. We compared the correlation and observed a median of 24.4% drop of the cross-population correlation coefficient (Fig. 5), strongly support that $G\times E$ interactions or mediations contribute to the marginal effect size heterogeneity across populations.

Fig. 5: Cross-population comparison of the LDL-C, HDL-C, and TG marginal effect sizes of the variants reported in Graham et al.³.

A EUR vs AFR, no ${{{{{\boldsymbol{G}}}}}}{{{{\times }}}}{{{{{\boldsymbol{E}}}}}}$ interaction or mediation. B EUR vs AFR, ${{{{{\boldsymbol{G}}}}}}{{{{\times }}}}{{{{{\boldsymbol{E}}}}}}$ interaction or mediation. C EUR vs HIS, no ${{{{{\boldsymbol{G}}}}}}{{{{\times }}}}{{{{{\boldsymbol{E}}}}}}$ interaction or mediation. D EUR vs HIS, ${{{{{\boldsymbol{G}}}}}}{{{{\times }}}}{{{{{\boldsymbol{E}}}}}}$ interaction or mediation. E EUR vs EAS, no ${{{{{\boldsymbol{G}}}}}}{{{{\times }}}}{{{{{\boldsymbol{E}}}}}}$ interaction or mediation. F EUR vs EAS, ${{{{{\boldsymbol{G}}}}}}{{{{\times }}}}{{{{{\boldsymbol{E}}}}}}$ interaction or mediation. G AFR vs HIS, no ${{{{{\boldsymbol{G}}}}}}{{{{\times }}}}{{{{{\boldsymbol{E}}}}}}$ interaction or mediation. H AFR vs HIS, ${{{{{\boldsymbol{G}}}}}}{{{{\times }}}}{{{{{\boldsymbol{E}}}}}}$ interaction or mediation. I AFR vs EAS, no ${{{{{\boldsymbol{G}}}}}}{{{{\times }}}}{{{{{\boldsymbol{E}}}}}}$ interaction or mediation. J AFR vs EAS, ${{{{{\boldsymbol{G}}}}}}{{{{\times }}}}{{{{{\boldsymbol{E}}}}}}$ interaction or mediation. K HIS vs EAS, no ${{{{{\boldsymbol{G}}}}}}{{{{\times }}}}{{{{{\boldsymbol{E}}}}}}$ interaction or mediation. L HIS vs EAS, ${{{{{\boldsymbol{G}}}}}}{{{{\times }}}}{{{{{\boldsymbol{E}}}}}}$ interaction or mediation. The variants with no ${{{{{\boldsymbol{G}}}}}}{{{{\times }}}}{{{{{\boldsymbol{E}}}}}}$ interaction or mediation are those not included in Supplementary Data S3a. The variants with ${{{{{\boldsymbol{G}}}}}}{{{{\times }}}}{{{{{\boldsymbol{E}}}}}}$ interaction or mediation are those in Supplementary Data S3a. We only included independent variants. The shadow error bands represent the 95% confidence intervals. Clearly the variants without ${{{{{\boldsymbol{G}}}}}}{{{{\times }}}}{{{{{\boldsymbol{E}}}}}}$ interactions or mediations have substantially larger cross-population correlations than the variants with ${{{{{\boldsymbol{G}}}}}}{{{{\times }}}}{{{{{\boldsymbol{E}}}}}}$ interactions or mediations, suggesting that ${{{{{\boldsymbol{G}}}}}}{{{{\times }}}}{{{{{\boldsymbol{E}}}}}}$ interactions or mediations contribute the marginal effect size heterogeneity across populations. (European (EUR), African (AFR), Hispanics (HIS), Eastern Asian (EAS)).

Discussion

In this study, we utilize marginal effects from GWAS to search for $G\times E$. We conceptually demonstrated the deep connection between detecting $G\times E$ and MR for causal inference. Although ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$ tests for the combined effect of $G\times E$ and mediation, the two-step procedure of ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$ followed by ${T}_{{Direct}}$ in fact tests for $G\times E$, and its statistical power is much improved because of the following reasons: (1) the step 1 ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$ can increase power when a genetic variant has a mediation effect through the environmental factor. In this case, we expect a larger difference between the marginal effect and the main effect (Eq. (3)) than no mediation. (2) The difference between the marginal and main effect can further increase when the environmental distributions between GWAS and GWIS cohorts are different (Eq. (3) and Fig. 1D). (3) At the two-step procedure, multiple comparison burden is significantly alleviated because only significant variants survived at step 1 need to be examined. As demonstrated in this study, the two-step procedure identified 8 independent signals in comparing with two by the direct test in GWIS. This is also consistent with when comparing with the direct test in GWIS, the two-step procedure identified more $G\times E$ signals for each of the three lipid traits and four environmental factors (Supplementary Data S3b). Detecting $G\times E$ using direct tests can be biased by unmeasured confounders due to omitting covariates in the regression models³⁵, but the two-step procedure is robust because ${T}_{{MR\_GxE}}$ is not affected by confounders such as population structure. Considering the advantages of the two-step procedure, we view it as a complement rather than a replacement of the direct test. This perspective arises from the fact that the two-step test necessitates additional GWAS summary statistics and may be less powerful than the direct test in some situations (Fig. 1D).

Our study demonstrated that the current heritability estimates based on marginal effects could also include contributions from $G\times E$ and mediation through the corresponding environment factors (Fig. 4 and Supplementary Data S4). We excluded cigarette smoking- or alcohol drinking-associated variants identified from a large cigarette smoking and alcohol consumption GWAS of 1.2 million individuals²⁸ in our analysis, which mitigates the potential mediation contribution in the ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$ analysis. However, among the 15 loci identified by ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$, only three were not reported in the much larger recent cigarette smoking and alcohol consumption GWAS of 3.4 million individuals, suggesting mediation through cigarette smoking and/or alcohol consumption is still present but with modest effects. Among the six $G\times E$ variants identified, 4 are associated with either cigarette smoking or alcohol drinking, suggesting that the $G\times E$ variants are also likely to be mediated through E and the mediation improves power to detect $G\times E$. Furthermore, we demonstrated that the current SNP heritability estimates based on marginal effects could also include significant contributions from $G\times E$ and mediation through the corresponding environment factors for LDL-C and TG (Fig. 4 and Supplementary Data S4). We did not observe significant contributions from $G\times E$ and mediation to heritability for HDL-C, potentially attributable to the relatively small sample sizes in our GWIS. Since the LDSC regression³⁴ cannot be used to estimate $G\times E$ heritability, our estimates reflect the low bound of the interaction and environmentally mediated heritability. We therefore suggest that the current SNP heritability estimates based on the marginal genetic effects be called marginal SNP heritability, to differentiate it from narrow-sense heritability³⁶ that is defined by additive genetic actions without the inclusion of $G\times E$ or mediation contributions. We believe this differentiation is important for correctly interpreting the current heritability estimates and understanding the genetic architecture of complex traits.

The 5 (6 independent signals) replicated loci interacting with cigarette smoking and alcohol consumption contain genes that are enriched in liver tissue, possibly reflecting the effect of alcohol drinking on aspartate amino transferase, alanine aminotransferase and γ-glutamyl transferase activities via the actions of numerous ingredients that alter the activities of enzymes found in the liver³⁷. Among them, the interaction between alcohol consumption and cholesteryl ester transfer protein (CETP) has been reported for HDL-C and coronary artery disease^38,39,40. The interaction between alcohol consumption and APOE on LDL-C has also been reported in a Mediterranean Spanish population⁴¹, while the interactions between APOA5 and cigarette smoking and alcohol drinking status associated with elevated TG and reduced HDL-C were observed in the Chinese and Korean populations^42,43. However, our study is the only well-powered study demonstrating significant evidence at the genome-wide level and the interaction loci are replicable. SMARCA4 was reported to be associated with LDL-C in the lipids GWAS in Africans⁴⁴ but not in the recent largest lipids GWAS which is predominantly European ancestry³. Overall, the marginal effect sizes of the variants are less correlated across populations for the variants with than without $G\times E$ interaction or mediation (Fig. 5), empirically verified that $G\times E$ and mediation contribute to marginal effect differences across different populations⁴⁵. We expect that including $G\times E$ interactions should improve polygenetic risk score prediction across populations.

It is well known that causal effect estimate in MR framework can be biased when the three IV assumptions are violated. However, our goal is to detect $G\times E$ rather than to estimate the causal effect. Detecting $G\times E$ based on MR is less likely to be biased for these reasons: (1) the effect sizes of IVs on the pseudo exposure are all highly significant in GWAS, which represent strong IVs. (2) It is less likely to have a confounding effect between a trait and its pseudo exposure, i.e., a polygenic score. (3) The iterative Mendelian randomization and pleiotropy test is a powerful method to detect pleiotropy when the two IV conditions are satisfied²³, in particular, it is expected that most of the IVs are not interacted with E. (4) Although the causal effect estimate can be affected if population structure is not well corrected, testing $G\times E$ by ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$ is not. The reason is that ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$ can be viewed as a weighted linear regression of the effect size of GWAS on the effect size of GWIS followed by searching the outlies of variants departing from the regression line. While the regression line (equivalent to causal effect estimate in MR) can be affected by population structure, the outlie detection is not.

In summary, our $G\times E$ approach is powerful and able to detect genetic loci interacting with environments that account for significant phenotypic variability. Our findings indicate that the contribution of $G\times E$ in lipids is not ignorable. Our study only focuses on the interactions of genes with cigarette smoking and alcohol consumption in lipids. The cumulative interaction contribution with many environmental factors can even be greater. Detecting individual genetic loci with environmental interactions facilitates a better understanding of the genetic architecture of complex traits and can improve phenotype prediction.

Methods

Summary statistics data

The marginal summary statistics of high-density lipoprotein cholesterol (HDL-C), low-density lipoprotein cholesterol (LDL-C), and triglycerides (TG) from the Global Lipids Genetics Consortium study (GLGC, n = 1.65 M)³ were downloaded at http://csg.sph.umich.edu/willer/public/glgc-lipids2021.

GLGC consists of GWAS results from 1.65 M subjects representing five genetic ancestry groups: European (N = 1.32 M); African or admixed African (N = 99 K); East Asian (N = 146 K); Hispanic (N = 48 K); and South Asian (N = 41 K). We did not perform South Asian specific analysis because there was no corresponding GWIS in the Cohorts for Heart and Aging Research in Genetic Epidemiology (CHARGE) consortium. The GWIS summary statistics from CHARGE gene-lifestyle (GLI) working group in this study are available via dbGaP (accession number phs000930). The CHARGE GWIS consists of 60 GWIS summary datasets: (LDL-C, HDL-C, and TG)-current smoking, (LDL-C, HDL-C, and TG)-ever smoking, (LDL-C, HDL-C, and TG)-current alcohol drinking, and (LDL-C, HDL-C, and TG)-regular alcohol drinking, for European, African or Admixed African, East Asian, Hispanic and multi-ancestry.

QCs for performing ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$ analysis

To perform MR analysis, we aligned the GWAS summary statistics HDL-C, LDL-C, and TG from the GLGL with the corresponding GWIS summary statistics from the CHARGE gene-lifestyle consortium. We flipped the effect size if the corresponding reference allele did not match. We dropped a genetic variant if the two alleles were either {A, T} or {C, G}. We also excluded any variants with minor allele frequency (MAF) difference larger than 0.15 between GWAS and GWIS study. If multiple variants fall on the same chromosome position, we required the matched variants with MAF difference less than 0.01. We further excluded any variants with the effective sample size in GLGC trans-ethics or European less than 100 K and the other populations (African, Hispanic, East Asian) less than 30 K. To reduce the effect by mediations through the smoking and alcohol drinking, we excluded all loci with P-value < 5E−7 identified by the GWAS of smoking status or alcohol drinking²⁸.

${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$ analysis

To perform ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$, we applied the Mendelian randomization (MR) software IMRP²³ to estimate the causal effect by considering the main effect sizes from the GWIS of the CHARGE gene-lifestyle consortium as the exposure effects, and the marginal effects from the GLGC as the outcome effects, respectively. To identify instrument variables, we first selected all the variants with the P-value < 5E−8 after GC-correction in the GLGC, and then pruned them using the window size 500 KB and ${r}^{2}$ value 0.1 by the Plink software⁴⁶. We standardized the effect sizes as in²⁷. IMRP requires the input of the correlation coefficient to account for the effect of sample overlapping between GWAS and GWIS cohorts and this correlation was calculated based on the unsignificant variants (P-value > 0.05) across the genome. After estimate the causal effect, we performed ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}},$ which is equivalent to the pleiotropy test in the IMRP, to all the genetic variants across the genome.

Independent locus definition

Independent loci were defined as the regions within 1 Mb of the most significant variants by the ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$ test. Independent signals were defined as the variants in a locus with ${r}^{2}$  <  0.1. The 1000 G data was used as the reference genetic data for LD calculation.

Choosing independent variants for replication in UK Biobank

By applying ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$, we observed that 15 genome-wide significant loci consisting of 17 independent signals (P-value < 5E−8), including 4 and 5 loci for LDL-C, 7 and 5 loci for HDL-C, and 5 and 6 loci for TG, interacted with alcohol drinking and cigarette smoking, respectively (Supplementary Data S3a). At a locus with the ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$ significant (P-value < 5E−8) for a lipid trait (LDL-C, HDL-C, or TG) and environment (smoking or alcohol drinking), we searched the variant with the smallest P-value of the direct test ${T}_{{Direct}}$ among the significant variants by the ${T}_{{MR}{{{{{\rm{\_}}}}}}{GxE}}$. The variants with ${T}_{{Direct}}$ P-value < 7.35E−4, which correct for the 17 tests and 4 environmental factors, were considered as significant for $G\times E$ interaction (two-step procedure). We obtained 8 independent variants in 6 loci among these 17 independent signals survived the threshold P-value = 7.35E−4 and these variants were further tested for the replication of the interaction effects in UK Biobank using ${T}_{{Direct}}$ test.

LD score regression

We applied the LD score regression to estimate heritability contributed by $G\times E$ interaction and mediation through the environment factor $E$. We estimated heritability by combining all chromosomes rather than chromosome specifically. We used the R package bigsnpr⁴⁷ to estimate LD scores in the corresponding populations from 1000 G reference data with default settings.

Functional mapping and annotation

We performed overall enrichment tests using the residual ${\hat{\alpha }}_{j}-{\hat{\beta }}_{j}\hat{\theta }$ as the effect size and ${se}({\hat{\alpha }}_{j}-{\hat{\beta }}_{j}\hat{\theta })$ as the corresponding standard error. We used MAGMA³⁰ (Multi-marker Analysis of GenoMic Annotation) and DEPICT⁴⁸ (Data-driven Expression Prioritized Integration for Complex Traits) to identify tissues and cells that are highly expressed at genes within the $G\times E$ loci. We also used DEPICT to test for enrichment in gene sets associated with gene ontology (GO) ontologies, mouse knockout phenotypes and protein-protein interaction networks. In addition, we reported significant enrichments with a false discovery rate 0.05. Analysis was done using the online platform FUMA GWAS.

Colocalization

We performed colocalization analysis by using the software ezQTL (https://analysistools.cancer.gov/ezqtl/#/home). We chose the public genotype-tissue expression (GTEx) v7 with eQTL³² as the QTL data and chose the public European reference panels for calculating the LD data. We performed colocalization analysis between GWIS and QTL results within a locus using eCAVIAR (eQTL and GWAS Causal Variant Identification in Associated Regions)⁴⁹, where the Colocalization Posterior Probability (CLPP) is used to describe the significance level of colocalization. We only recorded colocalization with CLPP > 0.01, as suggested by the authors of eCAVIAR.

UK Biobank individual level data for replication

The UK Biobank (UKBB)⁵⁰ individual-level data used for replications were available through Application ID: 81097. Quality Controls Participants in the UKBB were genotyped using a custom Affymetrix UK Biobank Axiom array⁵¹. Genotypes were imputed by the UKBB using the Haplotype Reference Consortium reference panel⁵² with imputation ${r}^{2}$ value greater than 0.3. Related individuals with pairwise kinship coefficient greater than 0.0884 (suggested by UKBB) were removed from analysis, resulting in N = 445,424 individuals of European, African, and East Asian ancestries. The principal components were calculated by UKBB with genotype data within each ancestry to account for population structure. These data were independent of GLI cohorts and consisted of European, African, and Asian individuals (race determined using UKBB field ID 21000-0.0) in UKBB who were unrelated (genetic kinship coefficient less than 0.0884; 22021-0.0). Linear regression model (1) in main text was performed. Covariates included age at assessment (21003-0.0), age², sex (31-0.0), the first 10 PCs (22009-0.1 to 22009-0.10), and environment exposure, a genetic variant and their interaction. Environmental exposures included ever/never smoking status (20116-0.0), current/non-current smoking status (20116-0.0), and alcohol intake frequency (1558-0.0).

Analogous to the $G\times E$ analysis in ref. ¹⁷, HDL-C (30760-0.0) and TG (30870-0.0) measurements were natural log transformed and LDL-C measurements (30780-0.0) were converted from mmol/L to mg/dl then multiplied by a factor of 0.7 if there was a history of lipid-lowering medication (6177-0.0) present. LDL-C measurements were therefore considered no medication if there were missing values for medication history. This introduced missing values in LDL-C for 248,419 individuals.

Theoretical properties of ${T}_{{MR}{GxE}}$

In MR analysis, the instrumental variables are independent and are genome-wide significant variants selected from GWAS. Let ${\hat{\beta }}_{1,j},{\hat{\beta }}_{2,j},{\hat{\beta }}_{3,j}$ and ${\hat{\alpha }}_{j}$, $j=1,\ldots,m$, be the corresponding effect size estimates in GWIS (model (1) and GWAS (model 2)) for the m instrument variables.

The causal effect $\theta$ of the inverse variance weighted (IVW) is estimated by

$$\hat{\theta }={{\arg }}\mathop{\min }\limits_{\theta }\left\{\frac{1}{m}\mathop{\sum }\limits_{j=1}^{m}\frac{{\left({\hat{\alpha }}_{j}-{\hat{\beta }}_{1,j}\theta \right)}^{2}}{{{{{\mathrm{var}}}}}\left({\hat{\alpha }}_{j}\right)}\right\}.$$

(4)

It is much simply to work on $\hat{\theta }$ by standardizing the IVs and this procedure does not change the conclusion. Thus, we let ${\sigma }_{G,j}^{2}=1$, j = 1$,$ …, m, in both GWAS and GWIS data. Further, we let the phenotype residue variance ${\sigma }^{2}=1.$ By equation (S15) in Supplementary Note, we have ${{{{\mathrm{var}}}}}({\hat{\alpha }}_{{{{{{\rm{j}}}}}}})={n}_{1}^{-1},j=1,\ldots m,{{{{{\rm{and}}}}}}$ $\hat{\theta }=\frac{\mathop{\sum }\limits_{j=1}^{m}{\hat{\alpha }}_{j}{\hat{\beta }}_{1,j}}{\mathop{\sum }\limits_{j=1}^{m}{\hat{\beta }}_{1,j}^{2}}.$

Since only the variants without either $G\times E$ interaction or mediation are valid in the MR analysis, we assume $\rho=0$ (no mediation) and ${\beta }_{3,{{{{{\rm{j}}}}}}}=0$ (no interaction). We have

$${\hat{\alpha }}_{j}={\hat{\beta }}_{1,{{{{{\rm{j}}}}}}}+{\mu }_{E1}{\hat{\beta }}_{3,{{{{{\rm{j}}}}}}}$$

(5)

By applying the Slutsky’s theorem, and let ${\beta }_{3,{{{{{\rm{j}}}}}}}=0$, we have:

$$E\left(\hat{\theta }\right)=\frac{\frac{1}{m}\mathop{\sum }\limits_{j=1}^{m}E({\hat{\alpha }}_{j}{\hat{\beta }}_{1,j})}{\frac{1}{m}\mathop{\sum }\limits_{j=1}^{m}E({\hat{\beta }}_{1,j}^{2})}=\frac{\frac{1}{m}\mathop{\sum }\limits_{j=1}^{m}{\beta }_{1,j}^{2}+\frac{{n}_{0}}{{n}_{1}{n}_{2}}\left(1+\frac{{\mu }_{E0}^{2}}{{\sigma }_{E0}^{2}}\right)}{\frac{1}{m}\mathop{\sum }\limits_{j=1}^{m}{\beta }_{1,j}^{2}+\frac{1}{{n}_{2}}\left(1+\frac{{\mu }_{E2}^{2}}{{\sigma }_{E2}^{2}}\right)}.$$

(6)

Because ${\sigma }_{G,j}^{2}=1$, $\frac{1}{m}\mathop{\sum }\limits_{j=1}^{m}{\beta }_{1,j}^{2}$ is the average phenotypic variance accounted by an IV. Define ${\sigma }_{\beta }^{2}=\frac{1}{m}\mathop{\sum }\limits_{j=1}^{m}{\beta }_{1,j}^{2}$, we have:

$$E\left(\hat{\theta }\right)=\frac{{\sigma }_{\beta }^{2}+\frac{{n}_{0}}{{n}_{1}{n}_{2}}\left(1+\frac{{\mu }_{E0}^{2}}{{\sigma }_{E0}^{2}}\right)}{{\sigma }_{\beta }^{2}+\frac{1}{{n}_{2}}\left(1+\frac{{\mu }_{E2}^{2}}{{\sigma }_{E2}^{2}}\right)},$$

(7)

which converges to 1 when n₁ and n₂→∞. However, when ${\sigma }_{\beta }^{2}$ is small (weak instrument in MR analysis), the converge of $E(\hat{\theta })$ to 1 is slow. We also note that $E(\hat{\theta })\le 1$.

Simulation settings without medication contribution (Fig. 2A–D, Supplementary Figs. 6–9)

For ith individual, we generated m = 102 number of independent variants through for $j=1,\ldots,m$ by ${G}_{{ij}}^{*}\sim {{\mbox{Binom}}}(2,{p}_{j})$, where p_j∼unifom (0.05,0.5). We standardized genotypes by ${G}_{{ij}}=\frac{{G}_{{ij}}^{*}}{2{p}_{j}\left(1-{p}_{j}\right)}.$ For the environment factor in the GWAS model, we generated ${E}_{i1}\sim {{{{{\mathscr{N}}}}}}\left({\mu }_{E1},1\right)$. For the environment factor in the GWIS model, we generated ${E}_{i2}\sim {{{{{\mathscr{N}}}}}}\left({\mu }_{E2},1\right)$. For the samples overlapped between the GWAS and GWIS, we generated their environment values through ${{{{{\mathscr{N}}}}}}\left({\mu }_{E2},1\right).$ We varied the values of ${\mu }_{E1}$, ${\mu }_{E2}$ and the proportion of overlapped samples.

The main effect size of the jth variant was generate by ${\beta }_{1j}{{{{{\mathscr{\sim }}}}}}{{{{{\mathscr{N}}}}}}(0,{\sigma }_{\beta }^{2}),$ where ${\sigma }_{\beta }^{2}$ is the trait variance accounted for by the IVs. For the first variant, we added its interaction effect with E. The phenotype ${Y}_{i}$ by generated by

$${Y}_{i}=\mathop{\sum }\limits_{j=1}^{m}{G}_{{ij}}{\beta }_{1j}+0.1{E}_{i}+0.05\left({G}_{i1}\times {E}_{i1}\right)+{\epsilon }_{i},$$

(8)

where ${\epsilon }_{i}{{{{{\mathscr{\sim }}}}}}{{{{{\mathscr{N}}}}}}\left(0,{\sigma }^{2}\right)$. The causal effect $\theta$ was estimated using the last 100 variants as the IVs. The power and type I error rate for ${T}_{{Direct}}$ and ${T}_{{MR\; GxE}}$ were calculated based on the first and second variants, respectively.

Simulation settings without medication contribution (Fig. 2E–F), Supplementary Fig. 10

We generated $20$ independent variants by ${G}_{j}$ ${{\mbox{Binom}}}\left(2,0.3\right)$ and standardized it but without mean correction. We simulated environment E according to mediation present or not present. If no mediation, $E$ is generated from ${{{{{\mathscr{N}}}}}}\left(1,1\right)$. If there is mediation, $E \sim 0.05G+{{{{{\mathscr{N}}}}}}\left(2,0.9975\right)$, or $G$ contributes 0.25% variation of E. The phenotype is generated according to the following models:

• 1. No mediation and no interaction: $Y \sim 0.1G+\gamma E+N(0,10)$, where $E \sim {{{{{\mathscr{N}}}}}}\left(1,1\right)$

• 2. Mediation but no interaction: $Y \sim 0.1G+\gamma E+N(0,10)$, where $E \sim 0.05G+{{{{{\mathscr{N}}}}}}\left(1,0.9975\right)$.

• 3. Mediation and interaction: $Y \sim 0.1G+\gamma E+0.1*G*E+N(0,10)$, where $E \sim$$0.05G+{{{{{\mathscr{N}}}}}}\left(1,0.9975\right)$.

We let $\gamma$ take values of 1 and $\sqrt{5}$. We also simulated data with environment mean 0.5 (Supplementary Fig. 10). We first simulated ${n}_{2}={{{{\mathrm{20,000}}}}}$ subjects for GWIS cohort (or main effect estimation). The sample size for marginal effect estimation varied from ${n}_{1}={{{{\mathrm{20,000}}}}}$ to 300,000, with the 20,000 subjects in GWIS cohort was always included. For the non-overlapped subjects, we let the environment mean to be 1.5 times of the environment mean in GWIS cohort. The type I error and power for ${T}_{{Direct}}$ and ${T}_{{MR\_GxE}}$ were calculate by correcting for 20 tests using the Bonferroni correction. For the two-step procedure, we first applied ${T}_{{MR\_GxE}}$ and Bonferroni correction. The variants survived after ${T}_{{MR\_GxE}}$ were further tested by ${T}_{{Direct}}$ and Bonferroni correction was further applied.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.