Major histocompatibility complex harbors widespread genotypic variability of non-additive risk of rheumatoid arthritis including epistasis

Genotypic variability based genome-wide association studies (vGWASs) can identify potentially interacting loci without prior knowledge of the interacting factors. We report a two-stage approach to make vGWAS applicable to diseases: firstly using a mixed model approach to partition dichotomous phenotypes into additive risk and non-additive environmental residuals on the liability scale and secondly using the Levene’s (Brown-Forsythe) test to assess equality of the residual variances across genotype groups per marker. We found widespread significant (P < 2.5e-05) vGWAS signals within the major histocompatibility complex (MHC) across all three study cohorts of rheumatoid arthritis. We further identified 10 epistatic interactions between the vGWAS signals independent of the MHC additive effects, each with a weak effect but jointly explained 1.9% of phenotypic variance. PTPN22 was also identified in the discovery cohort but replicated in only one independent cohort. Combining the three cohorts boosted power of vGWAS and additionally identified TYK2 and ANKRD55. Both PTPN22 and TYK2 had evidence of interactions reported elsewhere. We conclude that vGWAS can help discover interacting loci for complex diseases but require large samples to find additional signals.

Scientific RepoRts | 6:25014 | DOI: 10.1038/srep25014 thus give rise to a new avenue of mapping genetic variants 13,14 . Hence, vGWAS has a great advantage in prioritizing potentially interacting loci without requiring prior knowledge of interaction types and interacting factors 15,16 . Additional explicit tests of GxE and/or GxG interactions are needed but only for the identified vGWAS loci, leading to a power increase attributed to a much reduced number of multiple tests 17,18 . Furthermore, studying biology of the identified vGWAS loci may inform what environmental factors should be further investigated. Therefore, vGWAS could be a cost effective approach to boost studies of GxG and GxE interactions based on existing GWAS data. However, vGWAS loci may not necessarily interact because other factors such as overdominance and scaling (i.e. various transformations) could also generate apparent genotypic variability [19][20][21] . Therefore, explicit tests of interactions and careful interpretation are essential to claim interacting signals.
There have been several successful applications of vGWAS in quantitative traits using either GWAS 18,19,22 or gene expression data 17,23,24 . For example, Brown et al. 17 reported abundant vGWAS signals within the major histocompatibility complex (MHC) region using RNA sequence data from lymphoblastoid cell lines derived from the TwinsUK cohort and subsequently confirmed GxE and GxG interactions underlying the signals. Yang et al. 19 identified a well-known locus FTO in a vGWAS of body mass index in 170,000 samples from multiple GWAS cohorts, which is believed to interact with physical activity and/or lifestyle 25 . Intriguingly, the majority of the reported vGWAS loci are known major GWAS loci (i.e. with strong additive effects based on allelic means): FTO in body mass index 19 , LEPR in C-reactive protein levels and ICAM1 in soluble ICAM1 levels 18 , SLC2A9 in serum urate 22 , MHC in autoimmune diseases 5 . These observations led to a concern whether vGWAS might simply rediscover GWAS signals due to the mean-variance correlation, i.e. a distinct difference in allelic means likely accompanied by obvious genotypic variability in a locus [19][20][21] . Nevertheless, increasing evidence suggest that these major GWAS loci tend to have multiple roles and are indeed involved in interactions [26][27][28][29] , therefore it is biologically plausible to see them having signals in both GWAS and vGWAS analyses.
So far vGWAS has been little exploited in complex diseases probably because phenotypes classified either as a disease case or a healthy control provide little continuous variation. This problem can be remedied by applying a two-stage vGWAS approach: firstly using a mixed model method implemented in software such as GCTA 30 to partition polygenic additive risk and non-additive environmental residuals for each sample on the liability scale [31][32][33] and secondly using the Levene's (Brown-Forsythe) test to assess equality of the residual variances across genotype groups per marker. It is therefore of interest to find out if vGWAS can also identify interacting loci in complex diseases. Here we report a pilot vGWAS of RA using a large GWAS cohort for discovery and two independent cohorts for replication. We aim to demonstrate that vGWAS is also useful to study disease traits.

Materials and Methods
Study cohorts and quality control (QC). We used the Rheumatoid Arthritis Consortium International for Immunochip samples recruited in the UK (RACI-UK), US (RACI-US) and Sweden (RACI-SE, i.e. Swedish Epidemiological Investigation of Rheumatoid Arthritis) as described previously 2 . Each cohort was genotyped with the high density Immunochip 34 in accordance with Illumina protocols and has been described in detail elsewhere 2 . Briefly, all RA patients fulfilled the 1987 criteria of the American College of Rheumatology and were tested for anti-citrullinated peptide antibody (ACPA) with a status of either seropositive (ACPA + ), or seronegative (ACPA -) or unassigned. All participants provided written informed consent for participation. Further details of the sample collections in the three study cohorts can be found in the Supplementary Note by Eyre et al. 2 . This study was approved by the North West Research Ethics Committee (MREC 99/8/84). All experiments (e.g. genotyping) were performed in accordance with relevant guidelines and regulations.
We excluded SNPs on sex chromosomes and samples of non-European origin. RA cases that were either ACPA-or ACPA unassigned were also excluded in each cohort to avoid disease heterogeneity. A rigorous QC was conducted using PLINK 35 for each cohort as described previously 2 followed by additional QC based on the criteria suggested for accurate GCTA prediction of genetic risk 31,32 : minor allele frequency > 0.01, SNP call rate > 0.95, sample call rate > 0.99, deviation from Hardy-Weinberg Equilibrium P < 1.0e-04. Figure 1 illustrates the workflow of the proposed two-stage vGWAS approach with an extension of explicit interaction tests. In this study GCTA was used to compute the genetic relationship matrix (GRM) and subsequently the first ten principal components (PCs) and then polygenic liability risk for each unrelated individual was predicted by imposing a GRM relatedness threshold of 0.15 recommended for Immunochip 33 while setting the disease prevalence as 0.01 and fitting gender and the first 10 PCs as covariates in the mixed model. The resultant environmental residuals were tested for variance heterogeneity across three SNP genotypes using an R package VariABEL 36 implemented the Levene's (Brown-Forsythe) test that requires no assumption of normally distributed phenotypes:

Two-stage vGWAS and Statistical analysis.
where the environmental residuals were the trait y; Z i = |y i -ỹ gi | is the deviation of y of the i th sample (y i ) and the median of y in samples with genotype g (ỹ gi ); N is the sample size and k is the total possible genotypes; n j is the number of samples with genotype j; Z j . is the mean deviation from the median for genotype j and Z. is the mean deviation from the overall median. When N is large, T 2 is an approximate χ 2 test taking two degrees of freedom.
To be extra cautious about the impact of normality, we also compared the vGWAS results based on the raw environmental residuals with those based on fully normalized environmental residuals using rank transformation available in the GenABEL package 37 . We found the normalization led to overly conservative vGWAS results in RACI-US and RACI-SE (Supplementary Figure S1) and thus reported only the results based on the raw environmental residuals.
Scientific RepoRts | 6:25014 | DOI: 10.1038/srep25014 We derived the significance threshold as 2.5e-05 based on permutation with 1000 iterations where a vGWAS was carried out using the randomly permutated environmental residuals and the lowest P value was recorded in each iteration (Supplementary Figure S2). We performed vGWAS in each of the three cohorts and used RACI-UK as the discovery cohort and examined any vGWAS significant SNPs for direct replication in RACI-US and RACI-SE, i.e. the same SNP with a vGWAS P value less than 5.0e-02.
For comparison we also performed a conventional GWAS for each cohort using PLINK and following the previous study 2 that used a subset of independent common SNPs (i.e. minor allele frequency > 0.05 and minimum inter-SNP linkage disequilibrium (LD)) to compute GRM and subsequently the first 10 PCs. A forward selection approach was applied to the genome-wide significant (P < 5.0e-08) GWAS SNPs identified in RACI-UK iteratively to identify a set of independent signals to represent the additive background. In each iteration the logistic regression model fitting the covariates and any pre-selected SNPs was used to test all remaining SNPs and select the most associated SNP for next round until none had a P value less than 5.0e-02.
We tested interactions between the identified vGWAS SNPs in RACI-UK as follows: a) PLINK was used to compute pairwise LD (in r 2 ) for all the identified SNPs and then only pairs of SNPs in low LD (i.e. r 2 < 0.1) 9 were selected for interaction tests; b) the logistic regression model was used to test interactions for each of the k selected pairs of SNPs by fitting a pair of SNPs and their interaction terms as well as covariates of gender and the first 10 PCs derived from the GRM based on the subset of independent SNPs, with a significance of threshold defined as 0.05/(k * (k-1)/2); c) significant pairwise interactions were assessed while conditioning on the additive background of independent GWAS signals and unique pairs of interactions were identified using the forward selection approach above. Variance explained by the selected epistatic SNP pairs was calculated using a full logistic regression model fitting all the covariates, the additive background and the independent pairs of SNPs and their interactions. The identified epistatic SNP pairs were tested for replication individually in RACI-US and RACI-SE and the total variance explained was calculated following the approach above.
We further merged the three study cohorts into one and then performing the two-stage vGWAS as well as GWAS adjusted for cohort and the covariates described above to explore the impact of increased sample size. After the QC and GCTA analysis the combined cohort had 18,405 unrelated samples (5954 cases and 12,451 controls) and 107,144 SNPs in the vGWAS and GWAS. We quoted SNP genomic locations in the GRCh38/hg38 version throughout.

Results
The vGWAS analysis identified 1639 significant SNPs (P < 2.5e-05) in the RACI-UK cohort where the Immunochip-wide heritability of ACPA + RA was estimated as 24.7% (Table 1). About 92% of the identified vGWAS SNPs mapped to MHC and the remaining mapped to PTPN22 -both were major RA associated loci in the conventional GWAS that in addition identified TNFAIP3 and TYK2 at the genome-wide significance level (P < 5.0e-08) (Fig. 2). Forward selection of the significant GWAS SNPs led to an additive background composed of 20 independent SNPs (Supplementary Table S1). The quantile-quantile (QQ) plots of both the vGWAS  (Supplementary Figure S1) and GWAS (Supplementary Figure S3) analyses in RACI-UK showed no sign of inflation.
In contrast, the vGWAS analyses detected only 154 SNPs in RACI-US (heritability estimated as 24.6%) and 343 in RACI-SE (heritability estimated as 18.5%) under the same significance threshold, suggesting the two independent cohorts were underpowered probably due to the relatively small samples used ( Table 1). The QQ plots for RACI-US and RACI-SE indicated slight deflation in the vGWAS analyses (Supplementary Figure S1) but not in the GWAS analyses (Supplementary Figure S3). Nearly all the vGWAS SNPs identified in both RACI-US and RACI-SE also mapped to MHC, except for rs2322659 (P = 7.7e-06) mapped to LCT in RACI-US (Supplementary Figures S4 and S5).
We found that a quarter (i.e. 413 out of 1639) of the identified vGWAS SNPs in RACI-UK had direct replication in both the RACI-US and RACI-SE cohorts (Supplementary Table S2) and they were all within the MHC region spanning the Class III and Class II. In addition, the two top vGWAS SNPs rs2476601 (P = 6.7e-22) and rs6679677 (P = 1.4e-21) within PTPN22 in RACI-UK each had a direct replication only in RACI-SE with a P value of 3.3e-02 and 3.0e-02 respectively, suggesting genotypic variability within the locus was detectable only with a large sample size. The top vGWAS SNP rs2322659 within LCT in RACI-US had no replication in either RACI-UK or RACI-SE.
Interestingly, in RACI-UK the most strongly associated vGWAS SNP rs1964995 (intergenic between HLA-DRA and HLA-DRB5) was about 134 kb away from the most strongly associated GWAS SNP rs6931277 (near HLA-DRB1) (Fig. 3a), and the two SNPs were also the top vGWAS and GWAS signals respectively in both RACI-US and RACI-SE (Supplementary Table S2). These results (re)confirm that the MHC region harbors both strong additive and non-additive effects 29,38 and suggest that the top non-additive sites could differ from the top additive loci. Nevertheless, the top GWAS SNP rs6931277 could also be involved in interactions as it had a strong vGWAS signal (i.e. P = 1.2e-96, 2.9e-07 and 7.1e-11 in RACI-UK, RACI-US and RACI-SE respectively).
We tested interactions for 11,871 pairs of the RACI-UK significant vGWAS SNPs with r 2 < 0.1. Based on a significance threshold of 4.2e-06 derived from Bonferroni correction, we found a large number of significant interactions all within MHC, of which 2,962 pairs of SNPs had a conditional interaction P value between 3.4e-07 and 5.0e-02 independent of the additive background. Applying forward selection to the 2,962 pairs and requiring a conditional interaction P < 1.0e-03 in each step, we identified 10 independent epistatic pairs of SNPs that spanned the entire MHC region and were replicated in both RACI-US and RACI-SE except for the rs9366778 -rs2239707 pair in RACI-US (Table 2, Supplementary Figure S6). In contrast to the additive background that explained 17.8% of the phenotypic variance in RACI-UK, the 10 epistatic pairs jointly explained about 2.5% of the phenotypic variance mainly (i.e. 1.9%) by interactions, suggesting the existence of widespread GxG interactions within MHC that are individually moderate/weak (on average explaining <0.2% of phenotypic variance) but jointly appreciable. The additive background and the 10 epistatic pairs explained 17.1% and 1.7% (0.8% by interactions) of the phenotypic variance in RACI-US, and 14.5% and 1.7% (0.7% by interactions) in RACI-SE.
Combining the three study cohorts together boosted power for both the GWAS and vGWAS analyses as expected (Supplementary Figure S7), with no signs of inflation or deflation (Supplementary Figure S8). The vGWAS analysis identified two new significant loci ANKRD55 and TYK2 in addition to MHC and PTPN22, which were all genome-wide significant in the GWAS analysis. The top vGWAS SNPs rs71624119 for ANKRD55 (P = 4.9e-07) and rs34536443 for TYK2 (P = 7.0e-06) both had only weak or no signals in individual member     cohorts, i.e. the P values were 4.4e-03, 2.5e-02 and 7.9e-02 in RACI-UK, RACI-US and RACI-SE respectively for rs71624119; and 1.8e-04, 3.4e-01 and 2.1e-01 respectively for rs34536443. In contrast, the GWAS analysis identified additional 12 genome-wide significant loci that were insignificant in the vGWAS (Supplementary Figure S7). We plotted the vGWAS P values against the GWAS counterparts for the Combined cohort as well as each individual cohort to assess the impact of the mean-variance correlation (Fig. 3). In the Combined cohort, high concordance between GWAS and vGWAS was observed when the MHC and/or PTPN22 regions were included (Fig. 3a,b) but greatly diminished if excluding both MHC and PTPN22 (Fig. 3c). In individual cohorts, GWAS and vGWAS concordance was observed only in RACI-UK even when including PTPN22 (Fig. 3d,e,f). These results suggested that the mean-variance correlation might influence vGWAS directly mainly at the major GWAS loci. In addition, we reran the vGWAS analyses while fitting the additive effects of the pre-identified GWAS SNPs as additional covariates at either stage in RACI-UK and found little differences in the results, indicating the vGWAS signals were unlikely driven by the differences of allelic means or scaling factors.

Discussion
We showed that the two-stage vGWAS approach successfully identified MHC that harbored widespread genotypic variability of non-additive risk of RA across the three study cohorts. Such non-additive genotypic variability could be explained in part by GxG interactions across the entire MHC region that were independent of the strong additive effects, and had a moderate/weak effect individually but an appreciable joint effect. We also identified PTPN22 as a significant vGWAS locus in the discovery cohort that however was replicated in only one independent cohort. Combining the three study cohorts boosted the power of vGWAS and detected two additional loci ANKRD55 and TYK2 that are yet to be replicated. These results suggest that the two-stage vGWAS indeed works for complex diseases.
Our results of multiple GxG interactions within MHC are in line with the recent reports that focused on a few preselected MHC loci (e.g. HLA-DRB1) and used high resolution imputation and a much increased sample size to fine map within-locus interactions 29,39 . Similar fine mapping efforts are needed to accurately compute the additive background and subsequently refine GxG interactions across the entire MHC region. There is not information of environmental factors other than gender in the study cohorts to test GxE interactions. However, there have been several reports of GxE interactions (e.g. with smoking and alcohol) in RA involving MHC and/or PTPN22 [40][41][42][43][44] . Besides, there is also evidence of non-additive effects and heterozygous advantage in MHC variants from other autoimmune diseases 29,38,39 , all supporting the notion that MHC is highly diverse owing to strong natural selection pressure in order to cope with ever changing environmental pathogens 45 .
One major concern about vGWAS applications is potential conflation caused by the mean-variance correlation that is hardly disentangled, i.e. variance can be mathematically modelled as some function of the mean 20,46 . Here we showed that the impact of such a correlation was strong only in the major RA associated loci (i.e. MHC and PTPN22) that tended to be involved in GxG and/or GxE interactions. We argue that co-existence of strong additive and non-additive effects in the major associated loci may be biologically important 9 rendering these loci the key regulators of the study phenotype and even other related phenotypes (e.g. MHC regulates a number of autoimmune diseases). From this perspective, detection of the major GWAS loci should not be regarded as a rediscovery but in fact an assurance that vGWAS can at least pick up the key interacting loci for dichotomous disease phenotypes.
vGWAS may be more fruitful in meta-analysis of multiple cohorts as demonstrated in the Combined cohort analysis. We showed that the newly discovered vGWAS loci, like those in GWAS meta-analysis, had weaker association signals than MHC or PTPN22, suggesting that they are unlikely involved in strong interactions and thus require more samples to detect 9 . ANKRD55 is associated with multiple diseases including RA and multiple sclerosis 3,47 but its biological functions and involvement in statistical interaction are not yet clear. TYK2 is also a pleiotropic gene associated with multiple autoimmune diseases and was found interacting with IRF5 in lupus but the interactions seemed weak and inconsistent [48][49][50] . Further work is required to validate their non-additive roles in RA. We anticipate that future enlarged vGWAS meta-analyses may detect additional, especially novel loci that may be driven more by the variance rather than the mean.
vGWAS is also useful to increase power of detection of interactions via an effective reduction of multiple tests concerning only the prioritized loci. For example, testing GxG interactions benefited from a much relaxed significance threshold as 4.2e-06 after filtering on vGWAS significant SNPs and then their pairwise LD in this study. Caution is recommended when using such filtering because vGWAS could miss true interacting loci with limited genotypic variability 9,11,12 . Nevertheless, such a filtering is practically useful to prioritize the main non-additive loci particularly when a big sample size is available.
This pilot study has a few limitations. Firstly, the low power of detection in individual cohorts (particularly in RACI-US and RACI-SE). In addition to sample size and SNP density, power of a vGWAS is also determined by the effect sizes of the interactions involved and the unobserved interacting factors 11,36 . The MHC signals appeared significantly in every cohort probably because they were involved in multiple interactions as suggested previously 17 . Secondly, uncertainties associated with the estimated environmental residuals due to hidden or unaccounted factors. Factors influencing the accuracy of polygenic risk prediction (e.g. sample size, case control ratio, SNP density, population structure, sample relatedness and ascertainment) 31-33 will also influence the vGWAS power. Population histories and/or special environments may also cause biases if not unaccounted for properly. For example, RACI-SE was the only cohort using matched cases and controls 2 , which might explain in part the relatively low estimate of heritability (Table 1) and the slight deflation observed in the vGWAS QQ plot (Supplementary Figure S1). Thirdly, this study concerned only the Immunochip platform for an autoimmune disease. Further investigations using GWAS arrays and other types of diseases would be useful to provide a complete view of this approach. A carefully designed simulation study is also recommended to assess the factors above as well as to quantify the impact of the mean-variance correlation comprehensively. Such a simulation study may also help develop a consensus protocol to facilitate future vGWAS meta-analysis of diseases.