Meta-analysis of genome-wide association studies for body fat distribution in 694 649 individuals of European ancestry

1Approximately 39% of adults worldwide is either overweight or obese and is at increased risk of metabolic disease. While higher adiposity increases morbidity and mortality, epidemiological studies indicate that the location and distribution of excess fat within particular depots are more informative than general adiposity for predicting disease risk. Independent of their overall body mass index (BMI), individuals with higher central adiposity have increased risk of cardiometabolic diseases, including type 2 diabetes (T2D) and stroke; in contrast, individuals with higher gluteal adiposity have lower risk of such outcomes. Previous studies indicate that fat distribution, as assessed by waist-to-hip ratio (WHR), is a trait with a strong heritable component, independent of overall adiposity (measured by BMI), with twin-based heritability estimates ranging between 30% and 60%, and narrow-sense heritability estimates have been estimated at ∼50% in women and ∼20% in men. The most recent genome-wide association study (GWAS) in 224 459 samples implicated 49 loci associated with WHR adjusted for BMI (WHRadjBMI), and recent Mendelian randomization studies using known WHR-associated genetic variants showed putative causal effects of higher WHR on T2D and coronary artery disease independent of BMI.

2With the goal of pinpointing genetic variants associated to body shape and fat distribution and motivated by the recent release of genetic data from half a million individuals, we performed a meta-analysis of WHRadjBMI. WHRadjBMI is an easily-measured fat distribution phenotype that correlates well with imaging-based fat distribution measures. We performed GWAS of WHRadjBMI in the UK Biobank data set, a collection of 484 563 samples with densely imputed genotype data, using a linear mixed model (LMM) to account for relatedness and ancestral heterogeneity. We then combined the results with publicly available GWAS data generated by the Genetic Investigation of ANthropometric Traits (GIANT) consortium for the same phenotype (Table 1 and Methods), resulting in a meta-analysis of 694 649 samples (Table 1) and ∼27.4M single-nucleotide polymorphisms (SNPs) (Methods). As a sensitivity analysis and to evaluate the robustness of our results, we also performed a GWAS of WHR unadjusted for BMI (Table 1).

3We identified 346 loci (300 novel) containing 463 independent signals associated with WHRadjBMI [P < 5 × 10−9, to account for the denser imputation data; Methods, Supplementary Material, Table 1 and Supplementary Material, Fig. 1]. The linkage disequilibrium (LD) score regression intercept (1.035) of the meta-analysis results indicated that the observed enrichment in genomic signal was due to polygenicity and not confounding (Supplementary Material, Table 2). Of the 300 novel signals, 234 (78%, Pbinomial < 1 × 10−7) were directionally consistent in an independent dataset with a relatively small sample size (N = 7721), and signals were consistent in several sensitivity checks (Supplementary Material, Tables 3–5 and Supplementary Material, Figs 2–3). Combined, these variants explained ∼3.9% of the variance in WHRadjBMI in the independent study (Methods and Table 1). We constructed a weighted polygenic risk score (PRS) using the 346 index SNPs discovered in the combined meta-analysis and tested this score in the same independent study. The 5% of individuals carrying the most WHRadjBMI-raising alleles were 1.62 times more likely to meet the WHR threshold used to define metabolic syndrome than the 5% carrying the fewest (consistent with the results obtained from unweighted polygenic score; Methods). The WHRadjBMI of people in the top 5% of the PRS was 1.05 and 1.06 times greater in men and women, respectively, compared to those in the bottom 5% of the PRS.

4To investigate the potential for collider bias resulting from conditioning WHR on BMI, we investigated the behavior of WHRadjBMI-associated SNPs in GWAS of WHR (without adjustment for BMI) and BMI alone. We found that the majority of WHRadjBMI signals identified has genuine effect on body shape and that any bias caused by adjusting WHR for a correlated covariate (that is, BMI) was minimal. Of the 346 index variants, 311 associated with stronger standard deviation effect sizes for WHR (unadjusted) than with standard deviation effect sizes for BMI (Supplementary Material, Table 3 and Supplementary Material, Fig. 4). This observation also indicates that the WHR association is unlikely to be secondary to the known effect of higher BMI resulting in higher WHR. Furthermore, the common SNP associated with the largest known effect on BMI, that in the FTO gene, was not associated with WHRadjBMI (rs1421085, P = 0.40) despite a very strong association with WHR (P = 4 × 10−118)). Finally, carrying each additional (weighted) WHRadjBMI-raising allele was associated with an increase in WHRadjBMI of 0.0199 SD (P = 6 × 10−62; adjusted R2 = 4%), an increase in WHR of 0.011 SD (P = 3 × 10−20; adjusted R2 = 0.12%) and a decrease in BMI of 0.004 SD (P = 1.4 × 10−3; adjusted R2 = 0.13%) in our independent dataset, consistent with the results obtained from an unweighted polygenic score (Methods).

5Given the sex dimorphism of fat distribution in humans, previously shown to have a genetic basis, we next performed meta-analyses of WHRadjBMI in women and men separately (Table 1 and Supplementary Material, Fig. 5). We found SNP-based heritability of WHRadjBMI, estimated using the restricted maximum likelihood (REML) method implemented in BOLT-REML (Methods), to be stronger in women ( = 25.6%) compared to men ( = 16.7%, Pdifference = 9 × 10−85; Table 1, Supplementary Material, Table 6, and Equation 2). In addition to the heritability dimorphism, and in keeping with previous studies, we found signatures of sex dimorphism among associated loci: a total of 266 loci associated with WHRadjBMI in women, compared to 91 loci in men (P < 5 × 10−9). Genome-wide, SNP effects on WHRadjBMI were strongly correlated between men and women [LD Score rg = 0.514 (s.e. = 0.019), P = 3.43 × 10−159], but the consistency between the effect size of 266 female index SNPs on WHRadjBMI in women and men (adjusted R2 = 51%) was greater than the consistency between the effect size of 91 male index SNPs on WHRadjBMI in men and women (adjusted R2 = 9%). Of all associated index SNPs (P < 5 × 10−9 in the combined or sex-specific analyses), 105 SNPs were sex-dimorphic [Pdiff < 3.3 × 10−5; and Methods]. Variants discovered in the combined sex analysis will be enriched for those with similar effects in each sex, while variants discovered in sex-specific analyses will be enriched for those with differing effects between sexes. In the absence of any sex-specific effects, we would only expect a slight shift towards stronger associations in women due to the larger available sample size in that analysis. However, we observed that of the 105 sex-dimorphic signals, 97 (92.4%) showed stronger effects in women compared to men (Fig. 1, Supplementary Material, Fig. 6 and Methods). Scanning genome-wide for sex-dimorphic SNPs (Pdiff < 5 × 10−9), regardless of their association P-values in the sex-specific analyses, we identified 61 sex-dimorphic SNPs after LD-based clumping (r2 < 0.05). Of these, 19 (31.1%) overlapped with the sex-dimorphic and genome-wide significant loci, and 54 (88.5%) had stronger effect in women than in men (Supplementary Material, Information).

6Previous studies have shown that in addition to redistributing body fat, some WHRadjBMI variants are also associated with total body fat percentage (BF%). Of relevance to the biology of adipose tissue storage capacity, these studies have shown that these pleiotropic associations can occur in both directions: some alleles associated with higher WHRadjBMI are associated with higher total BF%, while others are associated with lower BF%. To test the hypothesis that alleles associated with higher WHRadjBMI could have pleiotropic effects on total BF% and that these effects could occur in both directions, we next investigated whether 346 index variants associated with WHRadjBMI also associated with BF%. Of the 59 of 346 variants associated with BF% in 443 001 European-ancestry UK Biobank individuals (P < 0.05/346 = 1.44 × 10−4), 25 SNPs associated with higher WHR and higher BF%, while 34 SNPs associated with higher WHR but lower BF% (Fig. 2). These findings indicate that WHR-increasing alleles do not strictly influence BF% in one direction but rather can associate with either higher or lower BF%, yielding biological insight beyond the known epidemiological correlation between BF% and WHR. Additionally, a large proportion (29%) of WHRadjBMI index SNPs with a stronger effect in women had a BF% phenotype in men: 28 of the 97 female-specific WHRadjBMI SNPs were associated with BF% in men, and 25 were associated with BF% in women (P < 0.05/105 = 4.8 × 10−4; Supplementary Material, Fig. 7). These variants appear to alter total BF% in men and women to a similar extent but distribute body fat between the upper and lower body to a much greater extent in women (Supplementary Material, Table 7–9 and Supplementary Material, Fig. 7). Finally, we tested the index SNPs from each of the meta-analyses (combined and sex-specific) in a recent GWAS of computed tomography (CT) and magnetic resonance imaging (MRI) image-based measures of ectopic and subcutaneous fat depots. Adjusting for the three sample groups and the eight depots examined in the imaging-based GWAS (P < 0.05/24 = 2.1 × 10−3), the alleles associated with higher WHRadjBMI were collectively associated with lower measures of subcutaneous fat and higher measures of visceral fat, including pericardial and visceral adipose tissue (Supplementary Material, Fig. 8).

7In a meta-analysis of nearly 700 000 individuals, we have increased the number of loci associated to WHRadjBMI by more than 7-fold. Of all the detected signals, 105 are sex-dimorphic, consistent with previous findings. While we have performed the largest meta-analysis of a measure of body fat distribution to date, a number of limitations remain. First, the substantially larger number of signals with a stronger effect in women compared to men may be influenced by the reduction in power (proportional to the product of sample size and SNP heritability) in the men-only analysis (Table 1) compared to the women-only analysis. Despite the power difference in the sex-specific analyses, we would not expect the difference to result in 92% of signals conferring a stronger effect in women. Second, our replication sample was too small (∼1% of the discovery) to formally replicate individual SNP associations, but the fact that 78% of the 300 previously unknown index associations showed consistent direction of effect suggests a low false-positive rate. Finally, our meta-analysis focused only on European-ancestry samples. Given the very different body fat distributions observed across ancestral groups, and the very different risks of adiposity-related disease across populations, studies in non-Europeans are urgently needed.

8In summary, the genetic variants and loci identified by this meta-analysis will likely provide starting points for further understanding the biology of body fat distribution and its relationship with disease.

9Code and data related to this project, including summary-level data from the meta-analyses, can be found online at https://github.com/lindgrengroup/fatdistnGWAS.

10To generate phenotypes for the WHR and WHRadjBMI analyses in the UK Biobank data (Supplementary Material, Table 10), we followed a phenotype conversion consistent with that performed in previous efforts investigating WHR and WHRadjBMI by the GIANT consortium.

11Using phenotype information from UK Biobank, we divided waist circumference by hip circumference to calculate the WHR measure, and then regressed the WHR measure on sex, age at assessment, age at assessment squared and assessment centre. To generate the WHRadjBMI phenotype, we followed the same procedure and included BMI as an additional independent variable in the regression. We performed rank inverse normalization on the resulting residuals from the regression (Supplementary Material, Fig. 9) and used these normalized residuals as the tested phenotype in downstream genome-wide association testing. To generate phenotypes for the sex-specific analyses, we followed this same procedure but ran the regressions in sex-specific groups.

12We conducted genome-wide association testing in the second release (June 2017) version of the UK Biobank data; this release did not contain the corrected imputation at non-Haplotype Reference Consortium [HRC] sites, and we therefore subset all of the SNP data down to HRC SNPs only. The UK Biobank applied quality control to samples and genotypes and imputed the resulting genotype data using sequencing-based imputation reference panels. We performed all of our genome-wide association testing and downstream analyses on the publicly available imputation data (released in bgen format).

13We excluded samples as suggested by the UK Biobank upon release of the data (Supplementary Material, Table 11). Sample exclusions included samples with genotype but no imputation information, samples with missingness >5%, samples with mismatching phenotypic and genotypic sex and samples that have withdrawn consent since the initiation of the project.

14We implemented all GWAS in BOLT-LMM, which performs association testing using an LMM. To run, BOLT-LMM requires three primary components: the (imputed) genotypic data for association testing; a reference panel of LD scores per SNP, calculated using LD Score Regression; and genotype data used to approximate a genetic relationship matrix (GRM), which is the best method available in this sample size to account for all forms of relatedness, ancestral heterogeneity in the samples and other (potentially hidden) structure in the data.

15We performed sensitivity testing (Supplementary Material, Information, Supplementary Material, Tables 12–13 and Supplementary Material, Fig. 10) using three LD Score reference datasets and four SNP sets to construct the GRM. For our final GWAS, we used LD scores calculated from a randomly selected 9748 unrelated UK Biobank samples (∼2% of the full UK Biobank sample set; Supplementary Material, Information) and a GRM constructed using imputed SNPs with imputation info score > 0.8, minor allele frequency (MAF) > 1%, Hardy Weinberg P-value > 1 × 10−8, genotype missingness < 1%, after converting imputed dosages to best-guess genotypes, LD pruned at a threshold (r2) of 0.2, and excluding the major histocompatibility complex, the lactase locus and the inversions on chromosomes 8 and 17 (Supplementary Material, Information).

16For genome-wide association testing, we used BOLT-LMM to run an LMM. We tested SNPs with imputation quality (info) > 0.3, MAF > 0.01% (equivalent to ∼100 copies of the minor allele in the full sample), and only those single-nucleotide variants and SNPs represented in the HRC imputation reference panel. We used only the standard LMM implementation (i.e. infinitesimal model, using l mm) in BOLT-LMM (Supplementary Material, Figs 11–12); we did not run association testing using a non-infinitesimal model. The only covariate used in the LMM was the SNP array used to genotype sample; we included no other covariates.

17After association testing, we looked at known SNPs already reported in WHR, WHRadjBMI and BMI. At the previously-described loci, we checked correlation of frequency, beta, standard error and -log10(P-value) between our UK Biobank GWAS and the previous GWAS results (Supplementary Material, Fig. 13). Additionally, we estimated genomic inflation (lambda) and the LD Score Intercept to check if the P-values were well calibrated (Supplementary Material, Table 2); calculations were performed using the LD Score software (https://github.com/bulik/ldsc).

18We downloaded summary-level results from previous meta-analyses of WHR and WHRadjBMI (https://portals.broadinstitute.org/collaboration/giant/index.php/GIANT_consortium_data_files and Supplementary Material, Information) performed by the GIANT consortium. Marker names in both the GIANT data and UK Biobank were lifted over to their dbSNP151 identifier. We additionally renamed markers as ‘rsID:A1:A2’ (where A1 was the tested allele in UK Biobank) to avoid ambiguity at multiallelic SNPs in the UK Biobank data. As the GIANT data was imputed with HapMap 2 data (hg18), we additionally lifted chromosomal positions to hg19 for this data. SNPs with a frequency difference > 15% between GIANT and UK Biobank were removed from the data (Supplementary Material, Fig. 14).

19We performed inverse variance-weighted fixed effects meta-analysis in METAL. To estimate LD score intercepts and genomic inflation (lambda) for the meta-analysis results, we first estimated LD scores from the same samples used to estimate the LD score reference for BOLT-LMM. LD scores were only estimated at high-quality SNPs (using the same criteria as used for SNPs included in the GRM in BOLT-LMM, but without applying a MAF threshold; Supplementary Material, Information). We then calculated LD Score Regression intercepts and lambda with the LDSC software.

20As an additional quality control check, we reran all of our GWAS using two different subsets of the UK Biobank samples: the unrelated samples only and the unrelated white British samples only. These subsamples were selected to test if our initial UK Biobank–wide GWAS was confounded by either relatedness or ancestral heterogeneity. After running these GWAS, we meta-analyzed the results with the existing GIANT summary-level data and checked the concordance of our signals (Supplementary Material, Figs 2–3).

21To identify genomic loci (i.e. genomic windows) containing independent association signals, we first constructed a reference dataset of best-guess genotypes from 20 275 unrelated UK Biobank samples (equivalent to 5% of the unrelated sample). We converted imputed dosages of SNPs with info score > 0.3 and MAF > 0.001% to best-guess genotypes using PLINK (version 1.9) and a conversion threshold (hard-call-threshold) of 0.1 (Supplementary Material, Information). SNPs with missingness > 5% after conversion or Hardy–Weinberg equilibrium P < 1 × 10−7 were removed.

22We then used the PLINK ‘clumping’ algorithm to select top-associated SNPs (P < 5 × 10−9) and identify all SNPs in LD (r2 > 0.05) with the top associated SNP and ±5 Mb away. We determined the genomic span of each LD-based clump and added 1 kb up- and downstream as buffer to the region. If any of these windows overlapped, we merged them together into a single (larger) locus. As a sensitivity analysis, we ran clumping also using a smaller genomic window to calculate LD (±2 Mb); the results were effectively unchanged, as <5 loci appeared independent using the ±2 Mb window but were found to correlate using ±5 Mb windows. Therefore, we report loci using the ±5 Mb window.

23To identify index and secondary signals within each of the clumping-based loci, we ran proximal joint and conditional analysis as implemented in the Genome-wide Complex Trait Analysis (GCTA) software. We ran this model (cojo-slct) using the summary-level data within each locus, the LD reference panel constructed from UK Biobank data and also used for the locus ‘clumping’, and setting genome-wide significance with P < 5 × 10−9.

24We used an independent dataset EXTEND (7721 individuals of European descent collected from South West England, Supplementary Material, Table 14) to validate our findings. We extracted the index SNPs from the HRC-imputed genotypes. To generate the WHRadjBMI variable, we regressed WHR on BMI, age, age-squared, sex and principal components 1–5. We then performed rank-based inverse normalization on the resulting residuals. We validated the findings in three steps:

25(1) Directional consistency. We checked for directional consistency between the effect of index SNPs on WHRadjBMI from the main meta-analysis and EXTEND. We performed linear regression of WHRadjBMI on each individual SNP. We ensured all alleles were aligned to the WHRadjBMI-increasing allele in the original meta-analysis. We compared directions between all 346 index SNPs and then split these into novel and known signals to determine the number of novel signals showing consistent directionality.

26(2) Variance explained. We evaluated the proportion of variance explained by including all the index SNPs into a linear regression model and calculated the adjusted R2. We performed these analyses using the lm function in R.

27(3) Polygenic scores. We created a weighted polygenic score based on the 346 index SNPs associated with WHRadjBMI. The weighted PRS was calculated by summing the dosage of the WHRadjBMI-increasing alleles (weighted by the effect size on WHRadjBMI from the meta-analysis). We then performed linear regression to test the association between WHRadjBMI and the PRS in our independent dataset.

28We sought to determine how likely the 5% of individuals carrying the most WHRadjBMI-increasing alleles were to meet the World Health Organization (WHO) WHR threshold used to diagnose metabolic syndrome (along with lipids and T2D status) compared to the 5% carrying the least. We used the WHR reference levels of >0.9 in men and >0.85 in women to define cases and WHR <0.9 in men and <0.85 in women to define controls. We excluded all individuals with missing data, leaving a sample size of 7513. We took 5% of individuals (7513 × 0.05 = 376) from the two ends of weighted PRS and coded them as 1 or 2. We tested for the likelihood of the top 5% meeting the WHR threshold to diagnose metabolic syndrome (WHO criteria) compared to the bottom 5% using a binomial logistic regression model adjusting for age, age-squared, sex and principal components 1–5.

29Given that we had conditioned WHR on the BMI phenotype for analysis (and BMI and WHR are correlated; r = 0.433 in the UK Biobank data; Supplementary Material, Fig. 15), we tested all index signals found in the WHRadjBMI analysis for evidence of collider bias. To do this, we ran meta-analyses of BMI and WHR using the UK Biobank samples and pre-existing summary-level data from GIANT (Supplementary Material, Methods). We performed these meta-analyses using identical methods to the meta-analysis of WHRadjBMI.

30Then, for each index SNP from the WHRadjBMI meta-analyses (combined as well as sex-specific), we extracted the association results from the BMI and WHR meta-analyses (Supplementary Material, Fig. 4). WHRadjBMI-associated SNPs with a stronger association for BMI than WHR show evidence of collider bias or pleiotropy. We additionally looked at the effect size and direction of effect in BMI and WHR, but whether the effects are from collider bias or pleiotropy cannot be determined from this data.

31We estimated correlation between WHRadjBMI in females and in males using bivariate LD score regression analysis.

32We performed sex-specific GWAS in UK Biobank and meta-analyzed the results with publicly available sex-specific data from the GIANT consortium. We identified the primary and secondary signals from these meta-analyses using methods identical to those performed in the combined analysis. We tested each primary and secondary signal for a sex-dimorphic effect by estimating the t-statisticwhere se is the standard error and r is the genome-wide Spearman rank correlation coefficient between SNP effects in females and males. We estimated the t-statistic and the resulting so-called Pdiff (P-value from a t-distribution) as implemented in the EasyStrata software.

33We tested a total of 2162 different index SNPs for sex dimorphism; we tested all of the secondary signals as well, but these signals are by definition in LD with the index SNPs (and therefore not independent). Given that we tested for sex dimorphism at index SNPs in not only WHRadjBMI but WHR and BMI as well, we performed a test at 1502 distinct genomic loci. Therefore, we set significance for sex dimorphism at a Bonferroni-corrected P = 0.05/1502 = 3.3 × 10−5.

34SNPs were determined to have a stronger effect in women if they fell into one of the following categories (abs, absolute value): betafemales ≤ 0 and betamales ≤ 0 and abs(betafemales) > abs(betamales)betafemales ≥ 0 and betamales ≥ 0 and abs(betafemales) > abs(betamales)betafemales ≤ 0 and betamales ≥ 0 and Pfemales < Pmales and abs(betafemales) > abs(betamales), orbetafemales ≥ 0 and betamales ≤ 0 and Pfemales < Pmales and abs(betafemales) > abs(betamales)

35betafemales ≤ 0 and betamales ≤ 0 and abs(betafemales) > abs(betamales)

36betafemales ≥ 0 and betamales ≥ 0 and abs(betafemales) > abs(betamales)

37betafemales ≤ 0 and betamales ≥ 0 and Pfemales < Pmales and abs(betafemales) > abs(betamales), or

38betafemales ≥ 0 and betamales ≤ 0 and Pfemales < Pmales and abs(betafemales) > abs(betamales)

39We implemented all heritability calculations in BOLT-LMM. We used the same GRM to estimate SNP-based heritability as we did to run our GWAS (see Genome-wide association analyses). This GRM included 790 000 SNPs. Heritability was estimated using only the UK Biobank samples, for which we had individual-level data; these estimates are likely more accurate than those resulting from only summary-level data. We used Restricted Maximum Likelihood Estimation, implemented as REML in BOLT.

40To test the impact of including lower-frequency SNPs in the heritability estimates, we constructed an additional GRM identically as we had for association testing but including no MAF threshold. This GRM included ∼1.7M SNPs. Heritability analyses were calculated identically using this GRM and REML in BOLT.

41To calculate whether heritability estimates in men and women were sex-dimorphic, we used the following equation to generate a z-score:

42We then converted the z-scores to P-values using the following formula in the statistical programming language and software suite R (version 3.4):

43Similarly to Shungin et al., we carried out analysis on the 346 index SNPs and their association with BF% and WHR. We obtained association statistics for the 346 SNPs on BF% and WHR from a GWAS of 443 001 unrelated European-ancestry UK Biobank individuals. We aligned all results to the WHR-increasing allele and used a Bonferroni-corrected P-value (0.05/346 = 1.44 × 10−4) to determine if an SNP was associated with BF% (Fig. 2). To determine whether sex-specific WHRadjBMI index SNPs have an adiposity phenotype, we took the 97 (female-specific) and 8 (male-specific) SNPs and independently compared their effects on WHRadjBMI and BF% in men and women. To identify which sex-dimorphic SNPs were strongly associated with BF% in men and women separately, we used a Bonferroni-corrected P-value of 0.05/105 (4.8 × 10−4) (Supplementary Material, Fig. 7 and Supplementary Material, Table 9). We obtained Pearson’s r correlations using the cor function in R for each comparison.

44Recently, Chu et al. performed a GWAS of subcutaneous and ectopic fat depots, as measured by CT and MRI, in a multi-ancestry sample. Since the meta-analysis results are publicly available (https://grasp.nhlbi.nih.gov/FullResults.aspx and Supplementary Material, Information for further details), we took the index SNPs from our WHRadjBMI meta-analyses (combined sample as well as sex-specific), checked for allele consistency, aligned effects to the reference allele and tested for associations with the imaging-based measures of subcutaneous and ectopic fat. We repeated these analyses in men and women separately. The depots investigated in the imaging-based GWAS were pericardial tissue (PAT), PAT adjusted for height and weight, subcutaneous adipose tissue (SAT), SAT Hounsfield units as measured by MRI, visceral adipose tissue (VAT), VAT Hounsfield units, ratio of VAT to SAT and VAT adjusted for BMI.

45We calculated Pearson’s r correlations between z-scores in WHRadjBMI (calculated by dividing the SNP beta by the standard error) and SNP z-scores reported in Chu et al.. We evaluated significance of the correlation by performing a t-test (implemented as cor.test in R). Correlations were considered significant if P-value < 0.05/3 sample groups/9 phenotypes = 1.9 × 10−3.