The causal role of thyroid hormones in bipolar disorders: a two-sample Mendelian Randomization study

Obtaining GWAS summary statistics for thyroid hormones and subtypes of bipolar disorder

GWAS summary statistics for TSH, FT4, FT3, TT3, FT3/FT4, and TT3/FT4 were obtained from Sterenborg et al.⁷ and were based on 271,040, 119,120, 59,061, 15,829, 51,095, and 15,510 individuals with European ancestry, respectively, who were participants in the ThyroidOmics Consortium; these individuals were above 18 years of age, not using thyroid medications, and did not have any history of thyroid surgeries. In generating these summary statistics, TSH, FT4, FT3, and TT3 were treated as continuous variables following an inverse normal transformation, and FT3/FT4 and TT3/FT4 underwent natural logarithm transformation⁷. Summary statistics for bipolar disorder, as well as Type 1 and 2 bipolar disorders were obtained from separate studies detailed in Mullins et al.⁸ involving a total of 41,917 bipolar disorder cases and 371,549 controls with European ancestry; of these cases, 25,060 individuals diagnosed with Type 1 bipolar disorder and 6,781 individuals diagnosed with Type 2 bipolar disorder were included for generating subtype-specific summary statistics. Diagnoses of bipolar disorder in most of these cases were based on the international consensus criteria for lifetime bipolar disorder (DSM-IV, ICD-9, or ICD-10), which were determined through structured diagnostic interviews, clinician-administered checklists, or medical record reviews⁸. All original studies contributing to these GWAS obtained ethical approval, and participants gave informed consent.

Study design of the two-sample Mendelian Randomization (MR)

For the causal estimates derived from a Mendelian Randomization (MR) study to be considered valid, three key assumptions must be fulfilled: (1) the genetic variants must be strongly correlated with the exposure, (2) the genetic variants should not be linked to any confounders of the exposure-outcome relationship, and (3) the variants must not influence the outcome independently of the exposure. In this study, we employed a two-sample MR approach utilizing GWAS summary statistics to assess the association between each thyroid hormone (exposures) on each bipolar subtype (outcomes), utilizing genetic variants as instrumental variables.

Selection of genetic variants as instruments

Separately for each thyroid hormone, we identified genetic variants associated with a p-value<5×10⁻⁸ and performed linkage disequilibrium clumping using PLINK 1.9⁹with a stringent threshold of r²<0.001 and a 10,000 kb window; European individuals from the 1000 genomes project¹⁰ were utilized as the reference panel for clumping. For thyroid hormones for which we identified very few genetic variants after clumping, we loosened the p-value threshold for clumping to 5×10⁻⁷ and 5×10⁻⁶ to obtain at least 10 initial genetic variants per thyroid hormone. Any genetic variants that were associated with the outcome (p<0.05) were then excluded from MR analyses. To assess whether pleiotropic pathways confounded the MR results, we queried all instruments in the NHGRI-EBI Catalog¹¹ to determine if any instruments were associated with potential confounders between thyroid function and bipolar disorders; variants which were associated with potential confounders (including psychiatric and neurological conditions, as well as properties related to neurobiology such as brain function) were excluded from our analyses. For the remaining variants, we calculated the F-statistics using the formula β²/SE² to assess the strength of the selected genetic instruments, ensuring that all instruments used in our MR analyses had sufficient strength (F>10). We determined the variance explained by a specific SNP using the formula 2β²×ƒ×(1-ƒ), where β represents the SNP effect on a given thyroid hormone metric and ƒ represents the effect allele frequency. We aligned the GWAS summary statistics for both exposure and outcome by matching them based on the effect alleles and excluded palindromic variants with ambiguous allele frequencies from the MR analyses^12,13.

Performing two-sample MR between thyroid hormone metrics and subtypes of bipolar disorder

The inverse-variance weighted (IVW) method with multiplicative random effects¹⁴ was then applied to combine SNP-specific causal estimates for each bipolar disorder using the TwoSampleMR v0.6.7¹⁵ package in R; the IVW method is comparable to a two-stage least squares or allele score analysis when using individual-level data, and is therefore regarded as a conventional approach for MR. Given the interrelatedness of thyroid hormone metrics, we utilized a threshold of p-value<8.33×10⁻³ from the IVW method to identifying causal associations after accounting for multiple testing, corresponding to an alpha of 0.05 divided by 6 tests for the six thyroid hormone metrics. We reported suggestive causal associations using a threshold of p-value<0.05.

Sensitivity analyses and assessment of pleiotropy

In addition to IVW, we utilized the TwoSampleMR package to implement MR-Egger¹⁶, weighted median¹⁷, and weighted mode¹⁸ MR methods to assess directional consistency of MR-IVW estimates. To confirm the robustness of significant findings, we conducted the Cochran’s Q-statistic test for heterogeneity, as well as tests for non-zero MR-Egger intercepts¹² and the MR-Pleiotropy RESidual Sum and Outlier (MR-PRESSO) global test to assess pleiotropy¹⁹. For both causal (p-value<8.33×10⁻³) and suggestive causal (p-value<0.05) relationships we identified, we evaluated the influence of individual or pleiotropic variants by both computing MR estimates separately for each instrument, as well as conducting a leave-one-out analysis, where each SNP was systematically removed one at a time¹⁶ as a sensitivity analysis. We also evaluated whether our study had sufficient power to detect the associations we observed between thyroid hormone metrics and bipolar disease by performing power calculations using the mRnd webtool²⁰. Lastly, we assessed whether there was publication bias in the summary statistics utilized for determining MR instruments by creating funnel plots for the causal associations we identified.

Genetic variants selected as instruments for each thyroid hormone metric

We identified genetic variants for each thyroid hormone metric (FT3, FT4, TSH, TT3, FT3/FT4, and TT3/FT4) that met the assumptions to serve as instruments in MR analyses including a total of 52 variants for FT4, 152 for TSH, 12 for FT3/FT4, and less than 10 variants for either FT3, TT3, or TT3/FT4 at a p-value threshold of 5×10⁻⁸ (Supplementary Table 1). When relaxing the p-value threshold to 5×10⁻⁷, we identified 13 suitable instruments for FT3 and 11 for TT3/FT4, and further relaxing the threshold to 5×10⁻⁶ provided 11 variants that qualified as instruments for TT3 (Supplementary Table 1). In total, we identified 8 genetic variants among all instruments across the six thyroid hormone metrics that were associated with pleiotropic pathways in the NHGRI-EBI GWAS Catalog¹¹ that may confound the MR results, including psychiatric and neurological conditions, as well as traits related to brain function/anatomy (Supplementary Table 2); these 8 variants are not included in the aforementioned numbers of instruments nor in any subsequent MR analyses. F-statistics for all valid instrumental variables were greater than 10, ranging from 20.84 to 2183.15, confirming the strength and reliability of these instruments for MR analyses (Supplementary Table 1). Additional characteristics for the instruments utilized in MR analyses are presented in Supplementary Table 1.

Two-sample Mendelian Randomization demonstrates thyroid hormones causally impact the risk of bipolar disorder subtypes

Utilizing genetic instruments for thyroid hormone metrics, we separately performed two-sample MR for overall bipolar disorder, as well as type 1 and type 2 bipolar disorders (Supplementary Table 4). We identified a causal negative relationship between FT4 hormone levels and the odds of bipolar disorder, overall, using the IVW method (OR: 0.92, 95% CI: 0.86, 0.97, p-value: 4.58×10⁻³) even after adjusting for multiple testing (Table 1). The directionality of this relationship between FT4 hormones and bipolar disorder identified by IVW was highly consistent across MR-Egger, weighted median, and weighted mode MR approaches (Figure 1). We also identified FT4 hormone levels had a suggestive causal negative relationship with Type 1 bipolar disorder (OR: 0.92, 95% CI: 0.86, 0.99, p-value: 3.21×10⁻²), which was also directionally consistent across other MR approaches (Table 1 and Figure 1). However, the MR association between FT4 and Type 2 bipolar disorder was not statistically significant (IVW p-value: 0.29), and directionality of this relationship differed between MR approaches (Supplementary Table 4). Furthermore, we identified two other thyroid hormone metrics, FT3 and FT3/FT4, had suggestive causal associations with Type 1 bipolar disorder yielding IVW ORs of 1.18 (95% CI: 1.03, 1.35, p-value: 1.55×10⁻²) and 1.97 (95% CI: 1.02, 3.82, p-value: 4.46×10⁻²), respectively (Table 1); these associations were directionally consistent across MR methods (Table 1 and Figure 1). While the associations between FT3/FT4 and Type 2 bipolar disorder were not statistically significant (p>0.05), the directionality of the association seemed to be opposite to that between FT3/FT4 and Type 1 bipolar disorder (Supplementary Table 4). MR associations between FT3 and Type 2 bipolar disorder also appeared to be opposite in directionality compared to the association with Type 1 bipolar disorder, but these models showed evidence of directional pleiotropy (MR-Egger intercept p-value: 0.048).

• Download figure
• Open in new tab

Figure 1.

Scatterplots for the effects of instruments on each thyroid hormone metric and bipolar disorder subtype for causal and suggestive associations identified in Mendelian Randomization (MR) analyses. A) Relationship between FT4 and bipolar disorder (overall), B) Relationship between FT4 and Type 1 bipolar disorder, C) Relationship between FT3 and Type 1 bipolar disorder, D) Relationship between FT3/FT4 and Type 1 bipolar disorder.

For the one causal and three suggestive causal relationships we identified between thyroid hormone metrics and bipolar disorder subtypes (Table 1), we rigorously assessed potential biases by conducting several sensitivity analyses. The instruments utilized in these relationships did not show evidence of being in any confounding or pleiotropic pathways when examining genome-wide significant associations between these instruments and diseases/traits cataloged in NHGRI-EBI GWAS Catalog (Supplementary Table 4). We did not find evidence for pleiotropy across instrument effects (Cochran’s Q test p-values>0.05) (Table 1). Next, there was no evidence of directional pleiotropy among these causal associations (MR-Egger intercept p-values>0.05) (Table 1). In addition, we computed the MR associations for each instrument separately and observed that they were generally consistent with the estimates inferred by each MR method collectively across all instruments (Figure 2). Moreover, MR-PRESSO did not detect any outlier genetic variants that needed to be removed (MR-PRESSO global test p-values>0.05) (Table 1). To further corroborate the absence of outlier instruments, we performed a leave-one-out analysis where we computed MR estimates but left each instrument one-by-one; we observed our findings were generally reliable and stable as none of the instruments heavily influenced the causal relationships we observed when excluded (Figure 3).

• Download figure
• Open in new tab

Figure 2.

Forest plots of the Mendelian Randomization (MR) estimates for each instrument individually compared to MR approaches using all instruments collectively. A) Relationship between FT4 and bipolar disorder (overall), B) Relationship between FT4 and Type 1 bipolar disorder, C) Relationship between FT3 and Type 1 bipolar disorder, D) Relationship between FT3/FT4 and Type 1 bipolar disorder.

• Download figure
• Open in new tab

Figure 3.

Forest plots of the inverse-variance weighted Mendelian Randomization (MR) method estimates utilizing the leave-one-out approach where instruments were excluded one-by-one. A) Relationship between FT4 and bipolar disorder (overall), B) Relationship between FT4 and Type 1 bipolar disorder, C) Relationship between FT3 and Type 1 bipolar disorder, D) Relationship between FT3/FT4 and Type 1 bipolar disorder.

Additionally, we conducted power calculations for sample sizes of 413,466 (41,917 cases and 371,549 controls) for bipolar disorder overall and 396,609 (25,060 cases and 371,549 controls) for Type 1 bipolar disorder, with an alpha of 0.05. The sum of variance in FT4 explained by instruments in the MR with bipolar disorder overall was 3.84%, yielding a power of 0.87 to detect the MR association we observed. For Type 1 bipolar disorder, the sum of variance in FT4, FT3, and FT3/FT4 explained by instruments was 3.20%, 0.96%, and 0.07%, respectively, yielding power values of 0.60, 0.76, and 0.97, respectively, for detecting the MR associations we observed. Funnel plots for each of the causal relationships we identified also suggested that there were no indications of publication bias in this study (Supplementary Figure 1).