A new 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 CHARGE Gene-Lifestyle Interactions Working Group has been spearheading efforts to investigate G X 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 identified and confirmed 5 loci (6 independent signals) interacting with either cigarette smoking or alcohol consumption for serum lipids, and empirically demonstrated 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 contribution ranges from 1.76–14.05% of SNP heritability of serum lipids in Cross-Population data. Our study improves the understanding of the genetic architecture and environmental contributions to 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 CHARGE Gene-Lifestyle Interactions Working Group has been spearheading efforts to investigate  ×  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 identified and confirmed 5 loci (6 independent signals) interacting with either cigarette smoking or alcohol consumption for serum lipids, and empirically demonstrated 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 contribution ranges from 1.76% to 14.05% of SNP heritability of serum lipids in Cross-Population data.Our study improves the understanding of the genetic architecture and environmental contributions to complex traits.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  ×  contribution to heritability because of different theoretical models and assumptions 10,11 .As pointed out in 12 , 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 5 , 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 models 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  ×  14 .
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  ×  at the genome-wide level are seldom considered as primary analyses.Instead, the joint evidence of main genetic and  ×  effects, in addition to the  ×  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  ×  alone [16][17][18][19] .The joint test limits our ability to determine to what degree the currently identified loci reflect evidence for  ×  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 MR approaches have been used to explain  ×  22 and assess risk factor interactions 23 , the underlying connection between testing pleiotropic variants in the MR framework and the detection of  ×  is currently unclear.Here, we conceptually connect  ×  with the MR framework, illuminate their similarities and demonstrate that the test of horizontal pleiotropy in MR 24,25 can be used for detecting  × .Based on this principle, one can identify novel  ×  using existing summary statistics without needing costly and time-consuming new analyses from all cohorts.We applied this approach to the summary statistics from the Global Lipids Genetics Consortium study (GLGC, n=1.65M) 3 and the summary statistics in the interaction analysis with cigarette smoking and alcohol drinking in the CHARGE GLI working group 17 , 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).

Testing 𝑮 × 𝑬 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: where  1 ,  2 and  3 correspond to the 'main' effect of  (in the presence of ), the main effect of  and the interaction effect of  × , respectively, and  is a random noise.Here , , and  ×  represent a genotype value, environmental factor, and their interaction, respectively.For simplicity, we do not include any covariates, but this does not affect the general conclusion.The interaction effect is evaluated by the direct test statistic   =  ̂3 2 ( ̂3) ⁄ , where  ̂3 refers to the estimate from the regression model (1).Theoretical work indicates that the test statistics for the main effect  1 = 0 and the interaction effect  3 = 0 are correlated, with the correlation coefficient equal to −  /√  2 +   2 , where   and   2 are the mean and variance of the environmental factor in the data 14 .However, the power of the direct test is usually low because of the collinearity between  and  ×  which induces a covariance between the estimates of  1 and  3 .This covariance produces uncertainty (i.e.larger standard error) which by itself reduces power for testing either  1 or  3 , even if the underlying true model includes  ×  alone (i.e.,  1 = 0 and  3 ≠ 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., where we refer to  as the 'marginal' effect from a GWAS (in the absence of ) to differentiate it from the main effect  1 in model (1).We show that ) 3 , where  is the mediation contribution of  through ,  1 ,  1 , and  1 represent the environmental mean, standard deviation, and genotype standard deviation in GWAS data, respectively, suggesting that testing the hypothesis H 0 :  −  1 = 0 for the difference between the marginal and main effects is equivalent to testing for the combined effect of  ×  and mediation, and further reduces to testing for the  ×  when  and  are independent (i.e.,  = 0, Supplementary Note).This hypothesis can be tested by the statistic   = ( ̂−  ̂1) 2 /var( ̂−  ̂1), where  ̂,  ̂1, and their corresponding standard errors are estimated from the GWAS and GWIS analyses, respectively.In fact,   and   are equivalent when GWAS and GWIS are performed in the same data.We verified this property using real data analysis in the GLI studies 17 , from which the summary statistics of the marginal, main, and interaction effects are available and the marginal effect was obtained after adjusting for  (Fig S1).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 26,27 .Therefore, directly using   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  ̂ and the main effect estimate  ̂1 have a linear relationship, where  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 (Online Methods).The genetic variants with no  ×  and no mediation will fall on the regression line but the variants with  ×  or mediation will depart from this line.This pattern will not be impacted by the systemic variation across studies.Therefore, we search the genetic variants that depart from this regression line to test the combined effect of  ×  and mediation, providing  can be correctly estimated.This idea is conceptually the same as the MR framework when we introduce a pseudo exposure  ̃, representing a polygenic score of the trait (Fig 1).We then estimate the causal effect  of the pseudo exposure  ̃ on trait  in the MR framework and the  ×  effect or mediation through  is tested in the same way as testing for horizontally pleiotropic variants 24 .In doing so, we first select a set of independent variants associated with trait  and perform the inverse variance weighted analysis to estimate , denoted as  ̂.Second, we test the  ×  or mediation of a genetic variant through  by the 2 .This test can be performed by the iterative Mendelian randomization and pleiotropy (IMRP) approach 24,28 .The statistic  _ is an asymptotically unbiased test for testing the combined effect of  ×  and mediation through  (Supplementary Note).Two-step procedure for Testing  ×  Note that  _ likely tests for the combined effect of  ×  and mediation unless  and  are independent (i.e.,  = 0).To test for  × , we propose a two-step procedure by applying  _ to screening the whole genome and then performing   on the variants surviving the  _ screen.This two-step procedure can increase power at the screening step when there is interaction and mediation and increases power at the testing step by substantially reducing the multiple comparison burden. _ and   are not independent (Supplementary Note), and therefore, the variants detected by the two-step procedure could still reflect the contribution of mediation and  × .To mitigate this problem, we can exclude the variants identified through GWAS of , which could have large mediation effects.

Type I error rate and power of 𝑻 𝑴𝑹_𝑮𝒙𝑬 and the two-stage procedure
We first performed a series of simulations to investigate the type-I error rate and power of  _ in the absence of mediation ( = 0).In simulations we observed that ( ̂) is close to 1 and the estimate  ̂ converges to 1 when sample size increases, which is expected by our theoretical prediction (Fig. 2A and Fig. S2).The interaction effect  3 can be estimated directly by model (1)  or by ( ̂−  ̂1 ̂)/  when  = 0. We observed that both ways are unbiased (Fig, 2B), although the standard error of ( ̂−  ̂1 ̂)/  is affected by the environmental means in GWAS and GWIS.When no mediation is present, the type-I error rates for both  _ and the direct test are well controlled (Fig. 2C and Fig. S2c)).The power of  _ depends on multiple parameters, including   and allele frequency in GWAS and GWIS and is less powerful than   when the environmental mean in GWAS is lower (Fig. 2D and Fig. S2d).Additional simulations for the estimates of  ̂, interaction effect ( ̂−  ̂1 ̂)/  , type-I error rate and power are presented in Fig. S3-S5.
We next investigated the performance of   ,  _ 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 Fig. S6).When mediation is present, the type-I error rate was still well controlled, although inflation can be observed for the two-step procedure and  _ when  contributes to 5% of the outcome variation and the samples between GWAS and GWIS are completely overlapped (Fig S6).This inflation was caused by mediation and quickly disappeared when the overlapping rate between GWAS and GWIS subjects was reduced.The statistical power of  _ and the two-step procedure for testing  ×  was much more improved than   when mediation was present (Fig. 2F and Fig. S6).

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.65M) and the CHARGE GLI (n=134K).To mitigate the effects of mediation through cigarette smoking or alcohol drinking, we excluded all loci with P-value <5×10 -7 reported in the early GWAS of cigarette smoking status or alcohol drinking 29 , which represent relatively large effect sizes of variants on cigarette smoking and alcohol drinking.We observed that  ̂ ranged from 0.92-1.33,0.95-1.62,0.83-1.25,0.87-1.37,and 0.95-1.28for European, African, Asian, Hispanic, and Cross-population data, respectively (Table S1).The departure of  ̂ 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  _ , with the genomic control  values being close to 1 (range 0.93-1.05,Table S2).
Using  _ to screen the genome, we observed 15 genome-wide significant loci consisting of 17 independent signals (P<5×10 -8 ,  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, Table S3).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 30 , suggesting the contribution of both  ×  and mediation.Since we already excluded the cigarette smoking and alcohol drinking variants identified from a relatively smaller study 29 , these detected variants should represent modest mediation effects.Locus specific plots of all significantly associated loci are presented in Fig. S7, 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  ×  or mediation loci are clearly departing from most of the lipids-associated variants (Fig. 3D-F, J-L).
The population-specific  _ results are presented in Fig. S8, which are also consistent with the Cross-population results, although the main contribution comes from the European population.
By applying the two-stage procedure, we observed that 8 of the 17 independent signals are significant when using the direct test   after correcting for the 17 tests and 4 environmental factors (Table 1, P<7.35×10 -4 ).Tissue enrichment analyses using GWAS-based pathway analysis tools, MAGMA 31 and FUMA 32 , suggest that these loci are enriched in liver, hippocampus, small intestine, and stomach tissues (Fig. S9).Multiple loci were colocalized with expression quantitative trait loci (eQTLs) in the corresponding liver, lung, and blood tissues in the genotypetissue expression database (GTEx) 33 (Fig. S10).

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   calculated in the UKBB is independent of  _ , so are the   test statistics calculated in UKBB and CHARGE GLI, thus qualifying this as an independent replication (Supplementary Note).Six of the 8 signals were replicated in the UKBB after adjusting for 32 tests (P<1.56×10 - ), 5 of which were genome wide significant by the   test in combined CHARGE GLI and UK Biobank data (Table 1).All 8 independent signals had the same interaction direction in CHARGE GLI and UKBB except LPL, which is not significant in UKBB (Table S3).The CETP and SMARCA4 loci were the only two loci with no reported mediation evidence through either cigarette smoking or alcohol drinking.
We next aimed to understand if the interaction evidence observed in this study had an alternative explanation 34 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 (Fig S11), 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.

𝑮 × 𝑬 Interaction and Mediation to SNP Heritability
Since  ̂−  ̂1 ̂ refers to the combined interaction and mediation contribution to the marginal effect, we can use  ̂−  ̂1 ̂ to estimate the heritability contributed by the interaction and mediation through the LD score (LDSC) regression 35 .Note that this heritability is a lower bound of the phenotype variance contributed by the  ×  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.We observed significant interaction and mediation heritability (P<1.93×10-3 ) with cigarette smoking and alcohol consumption to lipids traits, ranging from 1.76% (LDL-C, regular drinking) to 14.05% (TG, regular drinking) of SNP heritability estimated by the marginal effects (Fig. 4, Table S4).

𝑮 × 𝑬 Interaction and Mediation to heterogeneity of genetic effect sizes across populations
As noted in equation ( 3), the marginal effect estimate of a genetic variant in GWAS consists of the  ×  and mediation contribution when the  ×  or mediation occurs.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  ×  interaction or mediation than without.We calculated the marginal effect size correlations between European, African, Hispanic, and Eastern Asian populations for these variants reported in Graham et al 3 after excluding the variants in Table S3 where their  ×  interactions or mediations were observed in this study.Similarly, we calculated the marginal effect size correlations for the variants in Table S3 using the effect sizes reported in Graham et al 3 .We compared the correlations and observed a median of 24.4% drop of the cross-population correlation coefficient (Fig 5), strongly suggesting that  ×  interactions or mediations contribute to the marginal effect size hetergeniety across populations.

Discussion
To the best of our knowledge, this study is the first to utilize marginal effects from GWAS to search for  × .We conceptually demonstrated the deep connection between detecting  ×  and MR for causal inference.Although  _ tests for the combined effect of  ×  and mediation, the two-step procedure of  _ followed by   in fact tests for  × , and its statistical power is much improved when mediation is present and type I error rate is reasonably controlled.Detecting  ×  using direct tests can be biased by unmeasured confounders due to omitting covariates in the regression models 36 , but the new two-stage procedure is robust because  _ is not affected by confounders such as population structure.
Our study demonstrated that the current heritability estimates based on marginal effects also include significant contributions from  ×  and mediation through the corresponding environment factors (Fig. 4 and Table S4).We excluded cigarette smoking-or alcohol drinkingassociated variants identified from a large cigarette smoking and alcohol consumption GWAS of 1.2 million individuals 29 in our analysis, which mitigates the potential mediation contribution in the  _ analysis.However, among the 15 loci identified by  _ , 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  ×  variants identified, 4 are associated with either cigarette smoking or alcohol drinking, suggesting that the  ×  variants are also likely to be mediated through E and the mediation improves power to detect  × .Furthermore, we demonstrated that the current SNP heritability estimates based on marginal effects also include significant contributions from  ×  and mediation through the corresponding environment factors.We therefore suggest that the current SNP heritability estimates based on the marginal genetic effects be called marginal or broad-sense SNP heritability, to differentiate it from narrowsense heritability 37 that is defined by additive genetic actions without the inclusion of  ×  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 replicated loci (6 independent signals) interacting with cigarette smoking or 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 38 .Among them, the interaction between alcohol consumption and cholesteryl ester transfer protein (CETP) has been reported for HDL-C and coronary artery disease [39][40][41] .The interaction between alcohol consumption and APOE on LDL-C has also been reported in a Mediterranean Spanish population 42 , 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 43,44 .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 45 but not in the recent largest lipids GWAS which is predominantly European ancestry 3 .Overall, the marginal effect sizes of the variants are less correlated across populations for the variants with  ×  interaction or mediation than without (Fig 5), empirically verifying that  ×  and mediation contribute to marginal effect differences across different populations 46 .We expect that including  ×  interactions should improve polygenetic risk score prediction across populations.
It is well known that causal effect estimates in the MR framework can be biased when any of the three IV assumptions are violated.However, the MR-based  ×  approach 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 24 and, in particular, it is expected that most of the IVs do not interacted with the environmental factor E.
In summary, our novel  ×  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  ×  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.Instead of an explicit exposure, we create a pseudo exposure  ̃, which can be viewed as the total contribution of genetic variants to trait Y.The genetic variants associated with the pseudo exposure  ̃ but not through either the environment E or  ×  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  ̃ and Y, and thus have similar effects as  ×  and cannot be distinguished from  × .Thus, testing the combined effect of interaction and mediation is conceptually equivalent with testing the horizontally pleiotropic effect in the MR framework.Right panel: like the horizontal pleiotropic variants in the MR framework,  ×  variants depart from the regression line and can be detected assuming no medication.

Fig. 2. Simulation performance of 𝑻 𝑴𝑹_𝑮𝑿𝑬 and the two-step procedure. A-D:
No medication was present.The simulation details were described in Online Methods. A. ( ̂) is close to 1 as expected.B. The direct estimate of  3 in GWIS or by ( ̂−  ̂1 ̂)/  through MR analysis are both 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  _ and the direct test for different main and interaction effect directions.Both  _ and the direct test maintain the type I error rate well.D. Power comparison between  _ 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 Online    S4.Fig. 5. Cross population comparison of the LDL-C, HDL-C and TG marginal effect sizes of the variants reported in Graham et al 3 and the independent variants in Table S3 where their  ×  interactions or mediations were observed in this study.The correlations in the first and the third columns represent the variants reported in Graham et al while the correlations in the second and Online 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.65M) 3 were downloaded at http://csg.sph.umich.edu/willer/public/glgc-lipids2021.GLGC consists of GWAS results from 1.65M subjects representing five genetic ancestry groups: European (N = 1.32M);African or admixed African (N = 99K); East Asian (N = 146K); Hispanic (N = 48K); and South Asian (N = 41K).We did not perform South Asian specific analysis because there was no corresponding GWIS in the CHARGE consortium.The GWIS summary statistics from CHARGE GLI working group in this study are available via dbGaP (accession number phs000930).
QCs for performing  _ analysis.
To perform MR analysis, we merged the GWAS summary statistics for HDL-C, LDL-C, and TG from the GLGC with the corresponding GWIS summary statistics from the CHARGE GLI consortium.We flipped the effect size if the corresponding reference allele did not match.We dropped the genetic variant if its 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 GLGC and CHARGE GLI studies.If multiple variants fell 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 100K and the other populations (African, Hispanic, East Asian) less than 30K.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 29 .

𝑻 𝑴𝑹_𝑮𝒙𝑬 analysis
To perform  _ , we applied the Mendelian randomization (MR) software IMRP 24 to estimate the causal effect by considering the main effect sizes from the GWIS of the CHARGE genelifestyle consortium as the exposure effects, and the marginal effects from the GLGC as the outcome effects, respectively.To identify instrumental 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  2 value 0.1 by the Plink software 47 .We standardized the effect sizes as in 28 .
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 variants with P-value > 0.05) across the genome 48 .After estimating the causal effect, we performed  _ , 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 1Mb of the most significant variants by the  _ test.Independent signals were defined as the variants in a locus with  2 < 0.1.The 1000G data was used as the reference genetic data for LD calculation.

Choosing independent variants for replication in UK Biobank
By applying  _ , 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 (Table S3).At a locus with significant  _ test 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   among the significant variants by the  _ .The variants with   Pvalue < 7.35E-4, which correct for the 17 tests and 4 environmental factors, were considered as significant for  ×  interaction in the two-step procedure.We observed 8 independent signals in 6 loci among all 17 independent signals surviving the threshold P-value = 7.35E-4.These variants were further tested for the replication of the interaction effects in UK Biobank using   test.

LD Score Regression
We applied the LD score regression 35 to estimate heritability contributed by  ×  interaction and mediation through the environment factor . To account for potential heterogeneity, we first estimated chromosome specific heritability and then summed across chromosomes.We used the R package bigsnpr 49 to estimate LD scores with the 1000G Phase 3 reference data and then estimated the chromosome-specific heritability with default settings.We observed that some of the chromosome specific heritability estimates were negative.We only summed the non-negative chromosome specific heritability estimates.

Functional Mapping and Annotation.
We performed overall enrichment tests using the residual  ̂ −  ̂ ̂ as the effect size and ( ̂ −  ̂ ̂) as the corresponding standard error.We used MAGMA 31 (Multi-marker Analysis of GenoMic Annotation) and DEPICT 50 (Data-driven Expression Prioritized Integration for Complex Traits) to identify tissues and cells that are highly expressed at genes within the  ×  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 51 .We chose the public genotype-tissue expression (GTEx) v7 with eQTL 33 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) 52 , where the Colocalization Posterior Probability (CLPP) was 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) 53 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 54 .Genotypes were imputed by the UKBB using the Haplotype Reference Consortium reference panel 55 with imputation  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.For a pair of related individuals, one was randomly excluded.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).Linear regression model (1) in the main text was performed.Covariates included age at assessment (21003-0.0),age 2 , sex (31-0.0), the first 10 PCs (22009-0.1 to 22009-0.10), the environmental 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  ×  analysis in 17 , 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 lipidlowering medication (6177-0.0)present.LDL-C measurements therefore did not consider medication if there were missing values for medication history.This introduced missing values in LDL-C for 248,419 individuals.
Since only the variants without either  ×  interaction or mediation are valid in the MR analysis, we assume  = 0 (no mediation) and  3,j = 0 (no interaction).We have  ̂ =  ̂1,j +  1  ̂3,j By applying Slutsky's theorem, and let  3,j = 0, we have: , which converges to 1 when  1 and  2 → ∞.However, when   2 is small (weak instrument in MR analysis), the converge of ( ̂) to 1 is slow.We also note that ( ̂) ≤ 1.
We propose  _ ,  _ can be performed by the IMRP method in MR analysis 24 .
Our procedure to test for the  ×  interaction is: 1) We apply  _ to search variants with joint effect of mediation and interaction effect; 2) we apply   for the variants detected by  _ .To reduce the effect by the mediation, we can exclude the variants associated with E.
The main effect size of the th variant was generate by  1 ∼ (0,   2 ), where   2 is the trait variance accounted for by the IVs.For the first variant, we added its interaction effect with E. The phenotype   by generated by where   ∼ (0,  2 ).The causal effect  was estimated using the last 100 variants as the IVs.The power and type I error rate for   and    were calculated based on the first and second variants, respectively.

Fig. 1 .
Fig. 1.Illumination of Mendelian Randomization and  × . 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.Right panel: the horizontal pleiotropic variants depart from the regression line and can be detected.B. Left panel: the  ×  framework, with the goal of testing  × .Instead of an explicit exposure, we create a pseudo exposure  ̃, which can be viewed as the total contribution of genetic variants to trait Y.The genetic variants associated with the pseudo exposure  ̃ but not through either the environment E or  ×  are valid instrumental

10
Methods.E. Type I error comparison for   ,  _ and two-step precedure.F. Power comparison for   ,  _ and two-step precedure.

Fig. 3 .
Fig. 3. Manhattan plots for interactions, marginal and main effect size comparisons.A-C: The circle Manhattan plots of gene × alcohol drinking interactions for LDL-C, HDL-C, and TG.The genome wide significant loci are presented in red dots.D-F: The marginal and main effect size comparisons for LDL-C, HDL-C, and TG.The colored circles represent the genome-wide

Fig. 4 .
Fig. 4. The estimated heritability of Cross-Population (A) and European population (B) using the LDSC regression.Marginal effect heritability refers to the heritability estimates through the marginal effect  ̂ (red bars), and interaction effect heritability refers to the heritability estimated through  ̂−  ̂ ̂1 (blue bars).Only cross-population and European combined effect of  ×  interaction and mediation: 0.1  + 0.05( 1 *  1 ) +   ,

Table 1 . Interaction loci screened by 𝑻 𝑴𝑹_𝑮𝒙𝑬 and followed by the direct test 𝑻 𝑫𝒊𝒓𝒆𝒄𝒕 in GLI (two-step test) and replicated by the direct test 𝑻 𝑫𝒊𝒓𝒆𝒄𝒕 in UK Biobank.
The causal effect  of the inverse variance weighted (IVW) is estimated by