Association Between Telomere Length and Thyrotoxicosis: Insights from a Two-sample Mendelian Randomization Study

Summary

This protocol describes a two-sample Mendelian randomization pipeline assessing whether telomere length causally affects thyrotoxicosis risk using public genetic summary data. It covers instrument selection, harmonization, primary estimation, sensitivity analyses, and reproducible R code with figure-ready outputs to support transparent reporting.

Abstract

Thyrotoxicosis is an endocrine disorder characterized by excess thyroid hormones, yet its etiologic links to systemic aging biology remain incompletely defined. Telomere length (TL) reflects cellular senescence and genome stability and has been implicated in multiple complex diseases. We conducted a two-sample Mendelian randomization (MR) study to evaluate the causal effect of genetically predicted TL on the risk of thyrotoxicosis. Genetic instruments for TL were derived from a large genome-wide association study (GWAS) of European ancestry (n > 470,000). Thyrotoxicosis summary statistics were obtained from the latest FinnGen release (≈4,000 cases and >210,000 controls). Primary inverse-variance-weighted analyses indicated that longer genetically proxied TL is associated with a lower risk of thyrotoxicosis, and the direction and magnitude of the effect were consistent across complementary estimators (MR-Egger, weighted median/maximum likelihood, MR-PRESSO, and MR-RAPS). Sensitivity analyses showed no evidence of directional pleiotropy, and Cochran's Q was used to assess heterogeneity. A Steiger directionality test supported the causal flow from TL to thyrotoxicosis.

To our knowledge, this work is among the first MR analyses to assess the causal relationship between overall thyrotoxicosis risk and TL using contemporary GWAS resources, extending prior evidence focused on hyperthyroidism-related phenotypes. These findings suggest that cellular aging processes indexed by TL may contribute to thyrotoxicosis susceptibility and motivate future longitudinal and mechanistic studies on telomere biology in thyroid dysfunction.

Introduction

Thyrotoxicosis is a common endocrine condition characterized by excess circulating thyroid hormones. Typical manifestations include weight loss, tachycardia, irritability, anxiety, and tremor, with severe cases predisposing to arrhythmias, heart failure, and bone loss¹. Although clinical features and therapeutic strategies are relatively well established, the biological mechanisms that shape inter-individual susceptibility and progression remain incompletely defined. Etiologies range from autoimmune hyperfunction (e.g., Graves' disease) to toxic nodules, thyroiditis, exogenous hormone exposure, drug reactions, and tumor-related causes². Across these entities, immune activation, oxidative stress, and altered mitochondrial biogenesis recur as plausible mechanistic themes linking hormonal excess to tissue damage and systemic complications.

Telomere length (TL)-the tandem repeat DNA-protein structure that caps chromosome ends-integrates cellular aging, replicative history, and genome stability. TL shortens with cell division and is accelerated by oxidative stress and chronic inflammation, processes that are also prominent in thyroid dysfunction^3,4. Shorter TL has been associated with cardiometabolic disorders, certain cancers, and neurodegenerative conditions⁵, and emerging data suggest links to endocrine phenotypes such as obesity and type 2 diabetes, in which hormonal and immune-metabolic axes intersect with replicative stress⁶. Given that thyroid hormones regulate basal metabolic rate, mitochondrial function, and proliferative signaling, TL dynamics may be particularly relevant in thyrotoxic states, where systemic hypermetabolism and immune dysregulation could both influence and be influenced by telomere biology^7,8.

However, whether TL plays a causal role in thyrotoxicosis has remained uncertain. Conventional observational studies are vulnerable to confounding (e.g., lifestyle, comorbidities) and reverse causation (disease affecting TL). To address these limitations, we use two-sample Mendelian randomization (MR), which leverages germline genetic variants robustly associated with TL as unconfounded instruments to test its effect on disease risk⁹. Genome-wide association studies now provide sufficiently strong TL instruments and large, well-phenotyped outcome datasets to enable a rigorous causal assessment¹⁰. By integrating genetic epidemiology with telomere biology, this study evaluates whether genetically proxied TL influences the risk of thyrotoxicosis, thereby clarifying mechanistic pathways and exploring TL as a candidate biomarker for risk stratification and future intervention development.

To our knowledge, this analysis is among the first MR studies to interrogate the overall risk of thyrotoxicosis in relation to TL using contemporary large-scale GWAS resources, extending prior work that focused primarily on hyperthyroidism-related phenotypes or autoimmune thyroid disease¹¹.

Protocol

This study applies a two-sample Mendelian randomization (MR) design to evaluate whether genetically proxied telomere length (TL) causally influences the risk of thyrotoxicosis. Only de-identified, publicly available GWAS summary statistics were used, and no individual-level data were accessed. The Institutional Review Board of Qingdao Municipal Hospital determined that the analysis is exempt from further review because it relies exclusively on public summary data. All contributing genome-wide association studies obtained informed consent and ethical approvals as part of their original protocols. The analysis is structured to minimize confounding and reverse causation by using germline variants associated with TL as unconfounded instruments and by implementing sensitivity procedures to interrogate pleiotropy, heterogeneity, and the direction of effect.

Data sources
The exposure dataset for TL was obtained from IEU OpenGWAS under the identifier ieu-b-4879, comprising approximately 472,174 participants of European ancestry. The outcome dataset for thyrotoxicosis was taken from the FinnGen consortium, using the 2021 endpoint finn-b-thyROTOXICOSIS with 4,142 cases and 213,693 controls. These resources provide effective estimates and standard errors required to construct SNP-level instruments, harmonize alleles, and estimate causal effects with established MR estimators.

Software and computing environment
All analyses were conducted in R (version 4.3.1). The principal analytical package was TwoSampleMR (version 0.5.7), complemented by ieugwasr for programmatic access to GWAS resources, MRPRESSO for outlier detection and distortion testing, mr.raps for robust estimation under weak instruments and idiosyncratic pleiotropy, RadialMR for radial visualization, and general-purpose packages including tidyverse and data.table. Session information and package versions are written to a file to ensure strict reproducibility.

Core assumptions of Mendelian randomization
The MR framework presupposes that genetic instruments are strongly associated with TL, are independent of factors that confound the relationship between TL and thyrotoxicosis, and influence thyrotoxicosis only through TL rather than through alternative pathways. The analysis plan operationalizes these assumptions by quantifying per-variant strength, by testing for unbalanced horizontal pleiotropy using intercept-based methods and outlier screening, by evaluating between-variant heterogeneity, and by confirming the TL-to-thyrotoxicosis direction of effect with a formal directionality test¹².

Instrument selection and quality control
Genetic instruments were selected from the TL GWAS at a genome-wide significance threshold of p < 5×10^-8. To ensure independence, linkage disequilibrium was addressed by clumping at r² = 0.001 within a 10,000 kb window using European reference data from the 1000 Genomes Project; where variants were correlated, the variant with the smaller p value for association with TL was retained. Instrument strength was summarized for each variant using the statistic F_i = \beta_{E,i}^{2}/\mathrm{SE}(\beta_{E,i})^{2} and judged against the conventional F > 10 criterion^13,14. Where appropriate, the Sanderson-Windmeijer approximation was referenced to describe aggregate instrument strength across multiple variants in a sample of size N, with K instruments explaining exposure variance R^{2}: F_{\mathrm{SW}} = \{[R^{2}/(1-R^{2})]\,(N-K-1)\}/K ¹⁵. Ambiguous palindromic variants with allele frequencies near 0.5 were excluded unless allele frequency information permitted unambiguous alignment.

Harmonization
For each instrument, the corresponding association with thyrotoxicosis was extracted from FinnGen and aligned so that effect sizes represent the same effect allele across exposure and outcome. Harmonization removed allele mismatches, corrected strand issues, and excluded palindromic SNPs with unresolved ambiguity, producing a dataset suitable for valid Wald ratio construction.

Primary causal estimation
The primary analysis used the inverse-variance-weighted (IVW) estimator to meta-analyze SNP-specific Wald ratios into a pooled causal effect¹⁶. Estimates are reported as odds ratios per one standard-deviation increase in TL, with 95% confidence intervals derived under fixed- and random-effects models. Model choice was guided by heterogeneity diagnostics, and both specifications are presented to facilitate robust interpretation.

Sensitivity analyses and robustness checks
Robustness was assessed using MR-Egger regression with an intercept term to test for unbalanced directional pleiotropy, maximum likelihood estimation to improve efficiency under homogeneity while mitigating measurement error, MR-PRESSO to perform a global outlier test and to estimate distortion and outlier-corrected effects, and MR-RAPS to provide estimates resilient to weak instruments and idiosyncratic pleiotropy¹⁷. Between-variant heterogeneity was quantified using Cochran's Q statistic under the IVW framework¹⁸. Direction of effect was examined with the Steiger test, which compares the proportion of variance explained in exposure and outcome to determine whether the data are more compatible with TL causing thyrotoxicosis rather than the reverse. The MR-Egger intercept was used as a formal test for directional pleiotropy. Leave-one-out analyses were inspected to ensure that the overall association was not driven by any single variant.

Reverse Mendelian Randomization
To probe potential reverse causation, the analytic pipeline was repeated with thyrotoxicosis as the exposure and TL as the outcome. The same instrument selection criteria, harmonization procedures, IVW primary estimator, and sensitivity analyses were applied so that conclusions about directionality are made within an identical causal framework.

Multiple testing, power, and reporting
The primary hypothesis test pertains to the IVW estimator for the effect of TL on thyrotoxicosis. Sensitivity estimators and diagnostic tests are interpreted as supportive evidence; p values are reported in scientific notation for clarity, and conclusions emphasize consistency across methods rather than isolated significance thresholds. Instrument strength summaries and the proportion of variance explained inform approximate power under standard noncentrality formulations, acknowledging that power depends on sample size, instrument strength, and true effect magnitude.

Computational reproducibility and exact commands
Reproducibility is ensured by providing the full sequence of R commands that recreate instrument selection, outcome extraction, harmonization, primary and sensitivity analyses, reverse MR, diagnostic outputs, and export of analysis-ready files. The script writes stable, human-readable CSV files corresponding to the instrument list, the harmonized dataset, and the summary of MR estimates and diagnostics.

# R 4.3.1; TwoSampleMR 0.5.7
# Optional installation:
# install.packages(c("TwoSampleMR","ieugwasr","MRPRESSO","mr.raps",
"RadialMR","tidyverse","data.table"))

library(TwoSampleMR)
library(ieugwasr)
library(MRPRESSO)
library(mr.raps)
library(RadialMR)
library(tidyverse)
library(data.table)

# Exposure: telomere length (IEU OpenGWAS)
exposure_id <- "ieu-b-4879"

# Outcome: FinnGen 2021 thyrotoxicosis endpoint used for the reported analyses
outcome_id <- "finn-b-thyROTOXICOSIS"

# Instrument selection with genome-wide threshold and stringent LD clumping
exp <- extract_instruments(outcomes = exposure_id, p1 = 5e-8, clump = TRUE, r2 = 0.001, kb = 10000)
exp$F_stat <- (exp$beta.exposure^2) / (exp$se.exposure^2)
fwrite(exp, "S1_instruments_TL.csv")

# Outcome extraction and harmonization
out <- extract_outcome_data(snps = exp$SNP, outcomes = outcome_id)
dat <- harmonise_data(exp, out, action = 2)
fwrite(dat, "S2_harmonised_TL_vs_thyrotoxicosis.csv")

# Primary MR and sensitivity estimators
res <- mr(dat, method_list = c("mr_ivw","mr_ivw_fe","mr_egger_regression",
"mr_weighted_median","mr_raps","mr_maxlik"))

het <- mr_heterogeneity(dat) # Cochran's Q
pleio <- mr_pleiotropy_test(dat) # Egger intercept
steiger <- directionality_test(dat) # Steiger directionality

# MR-PRESSO global and outlier-corrected estimates
mrpresso <- mr_presso(BetaOutcome = "beta.outcome", BetaExposure = "beta.exposure",
SdOutcome = "se.outcome", SdExposure = "se.exposure",
OUTLIERtest = TRUE, DISTORTIONtest = TRUE,
data = dat, NbDistribution = 1000, SignifThreshold = 0.05) #

# Reverse MR: thyrotoxicosis (exposure) -> TL (outcome)
rev_exp <- extract_instruments(outcomes = outcome_id, p1 = 5e-8, clump = TRUE, r2 = 0.001, kb = 10000)
rev_out <- extract_outcome_data(snps = rev_exp$SNP, outcomes = exposure_id)
rev_dat <- harmonise_data(rev_exp, rev_out, action = 2)

rev_res <- mr(rev_dat, method_list = c("mr_ivw","mr_ivw_fe","mr_egger_regression",

"mr_weighted_median","mr_raps","mr_maxlik"))
rev_het <- mr_heterogeneity(rev_dat)
rev_pleio <- mr_pleiotropy_test(rev_dat)
rev_steiger <- directionality_test(rev_dat)

# Exports for archiving and figure/table generation
write.csv(bind_rows(res), "S3_mr_results_primary.csv", row.names = FALSE)
write.csv(het, "S3_mr_heterogeneity.csv", row.names = FALSE)
write.csv(pleio, "S3_mr_pleiotropy_egger.csv",row.names = FALSE)
write.csv(steiger, "S3_mr_steiger.csv", row.names = FALSE)
write.csv(bind_rows(rev_res), "S3_reverse_mr_results.csv", row.names = FALSE)
write.csv(rev_het, "S3_reverse_heterogeneity.csv", row.names = FALSE)
write.csv(rev_pleio, "S3_reverse_pleiotropy_egger.csv", row.names = FALSE)
write.csv(rev_steiger, "S3_reverse_steiger.csv", row.names = FALSE)

Results

Instrument characteristics and assumption checks
From the telomere-length GWAS, 129 genome-wide significant, LD-independent instruments (p < 5×10^-8; r² < 0.001; 10,000 kb window) were retained. Per-SNP F-statistics exceeded the conventional threshold for the large majority of variants, indicating limited risk of weak-instrument bias. Harmonization removed ambiguous palindromic variants and aligned alleles so that exposure and outcome effects referred to the same effect all...

Discussion

This Mendelian randomization (MR) study provides evidence that longer genetically predicted telomeres reduce the risk of thyrotoxicosis¹⁹. Directionally concordant estimates across IVW, MR-Egger, maximum likelihood, MR-PRESSO, and MR-RAPS indicate a robust protective association, and directionality checks support a TL → thyrotoxicosis pathway rather than the reverse²⁰. The biological plausibility is consistent with the role of telomeres in preserving chromosomal integ...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
1000 Genomes Project	Reference data for linkage disequilibrium (LD) clumping.	1000 Genomes Project	-
Cochran’s Q statistic	Statistical method used for heterogeneity testing in MR.	N/A	-
FinnGen dataset	GWAS summary statistics for thyrotoxicosis (cases and controls).	FinnGen Consortium	Release 2021, finn-b-thyROTOXICOSIS
IEU OpenGWAS dataset	GWAS summary statistics for telomere length.	IEU OpenGWAS	Dataset ID: ieu-b-4879
ieugwasr package	R package used for programmatic access to GWAS data.	CRAN	Version 0.1.5
mr.raps package	R package for robust estimation in Mendelian randomization under weak instruments.	CRAN	Version 0.1.2
MR-Egger regression	Sensitivity analysis method to test for directional pleiotropy.	N/A	-
MR-PRESSO outlier test	Statistical test used to identify and correct outliers in MR analyses.	N/A	-
MRPRESSO package	R package used for outlier detection and distortion testing in MR analysis.	CRAN	Version 1.0
R software	Statistical computing and graphics environment used for all analyses.	R Foundation for Statistical Computing	Version 4.3.1
RadialMR package	R package for radial visualization of MR results.	CRAN	Version 1.1.2
TwoSampleMR package	R package used for Mendelian randomization analyses.	CRAN	Version 0.5.7