We conducted a two-sample, bi-directional MR. Genetic associations for all traits were obtained from publicly available summary statistics of genome-wide association studies (GWAS)(Table 1). The data were restricted to European populations. To avoid potential weak instrument bias due to sample overlap [15], we sought GWAS with non-overlapping samples whenever possible. This was not possible for the analyses between age-related hearing impairment (ARHI) and body fat percentage (BFP), where we used the GWAS of UK Biobank for these two phenotypes, which entailed a maximum 70% participant overlap. Since both BMI and ARHI genetic associations have a strong signal within the FTO gene region, we conducted colocalisation analysis before the MR analysis to examine whether this signal is driven by a shared causal variant, or whether the association is due to confounding by linkage disequilibrium (LD).
Table 1
Sources of secondary summary statistics used to derive instrumental variables for adult body constitution and hearing loss
Phenotypes | Source | Year of publication | Codes1, references2 or ICD codes3 | Sample size | Cases (n) | Controls (n) | Cases (%) |
Body constitution-related phenotypes | | | | | | | |
Body mass index | UK Biobank | 2018 | 21001 | 359,983 | NA | NA | NA |
Body mass index | GIANT | 2015 | Locke et al. [16] | 322,154 | NA | NA | NA |
Waist circumference | UK Biobank | 2018 | 48 | 360,564 | NA | NA | NA |
Waist circumference | GIANT | 2015 | Shungin et al. [17] | 210,088 | NA | NA | NA |
Body fat percentage | UK Biobank | 2018 | 23099 | 354,628 | NA | NA | NA |
Hearing loss-related phenotypes | | | | | | | |
Sensorineural hearing loss | FinnGen | Release 6, 2022 | ICD 10: H90.3-5; ICD9: 3891 | 252,719 | 19,313 | 233,406 | 7.64 |
Noise-induced hearing loss | Finngen | Release 6, 2022 | ICD 10: H83.3; ICD 9: 3881 | 249,936 | 655 | 249,281 | 0.26 |
Age-related hearing impairment | UK Biobank | Release 6, 2022 | Wells et al. [19] | 250,389 | 87,056 | 163,333 | 34.77 |
Abbreviation: GIANT, Genetic Investigation of ANthropometric Traits; ICD, International Classification of Disease |
Data sources for body constitution-related phenotypes
The adult body constitution was measured by BMI (kg/m2), WC (cm) and BFP (%). We used sex-combined joint GWAS and Metabochip meta-analysis from the Genetic Investigation of ANthropometric Traits (GIANT) for BMI and WC involving 114 studies and 322,154 participants (mean age: 55.52, 46.02% men, mean BMI: 27.40 kg/m2) reported by Locke [16], and 101 studies and 210,088 participants (mean age: 55.25, 45.39% men, mean WC: 93.07 cm) reported by Shungin [17], respectively. Both GWAS were used in the MR analyses with ARHI from UK Biobank.
Apart from GIANT, the data of BMI (n = 359,983) and WC (n = 360,564) were additionally obtained from UK Biobank [18]. The two traits were assessed by anthropometric measurements. We also retrieved the GWAS of BFP from UK Biobank with 354,628 individuals. All traits used here from UK Biobank are continuous phenotypes that have been measured during the initial assessment visit (2006–2010) and rank normalised. The descriptions of BMI, WC and BFP from UK Biobank (data fields: 21001, 48 and 23099, respectively) are presented on its webpage (https://biobank.ctsu.ox.ac.uk/crystal/search.cgi).
Data sources for hearing loss-related phenotypes
Hearing loss traits included sensorineural hearing loss (SNHL), noise-induced hearing loss (NIHL) and ARHI. GWAS on risk of ARHI (87,056 cases and 163,333 controls, mean age: 60) was from Wells et al. [19] based on a self-reported hearing difficulty (HDiff) phenotype in the UK Biobank. The details of questionnaires and case-control assignments for the HDiff phenotype can be found in their published supplement material.
Summary-level genetic data for risk of SNHL and NIHL were obtained from the FinnGen [20] release 6 including 271,341 participants. The phenotypes were collected from nationwide hospital discharge and cause of death registries using the International Classification of Diseases (ICD) codes for SNHL: H90.3, H90.4, H90.5 (ICD-10) and 3891 (ICD-9) and NIHL: H83.3 (ICD-10) and 3881 (ICD-9) (Table 1). Genetic association estimates for SNHL liability were obtained from 19,313 cases (mean age: 63.37, 51.72% men) and 233,406 controls, and for NIHL liability from 655 cases (mean age: 52.43, 92.82% men) and 249,281 controls.
Genetic instrumental variables
Whenever possible, the instrumental variables (IVs) for the exposures were composed of genome-wide significant (P < 5×10− 8) single-nucleotide polymorphisms (SNPs) that are associated with the corresponding phenotypes in the aforementioned GWAS. This included 12 SNPs for SNHL and none for NIHL. Therefore, to include a sufficient number of IVs for MR sensitivity analyses, we included SNPs with P-value < 5×10− 5 as instruments for these traits [21].
We harmonised the effect alleles in the exposure and outcome datasets based on the guidelines provided by Hartwig [22] and excluded palindromic variants with minor allele frequency (MAF) > 0.4 and imputation info score < 0.9. We selected independent SNPs (LD) R2 ≤ 0.001, within a 10000 kb window using the ‘ld_clump’ function on the mrcieu/ieugwasr package. If a variant was unavailable in the outcome GWAS summary statistics, then proxy SNPs were searched for with a minimum LD R2 = 0.9.
F and R2 statistics were calculated for the individual IVs to quantify instrument strength using the method described by Burgess et al. [15]. The statistical power for continuous exposures was evaluated by calculating the approximate minimum detectable odds ratio (OR) for each exposure at a power of 80%, given the sample size of the exposure, the total variance explained by the instruments and the type 1 error rate of 0.05.
Colocalisation analysis
Variants at the FTO locus are associated with BMI and ARHI but it was unclear whether both phenotypes share the same causal variant. We conducted colocalisation analysis to examine it with the method ‘coloc’ developed by Giambartolomei et al. [23], assuming a maximum of one causal variant per genomic locus. The method outputs posterior probabilities (PP) for five models of (H0) no causal variants, (H1) causal variant for exposure, (H2) causal variant for outcome, (H3) distinct causal variants for exposure and outcome, and (H4) shared causal variant for exposure and outcome. A PP > 0.8 for model H4 was considered evidence of colocalization [23]. We examined the variants within the FTO gene region (chr16: 53,737,875—54,155,853 on hg19) and used the default prior probabilities of a variant associated with BMI only (P1) and ARHI only (P2) at 10− 4, and a variant associated with both traits (P12) at 10–5.
Statistical Analyses
We used the inverse variance weighted (IVW) method as the primary analysis in both forward and reverse MR to calculate the IVW estimator by meta-analysing the SNP-specific Wald estimates using a multiplicative random effects model. The Wald estimate for each SNP is calculated as a ratio between of the coefficients for the SNP-to-outcome and the SNP-to-exposure regressions presented in the summary statistics [24]. The random effects model was chosen to account for potential heterogeneity, evaluated by Cochran's Q (assuming balanced pleiotropy) [25]; Evidence of heterogeneity was defined with a P-value < 0.05 in the Cochran's Q. We used R version 4.2.0 [26] and TwosampleMR package [27] for all analyses.
Sensitivity analyses
There are 3 assumptions to consider before validating MR Wald ratios as a robust inference of causation: the IVs (i) are associated with the exposure, (ii) are independent of the confounders and (iii) have no effects on outcome other than via the exposure. In particular, horizontal pleiotropy, in which the IVs have an effect on other traits outside of the pathway of the exposure of interest and have an impact on the outcome of interest, or when the IVs have a direct effect on the outcome of interest, violates the third assumption [28]. We adopted complementary MR methods that make different assumptions about horizontal pleiotropy to assess the robustness of the main MR findings: MR-Egger [29], weighted median [30], weighted mode [31] and MR-PRESSO (Mendelian Randomisation Pleiotropy RESidual Sum and Outlier). MR-PRESSO was also used to detect outliers that may indicate pleiotropic effects [28]. The following additional analyses were conducted: single SNP analysis and leave-one-out analysis to provide visualisation and interpretation of two-sample MR [32].