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

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¹¹.

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)

Instrument characteristics and assumption checks
From the telomere-length GWAS, 129 genome-wide significant, LD-independent instruments (p < 5×10-8; r2 < 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 allele, yielding a dataset suitable for Wald-ratio construction and pooled estimation.

Primary causal estimates
All estimators were directionally consistent, indicating that longer genetically proxied telomere length is associated with a lower risk of thyrotoxicosis. The primary IVW analysis yielded OR = 0.689 with 95% CI 0.564-0.842 and p < 1.0×10-3. The maximum-likelihood estimator produced OR = 0.689 with 95% CI 0.576-0.823 and p < 1.0×10-3. MR-PRESSO (raw) gave OR = 0.689 with 95% CI 0.561-0.847 and p < 1.0×10-3. MR-RAPS estimated OR = 0.687 with 95% CI 0.574-0.821 and p < 1.0×10-3. For completeness, MR-Egger regression returned OR = 0.671 with 95% CI 0.471-0.955 and p = 2.9×10-2, which is concordant in direction with the IVW estimate. Together, these convergent results support a protective role of longer telomere length against thyrotoxicosis.

Heterogeneity, pleiotropy, and robustness
Between-variant heterogeneity assessed by Cochran's Q did not indicate material excess heterogeneity, supporting the stability of pooled estimates under fixed-effects assumptions, with similar results under random-effects IVW. There was no evidence of directional pleiotropy by the MR-Egger intercept (two-sided p > 1.0×10-1). MR-PRESSO global testing did not identify distortion of the causal estimate after outlier assessment; outlier-corrected estimates remained consistent with the primary result. Leave-one-out analyses did not identify a single SNP driving the association, indicating that results were not dominated by any individual instrument.

Directionality and reverse MR
The Steiger directionality test supported the exposure-to-outcome direction (TL→ thyrotoxicosis), arguing against reverse causation within the MR framework. When treating thyrotoxicosis as the exposure and TL as the outcome, reverse-direction analyses showed no evidence for a causal effect across MR estimators, and accompanying heterogeneity and pleiotropy diagnostics were unremarkable. These findings further reinforce the interpretation that genetically longer telomeres reduce thyrotoxicosis risk.

Visualizations and reporting
Forest plots summarize method-specific odds ratios with 95% confidence intervals, and SNP-level scatter plots display regression lines corresponding to IVW, MR-Egger, maximum likelihood, MR-PRESSO, and MR-RAPS. Density plots of per-SNP ratio estimates are centered near the overall effect with few outliers, and MR-Egger radial plots align along a negative slope (~-0.774), consistent with a protective association. All p-values are reported in scientific notation. Error bars denote 95% CIs. Axes are explicitly labeled in the statistical figures (e.g., β for TL on the x-axis, β for thyrotoxicosis on the y-axis, both on the log-odds scale). Physical scale bars are not applicable to statistical graphics generated from summary-level data.

Data Availability
This study used only de-identified, publicly available GWAS summary statistics. Telomere length data were obtained from IEU OpenGWAS (dataset ieu-b-4879, European ancestry). Thyrotoxicosis summary statistics were obtained from the FinnGen consortium (release 2021 endpoint finn-b-thyROTOXICOSIS; 4,142 cases and 213,693 controls).

To ensure full reproducibility and to satisfy journal requirements, we provide the raw analysis inputs and outputs as Supplementary Files attached to this article: S1_instruments_TL.csv (instrument list with per-SNP F-statistics), S2_harmonised_TL_vs_thyrotoxicosis.csv (allele-aligned dataset), and the S3 series (primary and reverse MR results and diagnostics: S3_mr_results_primary.csv, S3_mr_heterogeneity.csv, S3_mr_pleiotropy_egger.csv, S3_mr_steiger.csv, S3_reverse_mr_results.csv, S3_reverse_heterogeneity.csv, S3_reverse_pleiotropy_egger.csv, S3_reverse_steiger.csv). The exact analysis code is provided as Code S1 (MR_TL_thyrotoxicosis.R), which reproduces all results and generates the 600-dpi figures. No individual-level data were accessed. All datasets were used in accordance with the usage policies of their respective providers.

Figure 1: Study overview and instrument flow. Schematic of data sources and filtering steps: telomere length (IEU OpenGWAS ieu-b-4879) as exposure and thyrotoxicosis (FinnGen 2021 finn-b-thyROTOXICOSIS) as outcome. Instrument thresholds are p < 5×10-8, LD clumping r2 = 0.001, and kb = 10,000 using a European reference. The diagram indicates harmonization to the same effect allele and exclusion of ambiguous palindromes. Error bars and physical scale bars are not applicable to this schematic. Abbreviations: TL, telomere length; LD, linkage disequilibrium. Please click here to view a larger version of this figure.

Figure 2: Forest plot of causal estimates for TL on thyrotoxicosis. Odds ratios are shown per 1 SD increase in TL with 95% CIs for IVW (fixed and random effects), weighted median, maximum likelihood, MR-PRESSO (raw and outlier-corrected), and MR-RAPS. Error bars represent 95% CIs; p-values are printed in scientific notation. This figure visualizes effect-size consistency across methods. Please click here to view a larger version of this figure.

Figure 3: SNP-level scatter plots across MR estimators. Each point represents a SNP with x-axis β for TL and y-axis β for thyrotoxicosis, both on the log-odds scale. Regression lines correspond to IVW, MR-Egger, maximum likelihood, MR-PRESSO, and MR-RAPS; the slope of each line represents the corresponding causal estimator. Axis labels and units are shown on the plot; p-values are reported in scientific notation. The LD exclusion criterion (r2 = 0.001, 10,000 kb) used to define independent instruments is noted in the caption to aid reproducibility. Please click here to view a larger version of this figure.

Figure 4: Leave-one-out sensitivity analysis. Each point shows the inverse-variance-weighted (IVW) causal estimate of telomere length (TL) on thyrotoxicosis after removing one SNP at a time. The horizontal reference line indicates the overall IVW estimate using all instruments. Consistent estimates across removals suggest that the association is not driven by any single SNP (95% CIs shown). Please click here to view a larger version of this figure.

Figure 5: MR-Egger radial plot. Points denote individual SNPs, and the fitted regression line corresponds to the MR-Egger estimator with intercept. The overall negative slope (~ -0.774) indicates a protective association of longer TL. Outliers, if any, are annotated according to MR-PRESSO; p-values are in scientific notation. Please click here to view a larger version of this figure.

Table 1: Causal estimates and diagnostics across MR methods. For each estimator (IVW fixed/random, maximum likelihood, weighted median, MR-Egger, MR-PRESSO raw and outlier-corrected, MR-RAPS), the table lists the OR per 1 SD TL, 95% CI, p-value (scientific notation), nSNP, Cochran's Q with degrees of freedom and p-value, and the MR-Egger intercept p-value. Definitions of effect scale and units are provided in the table footnote.

Supplementary File 1 (S1). S1_instruments_TL.csv Instruments (SNPs) for telomere length (TL) used in the Mendelian randomization analyses, including SNP identifiers and instrument statistics. Please click here to download this File.

Supplementary File 2 (S2). S2_harmonised_TL_vs_thyrotoxicosis.csv Harmonised exposure-outcome dataset for TL and thyrotoxicosis after allele alignment and harmonization. Please click here to download this File.

Supplementary File 3 (S3 (1-8)). S3_MR_outputs.zip Output tables and diagnostics from the MR analyses (primary and reverse-direction), together with the reproducible R script(s) used to generate the results and figures. Please click here to download this File.

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 integrity and modulating cellular lifespan. Progressive telomere shortening accelerates senescence and cellular dysfunction²¹, whereas thyrotoxicosis is typified by heightened metabolic turnover and proliferative signaling²². Under this high-demand state, longer telomeres may buffer oxidative stress and limit DNA damage, thereby maintaining epithelial and follicular cell integrity and delaying pathological remodeling²³. Similar protective links between longer telomeres and reduced risk in other age-related conditions, including cardiovascular disease and some cancers, reinforce a systemic role for telomere dynamics in disease susceptibility²⁴.

This work proposes and implements a genetics-based causal framework to interrogate the TL-thyrotoxicosis relationship. By leveraging MR estimators with orthogonal sensitivity analyses, it addresses key sources of bias that complicate observational studies-most notably residual confounding and reverse causation²⁵. To our knowledge, no prior MR analysis has specifically tested the causal effect of telomere length on thyrotoxicosis; the present findings therefore add novelty by positioning TL as a putative protective factor in thyroid hormone excess states²⁶. Beyond the focal phenotype, the protocol demonstrates a transparent, scriptable analytic pathway that can be generalized to other endocrine traits where mechanistic uncertainty persists.

Identifying telomere length as a potential protective factor has several translational implications. TL could be explored as a biomarker for risk stratification and early detection, particularly in subgroups with fluctuating thyroid function or borderline biochemical profiles. If replicated, TL-informed risk models might complement conventional clinical predictors and support individualized monitoring intervals. Preventive strategies that stabilize telomere dynamics-such as structured physical activity, antioxidant-rich dietary patterns, sleep regularity, and stress reduction-are biologically plausible avenues, while experimental approaches that target telomerase require careful evaluation given context-dependent effects on proliferation and oncogenic risk.

MR mitigates but cannot abolish all bias. Unbalanced horizontal pleiotropy remains possible even after negative MR-Egger intercept findings²⁷ and outlier inspection by MR-PRESSO; canalization, residual population stratification, and measurement error in the exposure GWAS may also influence estimates. Instruments derived from leukocyte TL may not fully capture tissue-specific telomere dynamics in the thyroid, and effect generalizability is constrained by the predominantly European ancestry of the contributing GWAS²⁸. Complementary study designs could strengthen causal inference, including prospective cohorts with repeated TL assessments alongside thyroid hormone panels, tissue-level TL measurements using qPCR/Flow-FISH in surgical samples, and mechanistic experiments in cellular and animal models that manipulate telomere maintenance pathways²⁹. Multivariable MR that adjusts for correlated behavioral or metabolic traits, and phenome-wide scans to map off-target consequences of TL variation, represent additional analytic extensions.

Priority next steps include replication in ancestrally diverse populations; longitudinal profiling to delineate temporal coupling between TL trajectories and thyroid hyperfunction; and experimental interrogation of telomere maintenance, oxidative stress responses, and apoptosis/proliferation checkpoints in thyrocytes under hormone-rich conditions³⁰. Interventional studies that test whether lifestyle-based telomere stabilization reduces incident thyrotoxicosis or blunts relapse risk after treatment would provide clinically actionable evidence. Given the directionality support from Steiger testing and the lack of signal in reverse MR, future work can also examine whether TL interacts with specific autoimmune or nodular etiologies to shape disease onset and course.

In summary, the data support a protective causal association between longer telomeres and reduced risk of thyrotoxicosis. By integrating telomere biology with a robust genetic-instrumental framework, this study advances understanding of thyroid disease mechanisms and highlights telomere length as both a candidate biomarker and a potential-though still cautious-therapeutic lever. These results motivate replication and mechanistic follow-up aimed at translating telomere-informed strategies into preventive and clinical practice.