The combined effects of FADS gene variation and dietary fats in obesity-related traits in a population from the far north of Sweden: the GLACIER Study.

BACKGROUND: Recent analyses in Greenlandic Inuit identified six genetic polymorphisms (rs74771917, rs3168072, rs12577276, rs7115739, rs174602 and rs174570) in the fatty acid desaturase gene cluster (FADS1-FADS2-FADS3) that are associated with multiple metabolic and anthropometric traits. Our objectives were to systematically assess whether dietary polyunsaturated fatty acid (PUFA) intake modifies the associations between genetic variants in the FADS gene cluster and cardiometabolic traits, and to functionally annotate top-ranking candidates to estimate their regulatory potential.
METHODS: Data analyses consisted of the following: interaction analyses between the 6 candidate genetic variants and dietary PUFA intake; gene-centric joint analyses to detect interaction signals in the FADS region; haplotype-centric joint tests across 30 haplotype blocks in the FADS region to refine interaction signals; and functional annotation of top-ranking loci from the previous steps. These analyses were undertaken in Swedish adults from the GLACIER Study (N = 5,160); data on genetic variation and eight cardiometabolic traits were used.
RESULTS: Interactions were observed between rs174570 and n-6 PUFA intake on fasting glucose (Pint = 0.005) and between rs174602 and n-3 PUFA intake on total cholesterol (Pint = 0.001). Gene-centric analyses demonstrated a statistically significant interaction effect for FADS and n-3 PUFA on triglycerides (Pint = 0.005) considering genetic main effects as random. Haplotype analyses revealed three blocks (Pint < 0.011) that could drive the interaction between FADS and n-3 PUFA on triglycerides; functional annotation of these regions showed that each block harbours a number of highly functional regulatory variants; FADS2 rs5792235 demonstrated the highest functionality score.
CONCLUSIONS: The association between FADS variants and triglycerides may be modified by PUFA intake. The intronic FADS2 rs5792235 variant is a potential causal variant in the region, having the highest regulatory potential. However, our results suggest that multiple haplotypes may harbour functional variants in a region, rather than a single causal variant.

Introduction

The human fatty acid desaturase gene (FADS) cluster on chromosome 11 harbours multiple genes including FADS1, FADS2 and FADS3 (1). FADS1 and FADS2 encode the Δ-5 and Δ-6 desaturases, respectively. These desaturases introduce cis double bonds at the 5th and 6th positions from the carboxyl end in the biosynthesis of long-chain polyunsaturated fatty acids (LC-PUFAs) from shorter chain polyunsaturated fatty acid (PUFA) precursors during fatty acid metabolism. Genetic variants in FADS1/FADS2 are known to correlate with the variability of LC-PUFAs levels in serum/plasma (2–5), erythrocyte membranes (6, 7) and adipocytes (8), and are associated with inflammation (6, 9), cardiovascular disease (10, 11) and type 2 diabetes (12, 13). Genome-wide association studies have established associations between FADS1/FADS2 genetic variants and blood lipid traits (e.g., low-density lipoprotein cholesterol (LDL-C), high-density lipoprotein cholesterol (HDL-C), triglycerides, and total cholesterol) (14) and fasting glucose levels (15). FADS3 shares high sequence homology with FADS1 (62%) and FADS2 (70%), and has similar intron and exon structures (1). However, the function of FADS3 is yet to be determined.

The FADS region has been the focus of numerous gene-diet interaction studies (16–23). Fumagalli et al. recently detected a strong positive selection signal in the FADS region in a Greenlandic Inuit cohort (24). The authors identified six single nucleotide polymorphisms (SNPs) associated with multiple quantitative traits, with the strongest association observed for body weight and height. These findings were validated in an additional Inuit cohort and the association with weight was confirmed in further European cohorts. The Greenlandic Inuit is a population isolate that has been residing in Greenland for more than 1,000 years, subsisting predominately on marine animals that are rich in monounsaturated and n-3 polyunsaturated fatty acids; these specific dietary constraints are thought to have resulted in enrichment for genetic variants involved in fatty acid desaturation. We hypothesize that the genetic adaptation that drove the observed changes in allele frequencies at certain FADS loci in the Inuits might be partially due to gene × PUFA interactions on traits, such as height, weight, and blood lipid and glucose concentrations, factors directly related to evolutionary fitness. Thus, it would be reasonable to speculate that the enhanced genetic effects of FADS variants on height, weight and other metabolic traits shown in Inuit may be driven by high dietary PUFA concentrations in this population. Interactions may be evident in other Arctic and sub-Arctic populations whose diets vary in fatty acid composition.

The primary aim of this study was to test whether dietary PUFA intake modifies the associations of the six candidate FADS SNPs reported by Fumagalli et al. and eight obesity-related traits in the GLACIER Study, a population-based cohort from the far north of Sweden. Our secondary aim was to conduct gene- and haplotype-based interaction tests focused on the FADS locus and functionally annotate top-ranking variants from relevant interactions using state of the art in silico tools.

Materials and Methods

Study participants

The GLACIER Study (cohort registration number: ISRCTN35275922) is a prospective population-based cohort study of 19,547 Swedish adults from a sub-Arctic population in northern Sweden; roughly half the background population have ancestral lineages within the region (Västerbotten), dating back to the 1700s (25). Clinical and lifestyle data in GLACIER participants were collected within the on-going Västerbottens Hälsoundersökning (Västerbotten’s Health Survey, VHU; also known as the Västerbotten Intervention Programme) where residents in the county are invited to their primary care centre to complete a detailed health and lifestyle survey and clinical examination at ages of 40, 50, and 60 years (26). The Regional Ethical Review Board in Umeå approved the study protocol and all study participants provided written informed consent as part of VHU. The study sample is a sub-cohort of 5,160 participants, who received health examinations between 1985 and 2004 and have complete genotype and phenotype data for the analyses.

Clinical measurements

The clinical protocol is described in detail elsewhere (27). Anthropometric data were obtained by trained nurses: body weight was measured to the nearest 0.1 kg and height was measured to the nearest 1 cm. Capillary blood was drawn following an overnight fast and a second sample was drawn 2 hours after a 75-gram oral glucose load. Blood glucose and serum lipid levels were measured with a Reflotron bench-top analyzer (Roche Diagnostics Scandinavia AB). Approximately 5% of the participants reported being fasted for less than 8-hours and 15% of the participants did not report fasting duration. LDL-C levels were calculated using the Friedewald formula (28). For people on lipid lowering medications, lipid levels were corrected by adding a constant for statins (+0.208 mmol/l for triglycerids, +1.347 mmol/l for total cholesterol, -0.060 mmol/l for HDL-C and +1.290 mmol/l for LDL-C) (29), as previously described (30). Data outside the following trait-specific intervals were treated as missing: BMI, 15-70 kg/m2; weight, ≥ 35 kg; height, 130-210 cm; total cholesterol, 0.5-15 mmol/l; triglycerids, 0.15-20 mmol/l; HDL-C, 0.15-7 mmol/l; fasting glucose 1-25 mmol/l; 2-h glucose, 1-35 mmol/l. Food intake level (total energy intake divided by estimated basal metabolic rate, 0.63-2.3) was used to exclude participants with improbable self-reported total energy intake (31). Participants with >10% genotype data missing for the whole genotyping array were excluded.

Dietary intake assessment

A validated self-administered food frequency questionnaire (FFQ), designed to quantify habitual diet during the past year, was administered in connection with the health examination (32, 33). The FFQ initially consisted of 84 foods, food groups and beverage items, but it was shortened to 66 items in 1996 by collapsing similar items and removing others. Daily n-3 (calculated as α-linolenic acid + eicosapentaenoic acid [EPA] + docosahexaenoic acid [DHA]), n-6 (calculated as linoleic acid + arachidonic acid) and total PUFA intake (calculated as n-3 + n-6 PUFA) and total energy intake were then calculated using the Swedish Food Composition Database (34).

Variant selection and genotyping

Genomic DNA was extracted from peripheral leukocytes and diluted to 4 ng/μl. Genotyping was performed using Cardio-MetaboChip array (35) at the Wellcome Trust Sanger Institute (Hinxton, UK). Genotypes for all available variants in the FADS1-FADS2-FADS3 region (Chr.11: 61317028-61416099, build 36; n=290) were extracted and, under the hypothesis that the regulatory regions of these genes may be proximal to the exonic sites, genotypes for variants within 200 kb upstream and downstream of the FADS region (n=436) were also extracted (total n=726). Rare variants with minor allele count <10 were excluded (n=284). Variants with a Hardy-Weinberg equilibrium P-value <0.0001 (n=3) were flagged but not excluded from the analyses, as potential deviations from expected genotype frequencies might reflect evolutionary processes, which are of special interest for the FADS region (24). Haplotype blocks were constructed based on the algorithm developed by Gabriel et al. (36) using Haploview Java 1.0 (37).

Statistical analyses

Statistical analyses were performed using R software v3.2.2 (38), PLINK v1.07, gPLINK v2.050 (39). Statistical analyses were carried out in four steps. These steps represent distinct analytic procedures that generally build on the previous steps. The study steps are described below and illustrated in the project flowchart (Figure 1):

Figure 1

Project flowchart.

Step 1: generalised linear regression models were fitted to assess the pairwise interaction effects of six candidate SNPs (rs174570, rs174602, rs74771917, rs3168072, rs12577276, rs7115739) (24) and dietary PUFA intake in relation to eight quantitative traits (weight, BMI, HDL-C, LDL-C, triglycerides, total cholesterol, fasting glucose and 2-h glucose);

Step 2: gene-centric joint analyses were performed to assess gene × dietary PUFA interaction effects in the whole FADS region, thereby targeting the quantitative traits and PUFA type to be studied in the next step;

Step 3: to study the statistically significant results from Step 2, haplotype-based gene × dietary PUFA interaction analyses were performed to narrow down haplotype blocks that are likely to interact with PUFA intake in the FADS region;

Step 4: functional annotation of candidate SNPs from Steps 1 and 3 was undertaken.

All models were adjusted for age, age2, sex, FFQ version, and the first four genomic principal components to control for population stratification. To account for potential confounding by total energy intake, we added dietary fatty acid residuals (regressed on total energy intake) along with total energy intake as covariates (40). The models with body weight as the outcome were adjusted for height. We added fasting status known (1/0 for yes vs. no) and fasting hours (1/0 for >8 h vs. 4–8 h) as dummy covariates in the interaction tests with blood lipid and glucose levels as outcomes. In models where 2-h glucose was the outcome, fasting glucose was included as a covariate. In Step 1, where the six candidate SNPs were analysed, P<0.05 was considered statistically significant as there is a strong prior for these tests, given previous published analyses (16–23).

For the gene/haplotype based analyses in Steps 2 and 3, we applied the method proposed by Chen et al. (41), which is an extension of the sequence kernel association test (SKAT) (42), but further incorporates a joint test of main and interaction effects to maximise the power to detect associations/interactions. For each haplotype block we calculated P values for the interactions with fixed genetic main effects (INT-FIX), interactions with random genetic main effects (INT-RAN), and joint main effects and interaction (JOINT), whilst adjusting for covariates. Bonferroni correction was applied to correct for multiple testing in gene-centric and haplotype-based analyses. As the dietary PUFA exposures and metabolic traits were not independent, we calculated the total number of effective environmental factors (1.37 instead of 3), and effective phenotypic traits (7.38 instead of 8) by accounting for the collective correlation of the two sets (43). In Step 2, P<0.005 [0.05/(7.38 × 1.37)] was considered statistically significant. In Step 3, P<0.002 [0.05/30] for testing 30 haplotype blocks was considered statistically significant (only one metabolic trait and outcome used here).

Functional annotation

In Step 4, the six candidate SNPs from Step 1 and candidate variants located in the haplotypes showing the strongest statistical evidence for interactions from Step 3 were functionally annotated. ChromHMM chromatin state predictions (44) were used to classify the variants into enhancers, repressors, promoters, insulators or others based on chromatin states in nine human cell lines using ANNOVAR (45). The web based database 3DSNP (46) was used to narrow down the candidate variant list and identify potentially causal variants based on functionality. Initially, we included the six candidate SNPs (rs174570, rs174602, rs74771917, rs3168072, rs12577276, rs7115739) and all common variants (build 137) located in the haplotype blocks generated from the interaction analyses in Step 3. We then screened for variants in high LD (r2>0.8) with the candidate SNPs in the European population (imputed to 1000 Genomes Phase 3 data, using 3DSNP (46)). If a variant was in high LD and had a higher functionality score than its original counterpart, then it was determined to be the leading variant and was prioritized over the original candidate variant in subsequent analyses. All variants were assigned functionality scores based on six functional parameters: i) evidence for disruption of transcription factor binding sites (TFBS); ii) evolutionary conservation, iii) ability to alter sequence motifs; iv) being located in a promoter region; v) being located in an enhancer region; and vi) number of topological interactions with distant genomic regions (genes and/or variants) via chromatin loops (46). The candidate variants and leading variants were then ranked based on the functionality score. The overall top-ranking leading variant and the most significant candidate SNP were visualized according to six functional parameters and topological features using radar charts and Circos plots, respectively. We searched for quantitative trait loci (QTL) related to lipid concentrations and expression QTL (eQTL) using HaploReg (47). Gene ontology enrichment analysis was performed on distal interacting genes using the PANTHER Overrepresentation Test to identify relevant biological processes and pathways (48).

Results

Characteristics

Characteristics for the 5,160 GLACIER participants (38.6 % males and 61.4 % females) are shown in Table 1. Allele frequencies of the six SNPs in the GLACIER cohort are shown in S1 Table. With the exception of rs174570 and rs174602, the derived alleles are of low frequency in this northern Swedish population (DAF<5%). Dietary n-3, n-6 and total PUFA intake were highly correlated with each other (r2 > 0.9 for all pairwise Pearson correlations).

Table 1

Population characteristics of the GLACIER Study (N=5,160)

Variable	N	Mean	SD	Median
Age (year)	5160	49.0	8.6	50.1
Weight (kg)	5160	73.1	12.9	72.0
Height (cm)	5160	169.1	8.7	168.0
BMI (kg/m2)	5160	25.5	3.8	25.0
Fasting glucose (mmol/l)	5143	5.3	0.7	5.3
2-h glucose (mmol/l)	5012	6.6	1.5	6.5
HDL-C (mmol/l)	2721	1.4	0.5	1.4
LDL-C (mmol/l)	2706	4.2	1.2	4.2
Triglycerides (mmol/l)	4617	1.4	0.8	1.2
Total cholesterol (mmol/l)	5142	5.9	1.2	5.9
n-3 PUFA (g/day)	5160	1.7	0.8	1.6
n-6 PUFA (g/day)	5160	7.2	3.6	6.4
PUFA (g/day)	5160	9.9	4.7	8.8
TEI (cal/day)	5159	1799.4	569.1	1689.7

BMI - body mass index; HDL-C - high-density lipoprotein cholesterol; LDL-C - low-density lipoprotein cholesterol; PUFA - polyunsaturated fat intake; SD - standard deviation; TEI - total energy intake.

Dietary PUFA × SNP interactions (Step 1)

Several gene × dietary PUFA interactions were detected (P<0.05) for the six candidate SNPs across the eight traits. The strongest interaction effects were observed between rs174602 and n–3 PUFA intake on total cholesterol (P=0.001, Figure 2) and between rs174570 and n-6 PUFA intake on fasting glucose (P=0.005, S1 Figure). While rs174602 showed a negative association with total cholesterol (β=-0.09 mmol/l per C allele, 95% CI -0.17; -0.01, P=0.02) among those with low n-3 PUFA intake, no association was observed (β = -0.03 mmol/l per C allele, 95% CI -0.11; 0.05, P=0.50) among those with high n-3 PUFA intake. While rs174570 showed a negative association with fasting glucose (β = -0.06 mmol/l per T allele, 95% CI -0.11; -0.01, P = 0.01) among participants with a high n-6 PUFA intake, no association (β = 0.02 mmol/l per T allele, 95% CI 0.03; 0.95, P = 0.34) was observed in participants who reported low n-6 PUFA intake. Repeating the analysis with inverse normalized traits as dependent variables did not materially change the results. Results for all interaction tests are presented in S2 Table.

Figure 2

Total cholesterol levels stratified by median n-3 PUFA intake and rs174602 genotypes in the GLACIER Study. The bars represent adjusted means from the generalized linear regression models described in Methods. Error bars represent 95% confidence intervals.

Gene-centric analyses (Step 2)

As shown in Table 2, a significant interaction effect (P=0.005; q=0.05) was observed between the FADS gene cluster and n-3 PUFA intake on triglycerides with genetic main effects as random, and was followed up with haplotype-based analyses to refine the signal.

Table 2

Gene-centric interaction analysis (N=5,160)

Environment	Phenotype	PINT_FIX	PINT_RAND	PJOINT	qINT_FIX	qINT_RAND	qJOINT
PUFA	Weight (kg)	0.391	0.149	0.018	1.00	1.00	0.18
PUFA	BMI (kg/m2)	0.448	0.237	0.031	1.00	1.00	0.31
PUFA	HDL-C (mmol/l)	0.907	0.952	0.312	1.00	1.00	1.00
PUFA	LDL-C (mmol/l)	0.324	0.325	0.150	1.00	1.00	1.00
PUFA	Triglycerides (mmol/l)	0.168	0.143	0.164	1.00	1.00	1.00
PUFA	Total cholesterol (mmol/l)	0.337	0.470	0.068	1.00	1.00	0.69
PUFA	Fasting glucose (mmol/l)	0.070	0.203	0.332	0.71	1.00	1.00
PUFA	2-h glucose (mmol/l)	0.185	0.149	0.290	1.00	1.00	1.00
n-3 PUFA	Weight (kg)	0.691	0.229	0.024	1.00	1.00	0.24
n-3 PUFA	BMI (kg/m2)	0.724	0.342	0.033	1.00	1.00	0.33
n-3 PUFA	HDL-C (mmol/l)	0.888	0.783	0.310	1.00	1.00	1.00
n-3 PUFA	LDL-C (mmol/l)	0.549	0.591	0.154	1.00	1.00	1.00
n-3 PUFA	Triglycerides (mmol/l)	0.121	0.005	0.014	1.00	0.05	0.14
n-3 PUFA	Total cholesterol (mmol/l)	0.448	0.600	0.068	1.00	1.00	0.69
n-3 PUFA	Fasting glucose (mmol/l)	0.208	0.122	0.312	1.00	1.00	1.00
n-3 PUFA	2-h glucose (mmol/l)	0.218	0.165	0.399	1.00	1.00	1.00
n-6 PUFA	Weight (kg)	0.459	0.164	0.021	1.00	1.00	0.21
n-6 PUFA	BMI (kg/m2)	0.533	0.256	0.029	1.00	1.00	0.29
n-6 PUFA	HDL-C (mmol/l)	0.997	0.999	0.307	1.00	1.00	1.00
n-6 PUFA	LDL-C (mmol/l)	0.248	0.181	0.125	1.00	1.00	1.00
n-6 PUFA	Triglycerides (mmol/l)	0.155	0.050	0.079	1.00	0.51	0.80
n-6 PUFA	Total cholesterol (mmol/l)	0.277	0.137	0.042	1.00	1.00	0.42
n-6 PUFA	Fasting glucose (mmol/l)	0.047	0.166	0.281	0.48	1.00	1.00
n-6 PUFA	2-h glucose (mmol/l)	0.219	0.142	0.281	1.00	1.00	1.00

BMI - body mass index; HDL-C - high-density lipoprotein cholesterol; LDL-C - low-density lipoprotein cholesterol; PUFA - polyunsaturated fatty acids intake; P - P value for interaction effects treating genetic main effects as fixed; P - P value for interaction effects treating genetic main effects as random; P - P values for joint test of genetic main effects and gene-environment interactions; q - q value for interaction effects treating genetic main effects as fixed; q - q value for interaction effects treating genetic main effects as random; q - q values for joint test of genetic main effects and gene-environment interactions

Haplotype block interaction analyses (Step 3)

From over 700 variants in the FADS gene cluster, 442 variants passed quality control and were used to reconstruct haplotypes (S3 Table). Haploview inferred 30 haplotype blocks across the FADS cluster, and the number of variants in the haplotype blocks ranged from 2 to 63. Focusing on the finding from Step 2, we performed haplotype block based n-3 PUFA interaction analyses in relation to triglyceride levels (Table 3). Although none of the haplotypes surpassed the pre-defined P<0.002 threshold, haplotype blocks 12 (P = 0.01; q = 0.32), 16 (P = 0.01; q = 0.30) and 21 (P = 0.01; q = 0.34) showed tentative signals for interactions; we thus functionally annotated variants within these three haplotype blocks.

Table 3

Haplotype block × n-3 PUFA interactions in relation to triglyceride levels (N=5,160)

Haplotype	PINT_FiX	PINT_RAND	PJOINT	qINT_FiX	qINT_RAND	qJOINT
Haplotype 1	0.099	0.097	0.188	1.00	1.00	1.00
Haplotype 2	0.334	0.339	0.332	1.00	1.00	1.00
Haplotype 3	0.133	0.095	0.154	1.00	1.00	1.00
Haplotype 4	0.092	0.068	0.137	1.00	1.00	1.00
Haplotype 5	0.109	0.029	0.033	1.00	0.87	1.00
Haplotype 6	0.808	0.799	0.775	1.00	1.00	1.00
Haplotype 7	0.318	0.356	0.178	1.00	1.00	1.00
Haplotype 8	0.057	0.034	0.073	1.00	1.00	1.00
Haplotype 9	0.176	0.170	0.313	1.00	1.00	1.00
Haplotype 10	0.263	0.226	0.406	1.00	1.00	1.00
Haplotype 11	0.645	0.650	0.504	1.00	1.00	1.00
Haplotype 12	0.007	0.011	0.017	0.21	0.32	0.51
Haplotype 13	0.088	0.087	0.166	1.00	1.00	1.00
Haplotype 14	0.199	0.166	0.315	1.00	1.00	1.00
Haplotype 15	0.679	0.652	0.879	1.00	1.00	1.00
Haplotype 16	0.017	0.010	0.019	0.52	0.30	0.57
Haplotype 17	0.587	0.602	0.753	1.00	1.00	1.00
Haplotype 18	0.329	0.202	0.358	1.00	1.00	1.00
Haplotype 19	0.374	0.355	0.076	1.00	1.00	1.00
Haplotype 20	0.206	0.231	0.407	1.00	1.00	1.00
Haplotype 21	0.017	0.011	0.027	0.50	0.34	0.81
Haplotype 22	0.182	0.148	0.270	1.00	1.00	1.00
Haplotype 23	0.320	0.302	0.512	1.00	1.00	1.00
Haplotype 24	0.535	0.538	0.334	1.00	1.00	1.00
Haplotype 25	0.018	0.021	0.028	0.55	0.64	0.84
Haplotype 26	0.298	0.277	0.048	1.00	1.00	1.00
Haplotype 27	0.461	0.573	0.451	1.00	1.00	1.00
Haplotype 28	0.119	0.119	0.224	1.00	1.00	1.00
Haplotype 29	0.433	0.423	0.666	1.00	1.00	1.00
Haplotype 30	0.559	0.562	0.814	1.00	1.00	1.00

PUFA - polyunsaturated fat intake; P - P value for interaction effects treating genetic main effects as fixed; P - P value for interaction effects treating genetic main effects as random; P - P values for joint test of genetic main effects and gene-environment interactions; q - q value for interaction effects treating genetic main effects as fixed; q - q value for interaction effects treating genetic main effects as random; q - q values for joint test of genetic main effects and gene-environment interactions

Functional annotation (Step 4)

Variants from FADS haplotype blocks 12 (n=6), 16 (n=12) and 21 (n=19) and the six candidate SNPs (in total, 43 variants) were functionally annotated (S4 Table) using ANNOVAR. ChromHMM predictions based on nine human cell types demonstrated tissue-dependent functionality for the haplotype blocks and most variants. For instance, the overwhelming majority of variants (18 out of 19) in haplotype block 21 show an enhancer state in the K562 blood (leukemia) cell line, while in the other cell lines, the same variants demonstrate predominantly repressed or inactive chromatin states. As FADS1 and FADS2 are highly expressed in the liver, we assessed chromatin states in the HepG2 liver carcinoma cell line. Here a cluster of variants in haplotype block 21 was enriched with enhancer states, two neighbouring variants (rs187943834, rs117518711) showed promoter states and one variant (rs7115739; a previously reported variant (24)) showed an insulator state. We observed fewer regulatory states among variants in haplotype blocks 12 and 16 in HepG2 cells. Among the other candidate variants, rs174570 variant showed weak promoter, active promoter, or strong enhancer states across multiple cell lines, while the rs174602 variant showed weak transcription elongation state across all but two cell lines.

To validate the performance of our 3DSNP prioritization pipeline, we first showed that the established causal variant (rs1421085) at FTO (as reported by Claussnitzer et al. (49)) can be detected using the known tagging variant (rs9930506) (S1 Text). The 3DSNP functional annotation analyses of the candidate FADS SNPs is summarised in S5 Table and includes information about the functionality score, distal interacting genes, and the selected leading variants. From the list of tagging variants (prior to the selection of leading variants and ranking by functionality score), rs174570 ranked the highest (score=96.7, driven by promoter status, Figure 3). However, after re-ranking by the leading variant, the FADS2 intronic rs5792235 deletion (CA/C) showed the highest functionality score (218.8), which is 18x higher than that of its proxies, rs174599 (score=12.13; r2=0.83) and rs174601 (score=11.69; r2=0.82), driven by promoter and TFBS status (Figure 3). Eight variants in LD with at least one of the original variants ranked higher than rs174570 (score>96.7; range of 114.0-218.0), many of them driven by high promoter, TFBS and/or enhancer scores. The Circos plot based on rs5792235 (Figure 4) demonstrates 17 proximal (including the three FADS genes) and distal interacting genes (these genes are described in brief in S2 Text), many of which are known to be linked with cardiovascular traits. These findings are supported by HaploReg eQTL annotations (S6 Table), which demonstrate cis and trans eQTL evidence for rs5792235 and multiple other variants in the region. We curated the rs5792235 variant’s 65 distal interacting variants and their associated phenotypes in S7 Table; multiple variants show associations with glycaemic and lipid traits in previous GWAS, many of which are also QTLs for serum lipid traits (S6 Table). Gene ontology enrichment analysis of the 17 genes in distal interaction with rs5792235 revealed enrichment in the unsaturated fatty acid metabolic process (P=5.7×10-3); however, no pathways were identified when FADS1, FADS2 and FADS3 were removed from the analysis. The Circos plot of rs174570 was identical to the one generated for rs5792235, while the one based on rs174602 revealed less interacting nodes.

Figure 3

Radar charts of rs174570 and rs5792235. The six axes of the hexagon represent functionality levels (0-100) for enhancer status, promoter status, transcription factor binding site, motifs, evolutionary conservation and 3D interacting genes, as suggested by 3DSNP.

Figure 4

Circos plot based on rs5792235 (shown in red text). The plot shows rs5792235 and its proxies (shown in black text around rs5792235). From outer to inner, the circles represent ChromHMM chromatin states, annotated genes (green), histone modification set (red), transcription factor set (blue), current variant (rs5792235) and associated variants, and 3D chromatin interactions, respectively. The three circles in the histone modification set are H3K4me1, H3K4me3, H3K27ac, and the three circles in the transcription factor set are CTCF, CEBPB and CEBPD (adapted from [42]). Color schemes for the ChromHMM chromatin states and the 3D interaction loops can be found at http://biotech.bmi.ac.cn/3dsnp/documentation/tutorials/.

Discussion

We systematically assessed gene × dietary PUFA interactions at the FADS1-FADS2-FADS3 gene cluster across multiple obesity-related traits in a population isolate from the far north of Sweden. There is strong molecular and genetic data implicating FADS gene variation and action in lipid metabolism in several populations. In Greenlandic Inuit, the strongest associations were observed for FADS3 rs7115739 and FADS2 rs174570 and anthropometric traits (24). These signals persisted when conditioning on other variants in the FADS gene cluster, suggesting independent association signals and multiple causal variants in the region. This observation motivated us to study the entire FADS region in gene-centric and haplotype-centric analyses and undertake functional annotations. Indeed, results emanating from the genomic annotations reveal multiple potential causal variants in the FADS region.

Specifically, we investigated gene × environment interactions between PUFA intake (and its components, n-3 and n-6 PUFA intake) and the six candidate SNPs (rs74771917, rs3168072, rs12577276, rs7115739, rs174602, and rs174570); these candidate variants have been shown to be under selective pressure (24).

FADS1 and FADS2 represent rate limiting steps in the fatty acid metabolism pathway. Multiple genome-wide associations studies have consistently demonstrated that genetic variants at FADS1 and FADS2 are associated with plasma and tissue levels of arachidonic acid and EPA (50) and commonly measured blood lipids, such as triglycerides, total cholesterol, LDL-C and HDL-C (14). Here we show that the variant rs174570 at FADS2 appears to modify the association between n-6 PUFA intake and fasting glucose levels. Previous studies have investigated the role of FADS2 variants and glycaemic traits; for example, Corpeleijn et al. reported that the increased activity of the Δ-5 desaturase and decreased activity of the Δ-6 desaturase are associated with reduced insulin resistance in a fasted state, and this association is modified by total fat intake (51). In a gene × environment interaction analysis focusing on 18 FADS SNPs, Cormier et al. observed an interaction between rs174570 and n-3 PUFA intake on insulin sensitivity, but no interaction was observed on fasting glucose levels (52). Fumagalli et al. reported a putative association between the rs174570 T allele and decreased fasting glucose (24). Consistently, we found here that the minor T allele of the same SNP was associated with decreased fasting glucose levels, but only when n-6 PUFA intake was high. Cormier et al. also reported a statistically significant interaction between the rs174602 variant and n-3 PUFA intake on insulin sensitivity (52).

In our study, statistical interactions were observed between rs174602 and n-3 PUFA intake in levels of total cholesterol and LDL-C, but there was no evidence of interactions in glycaemic traits. Buckley et al. have also reported tentative evidence of interactions between the FADS1 rs174594 and lipid traits (53). Notably, in the case of the two observed interactions, the associations between the FADS variants and the outcomes are only apparent in the subpopulations with lower n-3 and higher n-6 PUFA intakes (less healthy diets), respectively. In these “unhealthful” strata of GLACIER, the major (compared with minor) alleles are associated with higher total cholesterol and fasting glucose levels. The fact that the major alleles in GLACIER are the minor alleles in the Inuit suggests the presence of a cross-over interaction, where the major allele raises blood cholesterol and glucose in the setting of an unhealthy diet (in GLACIER), and lowers these metabolites in the presence of a healthy diet (in Inuits), relative to the minor allele. Interactions of this nature are consistent with the hypothesis described by Neel (54) with regard to the storage and metabolism of metabolic substrate. Both rs174602 and rs174570 are reported with >10 within loop (proximal) or anchor-to-anchor (distal) interactions, which allows us to speculate on the putative functional background of the observed interactions. Although all these loci are curated in S2 Text, we highlight here that while some of these genes have yet unclear functions, some have putative associations with lipid (BEST1, FEN1, MIR1908, RAB3IL1), glycaemic (MIR1908, MYRF, SYT7), obesity traits (DAGLA, RPLP0P2) or basic cellular functions (INCENP). It is also possible, that following environmental, dietary triggers, different FADS variants interact with different target genes in the region, thereby exerting a complex metabolic phenotype.

Recently, results from two CHARGE Consortium reports indicated no interactions between FADS variants and dietary PUFAs in relation to two clinical endpoints, type 2 diabetes (55) and coronary heart disease (56). While these large-scale meta-analyses include larger sample sizes than the GLACIER cohort, the reported lack of interactions might be explained by the analyzed disease outcomes having diverse etiology, using ancestrally different populations, and the dilution of real effects/interactions due to within and between study heterogeneity, as discussed elsewhere (57).

Although the Inuit study indicates that variation within the FADS locus may modulate the effects of dietary fats in energy metabolism (24), the specific variants may not correspond across Greenlandic and northern Swedish populations, owing to their diverse evolutionary backgrounds. Thus, we extended the single variant analyses to include variation across the entire FADS1-FADS2-FADS3 gene cluster using the MetaboChip array. Although sequencing the region would likely provide better coverage and allow analyses of structural elements, such as indels and copy number variations, in terms of MetaboChip, there is an above-average coverage of the FADS region, as this locus was chosen for fine-mapping purposes by the developers of the array (35).

Due to the high multiple testing burden, the abundance of low-frequency variants and the small magnitude of interaction effects that were anticipated, we would likely have been underpowered to perform a region-wide single variant interaction analysis, and so this was not attempted. To address this, we undertook a top-down, sequencial approach (described above), where we first assessed the whole region for gene × environment signals to identify exposures and subsequently narrowed down our results to haplotype blocks and single variants through haplotype-based analyses and functional annotation, respectively. We detected a gene cluster interaction with n-3 PUFA intake on triglycerides concentrations, which we scrutinized further by formally testing which haplotype blocks might drive the observed gene-centric interaction effects. Three haplotype blocks (12, 16, and 21) showed suggestive interactions with n-3 PUFA intake on triglyceride levels.

We chose to employ 3DSNP to undertake detailed in silico functional annotations of selected FADS variants, as this software utilizes publically available HiC datasets to visualize long-range 3D chromosomal interactions, which are able to reveal distal regulatory potential of genomic regions—a key functional feature that many other tools do not include. Inspired by the findings of Claussnitzer et al. (49), we first tested our 3DSNP functional analysis pipeline on the GWAS-identified FTO rs9930506 variant and the subsequently reported causal FTO variant, rs1421085. Using 3DSNP, we successfully validated rs1421085 as the causal variant by demonstrating long-range interactions with IRX5 and a 15x higher functional score compared to rs9930506. We then applied the pipeline on 43 variants from FADS, selected from our previous analyses. Functional annotation of these variants revealed potentially important regulatory signals and topological interactions in the region. The variant showing the strongest evidence of functionality was rs5792235, an intronic FADS2 deletion, with 17 potential interacting genes and high evidence for localizing in a promoter region and a TFBS (Figures 3 and 4). Further annotations and visualizations suggested multiple regulatory variants in the region. Distal interactions of rs5792235 reveal multiple interacting loci with previously demonstrated associations with cardiometabolic traits, mainly lipids and glycaemia (58–60). Multiple annotated variants in the region demonstrate QTL evidence for serum lipid traits and eQTL evidence for genes in the proximity.

Functional annotation of the lead SNPs and haplotypes identified FADS2 rs5792235 as a probable causal SNP. Interestingly, the deletion rs5792235 is common in European (35%), East Asian (57%) and mixed American (60%) populations, while lower frequencies are seen in South Asian (14%) and African (10%) populations (46). rs5792235 is not captured by any genotyping array, which may explain why this has not previously been implicated in lipid variation.

A limitation of this study is that the dietary variables were obtained from FFQs. While validated FFQs are often used in large epidemiological studies such as the GLACIER Study, self-report methods are prone to recall bias. However, in the context of these gene-diet interaction analyses, where there is a strong a prior hypothesis, the most serious consequence of reporting bias may be the underestimation of the interaction effects. We also recognize that the in silico functional validations provide only suggestive evidence of functionality, and in vitro analyses, which go beyond the scope of this paper, would be necessary to determine function with very high certainty.

In conclusion, we provide evidence for gene × environment interactions on a gene variant level for fasting glucose and total cholesterol levels, and on a whole gene level for triglycerides. Through functional annotation and characterization of topological genomic interactions, we identified the intronic rs5792235 FADS2 variant as a potential causal variant (deletion) in the region, as well as multiple interacting genes proximal to FADS2. While contemporary studies often aim to identify single causal master regulators for GWAS-associated loci, our data suggest that it is more likely that multiple causative variants with various regulatory roles in the FADS1-FADS2-FADS3 gene cluster regulate the gene-diet interactions of relevance here.

Supplementary Material

Supplementary information is available at the International Journal of Obesity’s website.

Fasting glucose levels stratified by median n-6 PUFA intake and rs174570 genotypes in the GLACIER Study. The bars represent adjusted means from the generalized linear regression models described in Methods. Error bars represent 95% confidence intervals.