TWB was approved by the Institutional Review Board on Biomedical Science Research/IRB-BM, Academia Sinica, and the Ethics and Governance Council of Taiwan Biobank, Taiwan. TWB approved my application to access the data on 18 February 2020 (application number: TWBR10810-07). The current work further received approval from the Research Ethics Committee of the National Taiwan University Hospital (NTUH-REC no. 201805050RINB).

Since October 2012, TWB has recruited Taiwan residents aged 30 to 70 and collected their genomic and lifestyle factors (Chen et al., 2016). After signing informed consent, community-based volunteers took physical examinations and provided blood and urine samples. TWB health professionals collected lifestyle information through a face-to-face interview with each participant.

As of March 2022, 27,719 and 103,332 individuals (aged 30–70 years) have been whole-genome genotyped by the TWB 1.0 and TWB 2.0 genotyping arrays, respectively. The TWB 1.0 array was designed for Taiwan’s Han Chinese, running on the Axiom Genome-Wide Array Plate System (Affymetrix, Santa Clara, CA). The TWB 2.0 array was developed according to the experience of designing the TWB 1.0 array and the next-generation sequencing of ∼1,000 TWB individuals. These two arrays were released in April 2013 and August 2018, respectively. Because the sample size of “individuals genotyped by the TWB 2.0 array” (called “the TWB2 cohort”) was larger than that genotyped by the TWB 1.0 array (called “the TWB1 cohort”), the TWB2 cohort was treated as a discovery set. In contrast, the TWB1 cohort was used as a replication set.

A total of 1,462 individuals were genotyped by both arrays. To ensure that the replication set was independent of the discovery set, I removed these 1,462 individuals from the TWB2 cohort. I also tried to exclude individuals with more than 10% missing in their genotype calls, where 10% is a commonly adopted cutoff in quality control (Band et al., 2019). Nonetheless, no individuals were removed due to this criterion.

To remove cryptic relatedness, I calculated PI-HAT = Probability (IBD = 2) + 0.5 × Probability (IBD = 1) by PLINK 1.9 (Purcell et al., 2007), where IBD is the genome-wide identity by descent sharing coefficients between any two individuals. I excluded one individual from each pair with PI-HAT ≥ 0.2, a cutoff value chosen by some studies (Calabro et al., 2015; An et al., 2019; Lin et al., 2020a; Lin, 2022a; Lin, 2022b). After this step, the TWB2 (discovery) and TWB1 (replication) cohorts included 86,536 and 25,460 individuals, respectively.

TWB 2.0 and TWB 1.0 arrays comprised 648,611 and 632,172 autosomal SNPs, respectively. I excluded 17,419 SNPs with Hardy-Weinberg test p-values <5.7 × 10⁻⁷ (WTCCC, 2007) and 22,614 SNPs with genotyping rates <95% from the TWB2 cohort and removed 6,900 SNPs with Hardy-Weinberg test p values <5.7 × 10⁻⁷ (WTCCC, 2007) and 27,628 SNPs with genotyping rates <95% from the TWB1 cohort. The remaining 608,578 TWB2 SNPs and 597,644 TWB1 SNPs were used to construct ancestry principal components (PCs).

The Michigan Imputation Server (https://imputationserver.sph.umich.edu/index.html) was further used to impute genotypes. The East Asian population from the 1,000 Genomes Phase 3 v5 was chosen as the reference panel. I removed SNPs with low imputation information scores (R-square <0.8), with imputation rates <95%, or with Hardy-Weinberg test p-values <5.7 × 10⁻⁷ (WTCCC, 2007). The TWB2 and TWB1 individuals were finally genotyped (or imputed) on 6,546,183 and 7,433,014 autosomal SNPs, respectively.

With a larger sample size, TWB2 (n = 86,536) was treated as the discovery cohort. A total of 4,807,430 TWB2 SNPs with minor allele frequencies (MAFs) ≥ 5% were analyzed sequentially. I skipped the epistasis analysis for SNPs with MAFs <5% due to the inferior genotyping (or imputation) accuracy for low-frequency variants (Mitt et al., 2017). Moreover, GxG studies usually focus on common SNPs because of their better reproducibility (Wang et al., 2019). If the sample size of any genotype combination (from an SNP pair) is too small, this GxG evidence can hardly be replicated in another cohort. Significant SNPs would be further investigated using the TWB1 (n = 25,460) again, which was regarded as the replication cohort. Totally 29 traits in eight categories were investigated, including

(A) Six lung function traits: vital capacity, tidal volume, inspiratory reserve volume, expiratory reserve volume, forced vital capacity, and forced expiratory volume in 1 s.

The total number of pairwise interaction tests among 4,807,430 SNPs is very huge. Considering 29 phenotypes, conventionally, investigators need to perform 4 , 807 , 430 2 × 29 = 3.35 × 10 14 GxG tests. The power to detect GxG will be largely reduced due to the harsh penalty of multiple testing. In the following, I introduce a vQTL method to prioritize SNPs exhibiting GxG. Let G ₁ and G ₂ be the numbers of minor alleles at two SNPs. Without loss of generality, modeling a phenotype (denoted as “Y”) with the two SNPs (G ₁ and G ₂ ) can be expressed as follows, Y = β 0 + β G 1 G 1 + β G 2 G 2 + β INT G 1 × G 2 + ε , (1) where ε is the random error term. Conditional on the genotype of the first SNP, the variance of Y can be derived as V a r Y G 1 = g 1 = V a r β 0 + β G 1 G 1 + β G 2 G 2 + β INT G 1 × G 2 + ε G 1 = g 1 = V a r β 0 + V a r β G 1 G 1 | G 1 = g 1 + V a r β G 2 G 2 + V a r β INT G 1 × G 2 | G 1 = g 1 + 2 C o v β G 2 G 2 , β INT G 1 × G 2 | G 1 = g 1 + V a r ε = 0 + 0 + β G 2 2 V a r G 2 + β INT 2 g 1 2 V a r G 2 + 2 β G 2 β INT g 1 V a r G 2 + V a r ε = β G 2 + β INT g 1 2 V a r G 2 + V a r ε (2)

I obtained V a r Y G 1 = g 1 = β G 2 + β INT g 1 2 V a r G 2 + V a r ε . In the absence of GxG, β INT = 0 and V a r Y G 1 = g 1 = β G 2 2 V a r G 2 + V a r ε , representing the variance of Y remains constant across the three genotypes of the first SNP ( g 1 = 0, 1, and 2). Therefore, I investigated GxG by testing equal variance (homoscedasticity) of a phenotype across the three genotype groups of each SNP.

To provide results robust to outliers and the distributions of traits, I performed the “rank-based inverse normal transformation” (RINT) transformation (McCaw et al., 2020) on each trait before the vQTL analysis. RINT-trait was transformed to be normally distributed through this step. Subsequently, I obtained “genotypes-and-covariates adjusted RINT-trait” via the residuals of regressing RINT-trait on genotypes and covariates. Specifically, for each SNP, I adjusted RINT-trait with genotype effects as two dummy variables and covariates, including sex (male vs. female), age (in years), body mass index (BMI, in kg/m²), performing regular exercise (yes vs. no), educational attainment (an integer from 1 to 7), smoking status (yes vs. no), drinking status (yes vs. no), and the first 10 ancestry PCs.

The abovementioned covariates are commonly adjusted for continuous phenotypes, because each can influence the phenotypes to some extent (Lin et al., 2020b; Lin, 2022c). Current smoking indicated “having smoked cigarettes for at least 6 months and having not quit smoking when joining the TWB. Drinking was defined as “having a weekly intake of more than 150 mL of alcoholic beverages for at least 6 months and having not stopped drinking when joining the TWB.” Regular exercise was defined as “performing exercise lasting for 30 min thrice a week.” Educational attainment was an integer ranging from 1 to 7: 1 (illiterate), 2 (no formal education but literate), 3 (primary school graduate), 4 (junior high school graduate), 5 (senior high school graduate), 6 (college graduate), and 7 (Master’s or higher degree). When analyzing the five obesity traits (BMI, BFP, WC, HC, and WHR), BMI was excluded from the covariates.

Through the above step, I obtained the “genotypes-and-covariates adjusted RINT-trait,” denoted as “ e i ” for the ith individual (i = 1, 2, … , n). The dispersion of e i was then calculated by D i = e i − e ∼ 2 , where e ∼ was the sample median of e i across all n individuals. Because of the robustness, the sample median is adopted instead of the sample mean. Subsequently, I regressed D i on the two dummy variables for genotype coding to check whether the dispersion of e i was dependent on different genotypes. The significance of the F-statistic of this regression model meant that the dispersion of “genotypes-and-covariates adjusted RINT-trait” ( e i ) varied with different genotypes, which was a clue of GxG according to the derivation of Equation 2. This is called the “scale test” for detecting trait heteroscedasticity across genotypes.

This scale test for vQTL identification is not novel. Several studies have previously introduced it (Pare et al., 2010; Soave et al., 2015; Soave and Sun, 2017; Staley et al., 2022; Westerman et al., 2022). Some investigators used Levene’s statistics (Levene, 1960) to test for equal variance (Pare et al., 2010; Westerman et al., 2022), while others reformulated the concept to a regression framework (Soave et al., 2015; Soave and Sun, 2017; Staley et al., 2022). However, these two approaches are conceptually identical. The workflow shown here is based on the regression framework with the potential to allow for continuous exposures (Soave et al., 2015; Soave and Sun, 2017; Staley et al., 2022). The R code to implement this regression for vQTL testing can be downloaded from https://github.com/WanYuLin/Univariate-scale-test-UST.

I analyzed 4,807,430 SNPs with MAFs ≥ 5% in TWB2. p-values of the scale test <0.05/(4,807,430 × 29) = 3.6E-10 were considered significant. Significant SNPs identified from the TWB2 cohort were further analyzed using the TWB1 cohort. SNPs with p-values of the scale test <0.05/(the number of SNPs analyzed in TWB1) were considered successfully replicated. These SNPs were called the vQTLs.

Based on the replicated vQTLs, I then performed an over-representation analysis on the Reactome pathway webpage (Fabregat et al., 2017) (https://reactome.org/). Over-representation analysis assesses whether some Reactome pathways are enriched (over-represented) in the vQTLs identified within the same trait category. It answers the question, “Do the vQTLs from a trait category contain more proteins for a certain pathway than that is expected by chance?” This over-representation test generates a p-value based on the hypergeometric distribution, which is then corrected by the Benjamini–Hochberg false discovery rate (FDR) procedure (Benjamini and Hochberg, 1995) (https://reactome.org/userguide/analysis).

Although the scale test is not novel, I still conducted a simulation study to assess its power of detecting GxG, given the sample sizes of TWB2 and TWB1. Without loss of generality, four SNPs (rs150856817 on chromosome 2, rs2075291 and rs662799 on chromosome 11, and rs4144003 on chromosome 16) were, in turn (one by one), regarded as SNP2. The TWB2 and TWB1 SNPs on chromosome 1 were, in turn, treated as SNP1. As derived by Equation 2, V a r Y G 1 = g 1 = β G 2 + β INT g 1 2 V a r G 2 + V a r ε . I tested whether the phenotypic variance depended on the three genotypes of SNP1. The trait value was then generated as follows, Y i = β G 1 G 1 i + β G 2 G 2 i + β INT G 1 i × G 2 i + ε i , i = 1 , ⋯ , 86536   or  25460 , (3) where G 1 i , G 2 i = 0 , 1, 2 represented the numbers of minor alleles at SNP1 and SNP2, respectively. The random error term ε i was generated from a standard normal distribution with a mean of 0 and a variance of 1. I considered two situations, as follows,

(I) Without SNP main effects: By specifying β G 1 = β G 2 = β INT = 0 (scenario 1), I evaluated the type I error rates given the significance level of 3.6E-10. To demonstrate that this method is also capable of detecting “pure epistasis” (no individual SNP main effects) (Russ et al., 2022), power was assessed at β G 1 = β G 2 = 0 and β INT = 0.3 . Because the four SNPs were, in turn, regarded as SNP2, the power was evaluated under scenarios 2–5. A total of 10⁸ replications were conducted under this situation, and each scenario was investigated with 2 × 10 7 simulation replicates.

To sum up, a total of 4 × 10 8 simulation replicates were performed. Each of the following four situations was evaluated with 10 8 replications: TWB2 without SNP main effects; TWB2 with SNP main effects; TWB1 without SNP main effects; and TWB1 with SNP main effects. In each situation, scenario 1 (in the absence of GxG) and scenarios 2–5 (in the presence of GxG) were simulated to assess the type I error rates and power, respectively. Each scenario (under each situation) was investigated with 2 × 10 7 replicates.

For any two vQTL SNPs, RINT-trait was regressed on the numbers of minor alleles (0, 1, or 2) of the two SNPs and their product term (interaction term) while adjusting for the same covariates: sex (male vs. female), age (in years), BMI (in kg/m²), performing regular exercise (yes vs. no), educational attainment (an integer from 1 to 7), smoking status (yes vs. no), drinking status (yes vs. no), and the first 10 ancestry PCs. Similarly, when analyzing the five obesity traits (BMI, BFP, WC, HC, and WHR), BMI was excluded from the covariates. To avoid the multicollinearity problem, I calculated the variance inflation factor (VIF) of each regression model using the R package “car” (Fox and Weisberg, 2018). A VIF <5 was considered acceptable, which was commonly used as the VIF cutoff value (Rogerson, 2001).

Figures 1, 2 show the quantile-quantile (QQ) plots and power curves of the vQTL method, where the power was evaluated at α = 3.6E-10. With 4 × 2 × 10 7 simulation replicates for scenario 1 under the four situations (TWB2 without SNP main effects; TWB2 with SNP main effects; TWB1 without SNP main effects; and TWB1 with SNP main effects), no type I errors were observed given α = 3.6E-10. The QQ plots (Figures 1, 2A) showed that the observed p-values matched the expected p-values under the null hypothesis (no epistasis). The type I error rates were well protected given no p-value <3.6E-10.

As shown by Figures 1, 2, the presence of SNP main effects boosted the statistical power (D compared with C) while preserving the validity of the vQTL method (B). As expected, a larger MAF at SNP1 or SNP2 also increased the power. Moreover, the vQTL method is also capable of detecting “pure epistasis” (without individual SNP main effects) (Russ et al., 2022), as demonstrated by Figures 1, 2C.

I identified 1,668 vQTL SNPs from all 29 TWB2 phenotypes (p < 3.6E-10). I then performed the scale test on these 1,668 SNPs using the TWB1 cohort. If the p-value <0.05/1,668 = 3.0E-5, the SNP was considered to be successfully replicated in the TWB1. I further used the PLINK clumping command (Purcell et al., 2007) to find 18 nearly independent vQTL SNPs with linkage disequilibrium measure r ² < 0.01 (Table 1). Only 9 out of the 29 phenotypes demonstrated evidence of vQTLs, including three lipid traits: LDL (3 vQTLs), TCHO (1 vQTL), and TG (1 vQTL); three blood traits: RBC (4 vQTLs), WBC (1 vQTL), and HB (1 vQTL); a kidney trait: UA (1 vQTL); a liver trait: TB (3 vQTLs); and a diabetes trait: FG (3 vQTLs). Except for APOA5-rs2075291, all vQTLs were also QTLs (the last column of Table 1).

Five vQTLs identified from lipid traits are located in the APOB, APOA5, APOC1, and TOMM40 genes (Table 1). Many of these vQTL SNPs have been reported to be associated with complex diseases or to exhibit interactions with other SNPs. For example, APOA5-rs662799 presented a solid association with TG levels in various ethnicity samples (Liu et al., 2012), and this SNP was identified as a TG-vQTL and a TG-QTL through my analysis (Table 1). Korean data showed that the two vQTL SNPs in the APOA5 gene, rs2075291 and rs662799 (Table 1), were associated with increased arterial stiffness and decreased adiponectin levels (Kim et al., 2018). 13 Alzheimer’s disease (AD) GWAS cohorts demonstrated that the vQTL SNP in the APOC1 gene, rs438811, significantly interacted with the APOE-ε4 allele. Carrying one minor allele T of rs438811 increased the AD risk by 26.75% in APOE-ε4 carriers (Zhang et al., 2018). Data recruited from China’s Second Affiliated Hospital of Xi’an Jiaotong University showed that AXIN1-rs2301522 was significantly associated with the risk of osteoporosis (Cui et al., 2022).

While vQTLs were searched through the scale test, QTLs were identified by the location test (testing whether the phenotype mean was dependent on the genotypes) (Staley et al., 2022). Some previous searches for genome-wide epistasis were based on QTLs (Bocianowski, 2013; Laurie et al., 2014; Yang et al., 2018). As a comparison, I also searched for epistasis through the QTL approach as follows. Specifically, I regressed RINT-trait on the number of minor alleles (0, 1, or 2) for each SNP respectively (one SNP at a time) while adjusting for the same 17 covariates, including sex (male vs. female), age (in years), BMI (in kg/m²), performing regular exercise (yes vs. no), educational attainment (an integer from 1 to 7), smoking status (yes vs. no), drinking status (yes vs. no), and the first 10 ancestry PCs. The R code for QTL testing can also be downloaded from https://github.com/WanYuLin/Univariate-scale-test-UST.

I found 48,484 QTL SNPs (p < 0.05/(4,807,430 × 29) = 3.6E-10) from the 4,807,430 TWB2 common variants (MAF ≥ 5%). I then performed the location test on these 48,484 SNPs using the TWB1 cohort. If the p-value <0.05/48,484 = 1.0E-6, the SNP was successfully replicated in the TWB1. I further used the PLINK clumping command (Purcell et al., 2007) to find 281 nearly independent QTLs with linkage disequilibrium measure r ² < 0.01. Figure 3 compares the vQTL and QTL approaches to search for epistasis. When resorting to the QTL approach, I also found that the SNPs exhibiting GxG were the vQTL SNPs in Table 1 (or located in the same genes as the vQTL SNPs). Hence, the vQTL approach is more efficient as only 18 2 = 153 GxG tests were performed for each phenotype, compared with 281 2 = 39340 SNP pairs analyzed for the QTL method.

According to the 18 vQTLs listed in Table 1, ( 18 2 ) = 153 GxG tests were performed for each phenotype, and 4,437 (=153 × 29) GxG tests were implemented for 29 phenotypes. Among the 4,437 tests, 15 were significant at p < 0.05/4437 = 1.1E-5 in the discovery (TWB2) cohort (Table 2). These 15 GxG interactions were further analyzed using the replication (TWB1) cohort, and the GxG interaction p-values were all less than 0.02 (Table 2). When resorting to the false discovery rate control (FDR) procedure (Benjamini and Hochberg, 1995), the FDRs were all less than 2%. Controlling FDR at 5% is commonly adopted in genomic studies (Sun et al., 2006). Therefore, an FDR level of 2% was acceptable, and all 15 GxG were replicated using the TWB1 cohort. Furthermore, the directions of GxG interaction coefficients (synergistic or antagonistic interaction effects) were consistent across the two cohorts. Because the 18 vQTL SNPs were nearly independent with linkage disequilibrium measure r ² < 0.01, putting any two and their interaction (product) term into the model did not cause multicollinearity problems. The VIFs of all regression models were controlled under 2 (Table 2).

Figures 4–6 show the GxG interaction plots for TCHO (1 GxG), FG (1 GxG), TB (2 GxG), RBC (6 GxG), HB (3 GxG), and HCT (2 GxG). The y-axis represents the averages of RINT-trait of nine genotype combinations of two SNPs. Lines with different slopes suggest that the effect of an SNP depends on the genotype of another SNP, which is a clue of GxG. Nonetheless, these interaction plots may not completely correspond to the GxG interaction p-values (Table 2). Unlike the epistasis test results (Table 2), these plots are descriptive summaries without adjusting for any covariate. For example, lines converging at genotype GG of rs2301522 (Figure 4E) represented that individuals with rs2301522-GG had similar RBC, while individuals with rs2301522-AA (or rs2301522-AG) had divergent RBC depending on rs143660108’s genotypes. Lines showing crosscut (Figure 5E) indicated “cross-over interaction,” meaning that rs143660108’s genotypes with the larger mean HB switched over at rs2301522-GG.

All GxG were observed for gene pairs on the same chromosome, except for the APOA5 (chromosome 11)—TOMM40 (chromosome 19) interaction for TCHO. Both APOA5 and TOMM40 are involved in lipid metabolism (Hishida et al., 2014). The GxG analysis showed that minor alleles of APOA5-rs662799 and TOMM40-rs61679753 exhibited significant synergistic interaction on TCHO (Figure 4A).

Because all phenotypes were RINT-transformed before the analysis, the 15 GxG interaction coefficients (Table 2) could be directly compared. Both cohorts suggested that LUC7L-rs143660108 and DECR2-rs7197553 presented the most substantial interaction on RBC. The minor alleles of these two SNPs exhibited a notable synergistic interaction on RBC (Figure 4F).

Table 3 shows the Reactome pathway analysis results on genes identified within the same trait category. Many lipid-related pathways were enriched (FDR <0.05) in the vQTLs identified from lipid traits (APOB, APOA5, APOC1, and TOMM40). “Glucuronidation” and “Paracetamol ADME” (Absorption, Distribution, Metabolism, Excretion) pathways were over-represented in the vQTLs of liver traits (DGKD, UGT1A10, and MROH2A). Glucuronidation is a major metabolic reaction in the liver (Yang et al., 2017), while paracetamol is also extensively metabolized in this organ (Forrest et al., 1982). Pathways related to beta cells were enriched in diabetes’ vQTLs (SPC25 and GCK). Beta cells are critical to diabetes by producing insulin to control blood glucose levels (Cerf, 2013). These pathway analysis results showed that the vQTLs found in this study were highly relevant to the traits.