Large-scale association analysis identified novel differentiated thyroid carcinoma risk loci by integrating transcriptome and proteome

Thérèse Truong, See Hyun Park, Yazdan Asgari, Pierre-Emmanuel Sugier, Mojgan Karimi, Hauke Thomsen, Asta Försti, Rosella Elisei, Federica Gemignani, Stefano Landi, et al.

<figure><img src="image_1.png" /></figure>

To cite this version:

Thérèse Truong, See Hyun Park, Yazdan Asgari, Pierre-Emmanuel Sugier, Mojgan Karimi, et al. Large-scale association analysis identified novel differentiated thyroid carcinoma risk loci by integrating transcriptome and proteome. 2025. hal-05297682

HAL Id: hal-05297682

https://hal.science/hal-05297682v1

Preprint submitted on 22 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Large-scale association analysis identified novel differentiated thyroid carcinoma risk loci by integrating transcriptome and proteome

Therese Truong

therese.truong@inserm.fr

INSERM https://orcid.org/0000-0002-2943-6786

See Hyun Park

INSERM https://orcid.org/0000-0002-7278-7469

Yazdan Asgari

INSERM

Pierre-Emmanuel Sugier

INSERM

Mojgan Karimi

Université Paris-Saclay, UVSQ, Inserm, CESP, 94807, Villejuif, France

Hauke Thomsen

MSB Medical School Berlin

Asta Försti

German Cancer Research Center https://orcid.org/0000-0002-9857-4728

Rosella Elisei

Cisanello, University Hospital of Pisa, Italy

Federica Gemignani

Department of Biology, University of Pisa, Italy

Stefano Landi

University of Pisa https://orcid.org/0000-0001-8364-6357

Florent de Vathaire

Unit 1018 INSERM, Institut Gustave Roussy - Espace Maurice Tubiana https://orcid.org/0000-0002-8374-9281

Claire Mulot

INSERM Université Paris Descartes

Pierre Laurent-Puig

Centre de Recherche des Cordeliers https://orcid.org/0000-0001-8475-5459

Evgenia Ostroumova

International Agency for Research on Cancer, Lyon, France

Pascal Guénel

INSERM https://orcid.org/0000-0002-8359-518X

Jean-François Deleuze

Commissariat à l'énergie atomique et aux énergies alternatives https://orcid.org/0000-0002-5358-4463

Christina (C) Dahm

Aarhus University https://orcid.org/0000-0003-0481-2893

Federico Canzian

German Cancer Research Center (DKFZ)

Gabriella Iannuzzo

Dipartimento di Medicina Clinica e Chirurgia, Federico II University

Carlotta Sacerdote

Marcela Guevara

Department of Epidemiology, Murcia Regional Health Council-IMIB

Ana Jimenez-Zabala

Ministry of Health of the Basque Government

Karl Smith-Byrne

Cancer Epidemiology Unit, University of Oxford https://orcid.org/0000-0002-1932-7463

James Yarmolinsky

Imperial College of London

Sabina Rinaldi

International Agency for Research on Cancer, World Health Organization
Gianluca Severi

Inserm, Université Paris-Saclay, Université Pari-Sud, Institut Gustave Roussy

Fabienne Lesueur

Inserm https://orcid.org/0000-0001-7404-4549

Article

Keywords:

**Posted Date:** July 18th, 2025

DOI: https://doi.org/10.21203/rs.3.rs-6751995/v1

**License:** ©① This work is licensed under a Creative Commons Attribution 4.0 International License. Read Full License

**Additional Declarations:** There is **NO** Competing Interest.
Abstract

Differentiated thyroid cancer (DTC) is a prevalent malignancy with increasing global incidence, yet its genetic susceptibility remains poorly understood. Although previous genome-wide association studies (GWAS) have identified several susceptibility loci, the genetic, transcriptomic, and proteomic factors influencing DTC risk remain unclear. We conducted a large-scale GWAS of 7,681 cases and 963,550 controls of European ancestry. Transcriptome-wide association studies (TWAS) used the joint tissue imputation across thyroid, pituitary, blood, and hypothalamus tissues. Proteome-wide association studies (PWAS) integrated brain and plasma proteomic data to identify proteins influencing DTC risk. Mendelian randomization (MR) and Bayesian colocalization were conducted to infer causality. GWAS identified 18 novel loci associated with DTC risk, four of which were previously suggested and are now confirmed. TWAS identified 29 significant genes, including five genes (LRRC34, NRG1, HEMGN, PTCSC3, and SMAD3) located within known loci and three novel genes (SAMD4A, RAD51-AS1, and MPHOSPH6) validated as causal through MR and Bayesian colocalization. PWAS identified seven significant proteins, with three (MTHFR, KDEL C2, and SAMD4A) confirmed as causal, further highlighting 15q15.1 as a novel risk locus consistently emerging across all omics layers. This integrated multi-omics approach reveals novel genetic and molecular mechanisms underlying DTC, linking genomic variation to gene expression and protein abundance.

Introduction

Thyroid cancer is the most common endocrine malignancy, with its incidence rising rapidly in developed countries over the last three decades¹. Papillary thyroid carcinoma (PTC) and follicular thyroid carcinoma (FTC) are the most frequent types of differentiated thyroid carcinomas (DTC), comprising about 90% of all thyroid cancer. While the exact risk factors remain unclear, lifestyle, environmental, and genetic factors are strongly implicated.

DTC stands out notably high familial risk, suggesting a strong genetic contribution to its etiology². Previous genome-wide association studies (GWAS) have identified four well-established DTC susceptibility loci at 2q35, 8p12, 9q22.33, and 14q13.3³. A recent large-scale GWAS uncovered an additional five novel risk loci at 1q42.2, 3q26.2, 5p15.33, 5q22.1, 10q24.33, and 15q22.33⁴. However, most variants lie in non-coding regions, making their biological interpretation challenging.

In this study, we conducted a large-scale GWAS meta-analysis including 7,681 DTC cases and 963,550 controls of European ancestry. We integrated transcriptome-wide association studies (TWAS) using eQTL and proteome-wide association studies (PWAS) using pQTL to identify novel DTC risk genes and uncover biological mechanism driving DTC development. This multi-omics approach helps link genetic variants to disease, offering a more comprehensive understanding of DTC etiology.

We also used Mendelian Randomization (MR)⁵ to assess the causal effects of gene expression and protein abundance on DTC risk. Finally, we proposed an updated polygenic risk score (PRS) that integrates novel susceptibility loci to improve risk prediction.

Methods and material

Genome-wide association study meta-analysis of DTC

We conducted a GWAS meta-analysis of DTC using data from four studies with individual-level genotype and phenotype data, and two additional studies providing summary statistics. We performed logistic regression analyses adjusting for sex, age and 10 principal components in the two case-control studies, EPITHYR³ (1,552 cases and 1,954 controls) and ITALIAN⁶ (632 cases and 430 controls), and two nested case-control studies, EPIC⁷ (345 cases and 783 controls) and UK Biobank (518 cases and 358,640 controls) cohorts. Summary statistics for deCODE genetics⁴ (3,001 cases and 287,550 controls) and FinnGen⁸ (1,633 cases and 314,193 controls) were obtained from publicly available sources. Detailed information about the participants, genotyping, quality control, imputation and statistical analyses is summarized in the Supplementary Methods 1. We incorporated all summary statistics by fixed effect inverse-variance weighted meta-analysis using GWAMA software⁹ The heterogeneity of risk estimates was evaluated using Cochran's Q statistic and I². After imputation, 20.8 million variants were analyzed, based on Genome Reference Consortium GRCh38 (hg38).

Transcriptome-wide association study of DTC

We conducted a TWAS to identify genetically predicted gene expression associated with DTC. TWAS integrates GWAS data with cis-eQTL information, enhancing statistical power for identifying gene-disease associations. We used pre-computed gene expression with common cis-eQTL weights from the joint tissue imputation (JTI) model across four thyroid relevant tissues (thyroid, pituitary, whole blood, and brain hypothalamus) from the Genotype-Tissue Expression Project (GTEx, v8)¹⁰. The JTI model improves prediction accuracy by leveraging shared genetic regulation across tissues¹¹. The model uses similarity measures tuned through 5-fold cross-validation. Genes with strong predictive accuracy (r² > 0.1 and p-value < 0.05) from the 5-fold cross-validation were classified as imputable and included in subsequent analyses. Single tissue association tests for the four tissues were conducted using DTC GWAS summary statistics, JTI-derived gene expression weights and a SNP-correlation (LD) matrix from the 1,000 Genome Project through S-PrediXcan¹². Sensitivity analyses across 49 normal tissues accounted for variations in gene expression on different tissue. We also used the thyroid cancer tissue from The Cancer Genome Atlas (TCGA), as cancer significantly alters gene expression in surrounding tissue, providing a more comprehensive understanding. Cis-locus was defined within ± 500kb of gene boundaries. Statistical significance was adjusted for multiple testing using the Bonferroni correction, based on the number of genes tested.

Proteomics-wide association study of DTC
Proteins are highly effective biomarkers and therapeutic targets, as they are the primary functional components of cellular and biological processes and the end products of gene expression¹³. The hypothalamus-pituitary-thyroid axis suggests that proteomic profiles in blood plasma and brain tissue may influence DTC progression. We use proteomic data from two studies: the Religious Orders Study and Memory and Aging Project (ROS/MAP)¹⁴ for brain proteome and the Atherosclerosis Risk in Communities (ARIC)¹⁵ studies for plasma proteome.

The ROS/MAP study provided brain proteomes from 400 postmortem samples of the dorsolateral prefrontal cortex (dPFC). Proteomic analysis was conducted using isobaric tandem mass tag (TMT) peptide labeling and liquid chromatography coupled with mass spectrometry (LC-MS). Samples were randomized into 50 batches to minimize batch effects based on age, sex, post-mortem interval, cognitive diagnosis, and pathologies. Genotypes were derived from whole-genome sequencing or genome-wide genotyping using Illumina OmniQuad Express or Affymetrix GeneChip 6.0 platforms¹⁶. Of 8,356 proteins in 376 European descent subjects, 1,475 proteins were significantly *cis*-heritable and included in our PWAS.

The ARIC study, a cohort of 15,792 participants from four US communities, provided blood plasma proteome data from 7,213 European American. Plasma protein concentrations were measured using SomaLogic's V4 platform with an aptamer (SOMAmer)-based approach. Of 4,877 SOMAmers, 4,657 passed quality control, identifying 4,483 unique proteins encoded by 4,435 autosomal genes. 1,350 proteins were significantly *cis*-heritable and retained for PWAS analysis.

Protein weights were calculated using Functional Summary-based Imputation (FUSION)¹⁷, which estimates the association between brain and blood plasma proteins with DTC, utilizing the LD matrix from the 1,000 Genome Project. The genetic effect of DTC was computed by calculating the linear sum of GWAS Z-scores multiplied by the corresponding protein weights for independent SNPs. We applied false discovery rate (FDR) threshold of 0.05 to account for multiple testing.

## Mendelian randomization

MR uses genetic variants, typically, SNPs, as an instrumental variable (IV) to estimate the causal association between exposure and outcome. This approach minimizes reverse causation and confounding bias, as genetic variants are randomized at conception. For SNPs to be valid, they must satisfy three assumptions: i) strong association with the exposure; ii) no direct association with the outcome except through the exposure (i.e., no directional horizontal pleiotropy); and iii) no association with confounders of the IV-outcome relationship.

We performed LD clumping to obtain independent SNPs with threshold of r² < 0.1 (instead of r² < 0.001 in a standard MR analysis). This threshold was chosen based on a simulation study showing stable coverage rates up to a pruning threshold correlation of 0.16, while avoiding estimates due to multicollinearity emerging around a correlation of 0.36¹⁸. However, regulatory pleiotropy, especially within the same gene locus, can inflate type 1 error rates. To address this, we applied multivariable MR (MVMR) when a SNP was a *cis*-eQTL for multiple genes¹⁹. For each gene identified through TWAS/PWAS, we considered all independent *cis*-eQTLs/pQTLs with an FDR of 0.01 and included other associated genes in the model for each tissue (Supplementary Fig. 1). To avoid multicollinearity, SNPs were clumped to retain only independent variants. The exposure (*cis*-eQTLs/pQTLs) and outcome (DTC GWAS) data were harmonized to ensure the effect alleles aligned. For genes or proteins with a single independent *cis*-QTL, we used the Wald ratio to estimate the causal effect on DTC risk. For those with multiple independent *cis*-QTLs, we applied the inverse-variance weighted (IVW)²⁰ method, combining effect of multiple SNPs in a random-effects model. Finally, MVMR was used to account for pleiotropy when multiple *cis*-QTLs were found across genes or proteins.

## Colocalization Analysis

To identify shared causal variants between significant GWAS loci and eQTL or pQTL signals, we performed a Bayesian colocalization analysis²¹. We focused on variants within ± 500kb of lead variants in each locus and estimated the posterior probability (PP) for five hypotheses: no association (H0), causal variant for GWAS only (H1), causal variant for eQTL or pQTL only (H2), two independent causal variants (H3), and a shared causal variant (H4). A high PP4 (≥ 0.80) indicates that a single variant likely influences both traits within the same genomic region, suggesting a shared causal signal. Additionally, a gene-level colocalization was conducted using FOCUS (Fine-mapping Of CaUsal gene Sets)²² to identify the most likely causal genes associated with DTC risk. This Bayesian method integrates GWAS summary statistics with eQTL data to prioritize genes that mediate the association between genetic variants and DTC risk, providing credible sets of genes with posterior inclusion probabilities (PIP) that quantify the likelihood of each gene being causal for the trait.

## Gene Set Enrichment Analysis

We conducted enrichment analyses of the significant genes identified by GWAS, TWAS or PWAS using Enrichr across multiple resources, including Gene Ontology (GO) term and Kyoto Encyclopedia of Genes and Genomes (KEGG) to explore the biological functions and pathway to DTC²³,²⁴. The statistical significance of each pathway was calculated using the hypergeometric test and Fisher's exact test. Pathway with p-value ≤ 0.05 and at least two overlapping genes were considered significant.

## Polygenic Risk Score

A PRS was constructed to evaluate the cumulative genetic effects of DTC-associated variants. The PRS included GWAS significant SNPs, and its predictive performance was tested on EPIC and EPITHYR studies, with beta estimates derived from analyses excluding these studies to ensure independent validation. Logistic regression was applied to examine the association between the PRS and DTC, adjusting for age, sex, and the first 10 principal components. We also developed a gene/protein-based PRS using TWAS and PWAS to enhance the predictive capability (Supplementary Methods 2). The predictive performance of the PRS model for distinguishing DTC cases from controls was evaluated by calculating the area under the curve (AUC) of the receiver operating characteristic (ROC).
# Results

## GWAS analysis

The overall GWAS Manhattan plot is presented in Fig. 1. There was no indication of genomic inflation from QQ plot ($\lambda_{GC}$ = 1.07), suggesting the absence of significant cryptic population substructure and differences in genotypic variants between the cases and controls (Supplementary Fig. 2). A total of 27 significant loci were identified, including nine previously reported loci (1q42.2, 2q35, 3q26.2, 5q22.1, 8p12, 9q22.33, 10q24.33, 14q13.3, and 15q22.33). Table 1 presents the association of the lead variant with the lowest p-value from the meta-analysis. Details of the association for significant SNPs across all studies are provided in Supplementary Table 1.
<table><caption>Table 1
Association results for lead SNPs at significant loci from the differentiated thyroid cancer GWAS.</caption><thead><tr><th>Locus</th><th>chromosome</th><th>position (GRCh38)</th><th>SNP</th><th>Nearest gene</th><th>EA</th><th>OA</th><th>EAF</th><th>OR (95% CI)</th><th>P value</th></tr></thead><tbody><tr><td colspan="10"><strong>Novel loci</strong></td></tr><tr><td>1p31.3</td><td>1</td><td>61166357</td><td>rs334723<br>P<sub>het</sub> = 0.12; I<sup>2</sup> = 41%</td><td>NFIA</td><td>G</td><td>A</td><td>0.042</td><td>1.36 (1.24, 1.48)</td><td>3.66x10<sup>-11</sup></td></tr><tr><td>1q41</td><td>1</td><td>218424517</td><td>rs1342586<br>P<sub>het</sub> = 0.83; I<sup>2</sup> = 0%</td><td>TGFB2</td><td>C</td><td>T</td><td>0.778</td><td>1.21 (1.16, 1.27)</td><td>1.22x10<sup>-15</sup></td></tr><tr><td>1q43</td><td>1</td><td>243322950</td><td>rs2451668<br>P<sub>het</sub> = 0.82; I<sup>2</sup> = 0%</td><td>SDCCAG8</td><td>T</td><td>C</td><td>0.450</td><td>1.16 (1.11, 1.20)</td><td>1.26x10<sup>-13</sup></td></tr><tr><td>5p15.33</td><td>5</td><td>1280013</td><td>rs7734992<br>P<sub>het</sub> = 0.96; I<sup>2</sup> = 0%</td><td>TERT</td><td>C</td><td>T</td><td>0.440</td><td>1.18 (1.13, 1.22)</td><td>1.39x10<sup>-16</sup></td></tr><tr><td>5q31.1</td><td>5</td><td>134557286</td><td>rs56110108<br>P<sub>het</sub> = 0.37; I<sup>2</sup> = 7%</td><td>JADE2</td><td>T</td><td>A</td><td>0.905</td><td>1.30 (1.19, 1.42)</td><td>6.80x10<sup>-9</sup></td></tr><tr><td>6p21.1</td><td>6</td><td>43935540</td><td>rs1326141<br>P<sub>het</sub> = 0.66; I<sup>2</sup> = 0%</td><td>POLR1C</td><td>A</td><td>G</td><td>0.235</td><td>1.18 (1.13, 1.23)</td><td>3.91x10<sup>-13</sup></td></tr><tr><td>7q31.33</td><td>7</td><td>124757288</td><td>rs2299903<br>P<sub>het</sub> = 0.95; I<sup>2</sup> = 0%</td><td>GPR37</td><td>G</td><td>A</td><td>0.384</td><td>1.12 (1.07, 1.16)</td><td>2.89x10<sup>-08</sup></td></tr><tr><td>8q24.22</td><td>8</td><td>132871660</td><td>rs79676842<br>P<sub>het</sub> = 0.22; I<sup>2</sup> = 29%</td><td>TG</td><td>T</td><td>G</td><td>0.025</td><td>1.42 (1.26, 1.60)</td><td>2.02x10<sup>-08</sup></td></tr><tr><td>9p22.1</td><td>9</td><td>19064131</td><td>rs13287517<br>P<sub>het</sub> = 0.84; I<sup>2</sup> = 0%</td><td>HAUS6</td><td>C</td><td>G</td><td>0.379</td><td>1.14 (1.1, 1.18)</td><td>8.65x10<sup>-11</sup></td></tr><tr><td>10q22.3</td><td>10</td><td>77924587</td><td>rs1650149<br>P<sub>het</sub> = 0.56; I<sup>2</sup> = 0%</td><td>DLG5</td><td>A</td><td>G</td><td>0.708</td><td>1.15 (1.10, 1.20)</td><td>1.85x10<sup>-10</sup></td></tr><tr><td>12q14.3</td><td>12</td><td>65650504</td><td>rs12318900<br>P<sub>het</sub> = 0.76; I<sup>2</sup> = 0%</td><td>PCNPP3</td><td>A</td><td>G</td><td>0.030</td><td>1.37 (1.24, 1.52)</td><td>3.19x10<sup>-9</sup></td></tr><tr><td>15q15.1</td><td>15</td><td>40672837</td><td>rs62019923<br>P<sub>het</sub> = 0.78; I<sup>2</sup> = 0%</td><td>RAD51</td><td>T</td><td>C</td><td>0.843</td><td>1.17 (1.11, 1.23)</td><td>2.31x10<sup>-8</sup></td></tr><tr><td>16q22.2</td><td>16</td><td>72721676</td><td>rs59831429<br>P<sub>het</sub> = 0.72; I<sup>2</sup> = 0%</td><td>ZFHX3-AS1</td><td>T</td><td>C</td><td>0.938</td><td>1.26 (1.16, 1.36)</td><td>3.70x10<sup>-8</sup></td></tr></tbody></table>

Positions are based on Build 38 of the reference genome. Displayed data include the locus, nearest gene, effect allele (EA), the other allele (OA), the effect allele frequency (EAF), and the odds ratio (OR) with the upper and lower 95% confidence intervals (CI). The P value indicates the association between variants and disease, assessed using logistic regression. Meta-analysis results from various studies were combined using a fixed-effect model. P value for heterogeneity (Phet) among the study groups and the heterogeneity statistic (I²), which represents the proportion of variability attributed to heterogeneity between study groups.
<table><thead><tr><th>Locus</th><th>chromosome</th><th>position (GRCh38)</th><th>SNP</th><th>Nearest gene</th><th>EA</th><th>OA</th><th>EAF</th><th>OR (95% CI)</th><th>P value</th></tr></thead><tbody><tr><td colspan="10"><strong>Novel loci</strong></td></tr><tr><td>16q23.2</td><td>16</td><td>79678658</td><td>rs11645076</td><td>MAFTRR</td><td>T</td><td>C</td><td>0.335</td><td>1.18 (1.13, 1.24)</td><td>1.56x10<sup>-11</sup></td></tr><tr><td></td><td></td><td></td><td colspan="5">P<sub>het</sub> &lt;0.01; I<sup>2</sup> = 72%</td><td></td><td></td></tr><tr><td>18p11.32</td><td>18</td><td>799487</td><td>rs143705522</td><td>YES1</td><td>A</td><td>C</td><td>0.988</td><td>2.00 (1.56, 2.57)</td><td>4.22x10<sup>-8</sup></td></tr><tr><td></td><td></td><td></td><td colspan="5">P<sub>het</sub> =1.00; I<sup>2</sup> = 0%</td><td></td><td></td></tr><tr><td>19p13.2</td><td>19</td><td>7223837</td><td>rs4804416</td><td>INSR</td><td>T</td><td>G</td><td>0.586</td><td>1.13 (1.09, 1.18)</td><td>4.19x10<sup>-10</sup></td></tr><tr><td></td><td></td><td></td><td colspan="5">P<sub>het</sub> = 0.28; I<sup>2</sup> = 19%</td><td></td><td></td></tr><tr><td>19p12</td><td>19</td><td>22032639</td><td>rs8105767</td><td>ZNF257</td><td>G</td><td>A</td><td>0.299</td><td>1.15 (1.1, 1.2)</td><td>5.09x10<sup>-11</sup></td></tr><tr><td></td><td></td><td></td><td colspan="5">P<sub>het</sub> = 0.39; I<sup>2</sup> = 5%</td><td></td><td></td></tr><tr><td>22q12.1</td><td>21</td><td>28707610</td><td>rs186430430</td><td>CHEK2</td><td>C</td><td>T</td><td>0.006</td><td>3.36 (2.65, 4.31)</td><td>1.40x10<sup>-22</sup></td></tr><tr><td></td><td></td><td></td><td colspan="5">P<sub>het</sub> = 0.11; I<sup>2</sup> = 55%</td><td></td><td></td></tr><tr><td colspan="10"><strong>Known loci</strong></td></tr><tr><td>1q42.2</td><td>1</td><td>233276132</td><td>rs6697791</td><td>PCNX2</td><td>C</td><td>T</td><td>0.796</td><td>1.23 (1.17, 1.29)</td><td>1.34x10<sup>-16</sup></td></tr><tr><td></td><td></td><td></td><td colspan="5">P<sub>het</sub> = 0.07; I<sup>2</sup> = 49%</td><td></td><td></td></tr><tr><td>2q35</td><td>2</td><td>217431651</td><td>rs57481445</td><td>DIRC3</td><td>G</td><td>A</td><td>0.286</td><td>1.40 (1.34, 1.46)</td><td>4.80x10<sup>-59</sup></td></tr><tr><td></td><td></td><td></td><td colspan="5">P<sub>het</sub> = 0.91; I<sup>2</sup> = 0%</td><td></td><td></td></tr><tr><td>3q26.2</td><td>3</td><td>169793001</td><td>rs9868000</td><td>LRRC34</td><td>G</td><td>A</td><td>0.735</td><td>1.17 (1.12, 1.22)</td><td>2.76x10<sup>-12</sup></td></tr><tr><td></td><td></td><td></td><td colspan="5">P<sub>het</sub> = 0.70; I<sup>2</sup> = 0%</td><td></td><td></td></tr><tr><td>5q22.1</td><td>5</td><td>112150207</td><td>rs73227498</td><td>EPB41L4A</td><td>A</td><td>T</td><td>0.865</td><td>1.32 (1.25, 1.40)</td><td>2.80x10<sup>-20</sup></td></tr><tr><td></td><td></td><td></td><td colspan="5">P<sub>het</sub> = 0.80; I<sup>2</sup> = 0%</td><td></td><td></td></tr><tr><td>8p12</td><td>8</td><td>32549084</td><td>rs4733128</td><td>NRG1</td><td>T</td><td>C</td><td>0.494</td><td>1.29 (1.24, 1.34)</td><td>3.80x10<sup>-38</sup></td></tr><tr><td></td><td></td><td></td><td colspan="5">P<sub>het</sub> = 0.49; I<sup>2</sup> = 0%</td><td></td><td></td></tr><tr><td>9q22.33</td><td>9</td><td>97782586</td><td>rs7847663</td><td>PTCSC2</td><td>C</td><td>T</td><td>0.356</td><td>1.65 (1.59, 1.71)</td><td>2.47x10<sup>-144</sup></td></tr><tr><td></td><td></td><td></td><td colspan="5">P<sub>het</sub> = 0.43; I<sup>2</sup> = 0%</td><td></td><td></td></tr><tr><td>10q24.33</td><td>10</td><td>103916188</td><td>rs9419958</td><td>STN1</td><td>T</td><td>C</td><td>0.131</td><td>1.24 (1.18, 1.31)</td><td>1.88x10<sup>-15</sup></td></tr><tr><td></td><td></td><td></td><td colspan="5">P<sub>het</sub> = 0.37; I<sup>2</sup> = 7%</td><td></td><td></td></tr><tr><td>14q13.3</td><td>14</td><td>36114681</td><td>rs56400346</td><td>PTCSC3</td><td>A</td><td>G</td><td>0.605</td><td>1.32 (1.27, 1.38)</td><td>6.37x10<sup>-42</sup></td></tr></tbody></table>

Positions are based on Build 38 of the reference genome. Displayed data include the locus, nearest gene, effect allele (EA), the other allele (OA), the effect allele frequency (EAF), and the odds ratio (OR) with the upper and lower 95% confidence intervals (CI). The P value indicates the association between variants and disease, assessed using logistic regression. Meta-analysis results from various studies were combined using a fixed-effect model. P value for heterogeneity (Phet) among the study groups and the heterogeneity statistic (I²), which represents the proportion of variability attributed to heterogeneity between study groups.
<table><thead><tr><th>Locus</th><th>chromosome</th><th>position (GRCh38)</th><th>SNP</th><th>Nearest gene</th><th>EA</th><th>OA</th><th>EAF</th><th>OR (95% CI)</th><th>P value</th></tr></thead><tbody><tr><td colspan="10"><strong>Novel loci</strong></td></tr><tr><td></td><td></td><td></td><td>P<sub>het</sub> = 0.22; I<sup>2</sup> = 29%</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>15q22.33</td><td>15</td><td>67150258</td><td>rs17293632</td><td>SMAD3</td><td>T</td><td>C</td><td>0.250</td><td>1.23 (1.17, 1.28)</td><td>2.41x10<sup>-20</sup></td></tr><tr><td></td><td></td><td></td><td>P<sub>het</sub> = 0.87; I<sup>2</sup> = 0%</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td colspan="10">Positions are based on Build 38 of the reference genome. Displayed data include the locus, nearest gene, effect allele (EA), the other allele (OA), the effect allele frequency (EAF), and the odds ratio (OR) with the upper and lower 95% confidence intervals (CI). The P value indicates the association between variants and disease, assessed using logistic regression. Meta-analysis results from various studies were combined using a fixed-effect model. P value for heterogeneity (Phet) among the study groups and the heterogeneity statistic (I<sup>2</sup>), which represents the proportion of variability attributed to heterogeneity between study groups.</td></tr></tbody></table>

## Novel loci

We discovered 18 novel loci significantly associated with DTC risk at genome-wide significance level (p < 5.0x10<sup>-8</sup>), with regional plots shown in supplementary Fig. 3. Among these, 14 had more than one significant SNP within ± 500kb of the lead variant (1p31.3, 1q41, 1q43, 5p15.33, 6p21.1, 7q31.33, 8q24.22, 9p22.1, 10q22.3, 15q15.1, 16q23.2, 19p13.2, 19p12 and 22q12.1). Four loci at 1p31.3 (NFIA, rs334723[G], OR = 1.36; P = 3.66x10<sup>-11</sup>), 5p15.33 (TERT, rs7734992[C], OR = 1.17; P = 1.40x10<sup>-16</sup>), 16q23.2 (MAFTRR, rs11645076[T], OR = 1.18; P = 1.56x10<sup>-11</sup>) and 19p12 (ZNF257, rs8105767[G], OR = 1.15; P = 5.10x10<sup>-11</sup>) have been suggested as susceptibility loci from previous studies<sup>3,4</sup>, but are now confirmed at the genome wide significant level in our analysis (Fig. 2).

## TWAS analysis

A total of 15,814 genes were analyzed across thyroid relevant tissues through TWAS (Supplementary Table 2), and we identified 29 genes across 13 distinct loci significantly associated with DTC risk at the Bonferroni-corrected significant (P < 0.05/15,814 = 3.16x10<sup>-6</sup>), including 26 protein-coding and 3 non-coding RNA genes (Fig. 1). Sixteen genes were located within six previously identified loci (3q26.2, 8p12, 9q22.33, 10q24.33, 14q13.3 and 15q22.33) while 13 genes were located at seven novel loci (1q41, 1q43, 9p22.1, 14q22.2, 15q15.1, 16q23.2 and 19p12) (Table 2). Notably, all these loci, except 14q22.2, were also reported in the GWAS analysis above. Thirteen genes (7 known and 6 novel genes) had predicted expression negatively associated with DTC risk, while 16 genes (9 known and 7 novel) had predicted expression positively associated to DTC risk. These associations were further confirmed by an eQTL analysis, which explores the relationship between genetic variants affecting gene expression and those with DTC risk in corresponding tissues (Supplementary Fig. 4).
<table><caption>Table 2
Results of the 29 genes identified by TWAS followed by Transcriptome-wide Mendelian randomization and Bayesian posterior probability of colocalization</caption><thead><tr><th>Locus</th><th>Gene</th><th>Gene Name</th><th>TWAS</th><th>Tissue</th><th>R<sup>2</sup></th><th>Z-score</th><th>aTWAS P</th><th>Method</th><th>nsnp</th><th>OR (95% CI)</th><th>bMR P</th><th>COLOC</th></tr></thead><tbody><tr><td>1q41</td><td>TGFB2</td><td>Transforming growth factor beta 2</td><td>Thyroid</td><td></td><td>0.145</td><td>-11.25</td><td>2.44x10<sup>-29</sup></td><td>IVW</td><td>4</td><td>0.57 (0.52-0.64)</td><td>1.19x10<sup>-26</sup></td><td>29.70%</td></tr><tr><td rowspan="2">1q43</td><td>SDCCAG8</td><td>SHH signaling and ciliogenesis regulator</td><td>Thyroid</td><td></td><td>0.129</td><td>6.52</td><td>6.83x10<sup>-11</sup></td><td>IVW</td><td>4</td><td>1.34 (1.12-1.61)</td><td>1.81x10<sup>-03</sup></td><td>56%</td></tr><tr><td>CEP170</td><td>Centrosomal protein 170</td><td>Whole Blood</td><td></td><td>0.085</td><td>5.91</td><td>3.39x10<sup>-09</sup></td><td>cMVMR</td><td>5</td><td>2.23(1.15-4.34)</td><td>0.018</td><td>0.23%</td></tr><tr><td rowspan="2">9p22.1</td><td>SAXO1</td><td>Stabilizer of axonemal microtubules 1</td><td>Pituitary</td><td></td><td>0.042</td><td>-5.76</td><td>8.45x10<sup>-09</sup></td><td>h<sub>-</sub></td><td>-</td><td>-</td><td>-</td><td>75%</td></tr><tr><td>HAUS6</td><td>HAUS augmin like complex subunit 6</td><td>Thyroid</td><td></td><td>0.068</td><td>-5.13</td><td>2.85x10<sup>-07</sup></td><td>IVW</td><td>2</td><td>0.78 (0.58-1.04)</td><td>0.096</td><td>0.03%</td></tr><tr><td>14q22.2</td><td>SAMD4A</td><td>Sterile alpha motif domain containing 4A</td><td>Thyroid</td><td></td><td>0.064</td><td>5.90</td><td>3.75x10<sup>-09</sup></td><td>IVW</td><td>2</td><td>1.65 (1.25-2.17)</td><td>4.02x10<sup>-04</sup></td><td>92%</td></tr><tr><td rowspan="5">15q15.1</td><td>RAD51-AS1</td><td>RAD51 antisense RNA 1</td><td>Thyroid</td><td></td><td>0.037</td><td>-5.81</td><td>6.25x10<sup>-09</sup></td><td>dMVMR</td><td>6</td><td>0.41(0.25-0.66)</td><td>2.64x10<sup>-04</sup></td><td>84%</td></tr><tr><td>CHST14</td><td>Carbohydrate sulfotransferase 14</td><td>Thyroid</td><td></td><td>0.021</td><td>-5.69</td><td>1.26x10<sup>-08</sup></td><td>h<sub>-</sub></td><td>-</td><td>-</td><td>-</td><td>1.53%</td></tr><tr><td>RMDN3</td><td>Regulator of microtubule dynamics 3</td><td>Hypothalamus</td><td></td><td>0.023</td><td>5.13</td><td>2.95x10<sup>-07</sup></td><td>h<sub>-</sub></td><td>-</td><td>-</td><td>-</td><td>30.40%</td></tr><tr><td>CCDC32</td><td>Coiled-coil domain containing 32</td><td>Thyroid</td><td></td><td>0.338</td><td>5.37</td><td>7.73x10<sup>-08</sup></td><td>eMVMR</td><td>21</td><td>1.09(0.96-1.23)</td><td>0.206</td><td>72.50%</td></tr><tr><td>MAFTRR</td><td>MAF transcriptional regulator RNA</td><td>Hypothalamus</td><td></td><td>0.427</td><td>-5.67</td><td>1.44x10<sup>-08</sup></td><td>h<sub>-</sub></td><td>-</td><td>-</td><td>-</td><td>1.67%</td></tr><tr><td rowspan="2">16q23.2</td><td>MPHOSPH6</td><td>M-phase phosphoprotein 6</td><td>Pituitary</td><td></td><td>0.287</td><td>4.74</td><td>2.15x10<sup>-06</sup></td><td>IVW</td><td>2</td><td>1.15 (1.08-1.23)</td><td>2.05x10<sup>-05</sup></td><td>87%</td></tr><tr><td>19p12</td><td>ZNF257</td><td>Zinc finger protein 257</td><td>Whole Blood</td><td></td><td>0.173</td><td>5.52</td><td>3.39x10<sup>-08</sup></td><td>fMVMR</td><td>10</td><td>0.93(0.82-1.06)</td><td>0.293</td><td>93.80%</td></tr><tr><td colspan="12">Known loci</td></tr><tr><td rowspan="4">3q26.2</td><td>ACTRT3</td><td>Actin related protein T3</td><td>Pituitary</td><td></td><td>0.034</td><td>7.07</td><td>1.57x10<sup>-12</sup></td><td>h<sub>-</sub></td><td>-</td><td>-</td><td>-</td><td>82.60%</td></tr><tr><td>LRRC34</td><td>Leucine rich repeat containing 34</td><td>Thyroid</td><td></td><td>0.160</td><td>6.44</td><td>1.22x10<sup>-10</sup></td><td>IVW</td><td>4</td><td>1.23 (1.09-1.39)</td><td>6.52x10<sup>-04</sup></td><td>0.00%</td></tr><tr><td>LRRIQ4</td><td>Leucine rich repeats and IQ motif containing 4</td><td>Pituitary</td><td></td><td>0.039</td><td>-5.72</td><td>1.07x10<sup>-08</sup></td><td>h<sub>-</sub></td><td>-</td><td>-</td><td>-</td><td>4.62%</td></tr><tr><td>MECOM</td><td>MDS1 and EVI1 complex locus</td><td>Thyroid</td><td></td><td>0.054</td><td>5.02</td><td>5.05x10<sup>-07</sup></td><td>IVW</td><td>2</td><td>1.52 (0.82-2.82)</td><td>1.79x10<sup>-01</sup></td><td>18.20%</td></tr></tbody></table>

TWAS: Transcriptome-wide association studies; R²: R² of tissue model's correlation to gene's measured transcriptome; MR: Mendelian randomization; IVW: i variance weighted-random effect; MVMR: Multivariable Mendelian randomization; OR: Odds Ratio; CI: Confidence Interval; COLOC: Posterior probability (PP4 Bayesian colocalization; FOCUS: Fine-mapping Of CaUsal gene Sets; PIP: Posterior inclusion probability. aSignificance threshold (0.05/total number of gene 3.16x10⁻⁶); bSignificance threshold (0.05/total number of gene significant after TWAS = 2.78x10⁻³); cadjusted for SDCCAG8 in whole blood tissue; dadjusted for CCDC32, BAHD1 and ZFYVE19 in thyroid tissue; eadjusted for RAD51-AS1, BAHD1 and ZFYVE19 in thyroid tissue; fadjusted for ZNF626, ZNF430, ZNF431, ZNF738, ZNF493, ZNF429, ZNF100 and ZNF43 in whole blood tissue; gadjusted for TBC1D2 in thyroid tissue; hViolation of 1st assumption of Mendelian Randomization; jNot found in the FOCUS eQTL dataset.
<table><thead><tr><th>Locus</th><th>Gene</th><th>Gene Name</th><th>TWAS</th><th>MR</th><th>COLOC</th></tr><tr><th>Novel loci</th><th></th><th></th><th>Tissue</th><th>R<sup>2</sup></th><th>Z-score</th><th>aTWAS P</th><th>Method</th><th>nsnp</th><th>OR (95% CI)</th><th>bMR P</th><th>PP4</th></tr></thead><tbody><tr><td>8p12</td><td>NRG1</td><td>Neuregulin 1</td><td>Thyroid</td><td>0.495</td><td>13.48</td><td>1.94x10<sup>-41</sup></td><td>IVW</td><td>23</td><td>1.38 (1.3-1.46)</td><td>5.36x10<sup>-28</sup></td><td>98.30%</td></tr><tr><td>9q22.33</td><td>PTCSC2</td><td>PTC susceptibility candidate 2</td><td>Thyroid</td><td>0.128</td><td>-10.53</td><td>6.59x10<sup>-26</sup></td><td>h<sub>-</sub></td><td>-</td><td>-</td><td>-</td><td>37.10%</td></tr><tr><td></td><td>TRMO</td><td>tRNA methyltransferase O</td><td>Hypothalamus</td><td>0.291</td><td>-9.87</td><td>5.59x10<sup>-23</sup></td><td>IVW</td><td>6</td><td>0.8 (0.66-0.97)</td><td>2.45x10<sup>-02</sup></td><td>0.00%</td></tr><tr><td></td><td>TRIM14</td><td>Tripartite motif containing 14</td><td>Thyroid</td><td>0.087</td><td>6.78</td><td>1.21x10<sup>-11</sup></td><td>h<sub>-</sub></td><td>-</td><td>-</td><td>-</td><td>1.44%</td></tr><tr><td></td><td>TMOD1</td><td>Tropomodulin 1</td><td>Thyroid</td><td>0.015</td><td>-5.82</td><td>5.77x10<sup>-09</sup></td><td>h<sub>-</sub></td><td>-</td><td>-</td><td>-</td><td>1.57%</td></tr><tr><td></td><td>HEMGN</td><td>Hemogen</td><td>Thyroid</td><td>0.021</td><td>4.83</td><td>1.39x10<sup>-06</sup></td><td>Wald</td><td>1</td><td>1.61 (1.39-1.86)</td><td>3.32x10<sup>-10</sup></td><td>0.13%</td></tr><tr><td></td><td>COR02A</td><td>Coronin 2A</td><td>Thyroid</td><td>0.168</td><td>4.76</td><td>1.94x10<sup>-06</sup></td><td>gMVMR</td><td>8</td><td>1.22(1.05-1.42)</td><td>1.10x10<sup>-02</sup></td><td>0.00%</td></tr><tr><td></td><td>ANP32B</td><td>Acidic nuclear phosphoprotein 32 family member B</td><td>Whole Blood</td><td>0.025</td><td>4.71</td><td>2.53x10<sup>-06</sup></td><td>Wald</td><td>1</td><td>1.89 (1.09-3.3)</td><td>2.42x10<sup>-02</sup></td><td>0.09%</td></tr><tr><td>10q24.33</td><td>STN1</td><td>STN1 subunit of CST complex</td><td>Pituitary</td><td>0.017</td><td>-7.84</td><td>4.55x10<sup>-15</sup></td><td>h<sub>-</sub></td><td>-</td><td>-</td><td>-</td><td>41.90%</td></tr><tr><td>14q13.3</td><td>PTCSC3</td><td>PTC susceptibility candidate 3</td><td>Thyroid</td><td>0.143</td><td>9.79</td><td>1.19x10<sup>-22</sup></td><td>IVW</td><td>5</td><td>1.93 (1.35-2.77)</td><td>3.58x10<sup>-04</sup></td><td>0.00%</td></tr><tr><td>15q22.33</td><td>SMAD3</td><td>SMAD family member 3</td><td>Thyroid</td><td>0.215</td><td>-9.54</td><td>1.45x10<sup>-21</sup></td><td>IVW</td><td>5</td><td>0.63 (0.59-0.68)</td><td>2.21x10<sup>-34</sup></td><td>99.60%</td></tr><tr><td></td><td>AAGAB</td><td>Alpha and gamma adaptin binding protein</td><td>Whole Blood</td><td>0.019</td><td>-7.65</td><td>1.95x10<sup>-14</sup></td><td>h<sub>-</sub></td><td>-</td><td>-</td><td>-</td><td>98.40%</td></tr></tbody></table>

TWAS: Transcriptome-wide association studies; R<sup>2</sup>: R<sup>2</sup> of tissue model's correlation to gene's measured transcriptome; MR: Mendelian randomization; IVW: i variance weighted-random effect; MVMR: Multivariable Mendelian randomization; OR: Odds Ratio; CI: Confidence Interval; COLOC: Posterior probability (PP4 Bayesian colocalization; FOCUS: Fine-mapping Of CaUS gene Sets; PIP: Posterior inclusion probability. <sup>a</sup>Significance threshold (0.05/total number of gene 3.16x10<sup>-6</sup>); <sup>b</sup>Significance threshold (0.05/total number of gene significant after TWAS = 2.78x10<sup>-3</sup>); <sup>c</sup>adjusted for SDCCAG8 in whole blood tissue; <sup>d</sup>adjusted for CCDC32, BAHD1 and ZFYVE19 in thyroid tissue; <sup>e</sup>adjusted for RAD51-AS1, BAHD1 and ZFYVE19 in thyroid tissue; <sup>f</sup>adjusted for ZNF626, ZNF430, ZNF431, ZNF738, ZNF493, ZNF429, ZNF100 and ZNF43 in whole blood tissue; <sup>g</sup>adjusted for TBC1D2 in thyroid tissue; <sup>h</sup>Violation of 1st assumption of Mendelian Randomization; <sup>i</sup>Not found in the FOCUS eQTL dataset.

We then applied MR to the 29 significant genes to identify the mostly causal genes at each 13 loci (Table 2). Eleven genes (PTCSC2, STN1, AAGAB, ACTRT3, TRIM14, TMOD1, SAXO1, LRRIQ4, CHST14, MAFTRR, and RMDN3) were excluded due to a violation of the first assumption of MR (FDR ≥ 0.01). Among the 18 remaining genes (located at 12 loci), nine were located at 6 previously identified loci with causal associations for LRRC34 (3q26.2), NRG1 (8p12), HEMGN (9q22.33), PTCSC3 (14q13.3), and SMAD3 (15q22.33). Additionally, MR analysis on 9 genes located at 7 novel loci reported that predicted expression levels of TGFB2 (1q41), SDCCAG8 (1q43), SAMD4A (14q22.2), RAD51-AS1 (15q15.1) and MPHOSPH6 (16q23.2) were significantly associated with DTC risk.

Integrative analyses using Bayesian colocalization and FOCUS also revealed several genes to prioritize at the 13 TWAS-identified loci (Supplementary Table 3-4). Strong evidence of colocalization were reported for SAMD4A (14q22.2), RAD51-AS1 (15q15.1), MPHOSPH6 (16q23.2), ZNF257 (19p12), ACTRT3 (3q23.2), NRG1 (8p12), and SMAD3 (15q22.3). Moderate colocalization signals were detected for TGFB2 (1q41), SDCCAG8 (1q43), SAXO1 (9p22.1), LRRC34 (3q26.2), PTCSC2 (9q22.33), and STN1 (10q24.33). Notably, SAMD4A, RAD51-AS1, MPHOSPH6, NRG1 and SMAD3 showed strong evidence of causality, supported by both MR and Bayesian colocalization. Despite relatively low PP4 values, genes such as TGFB2 and SDCCAG8 exhibited a strong visual correlation between gene expression and DTC risk, with a coefficient of -0.54 and 0.63, respectively (Supplementary Fig. 4).

Expanding to all available 49 GTEx tissues, we identified 66 significant genes at the Bonferroni-corrected significant level (P = 0.05/20,438 = 2.45x10<sup>-6</sup>) (Supplementary Table 5 and Supplementary Fig. 5). The Z-scores for most of the genes identified were concordant in effect size and p-value compared to thyroid relevant tissues. Apart from thyroid relevant tissues, the tissues with the highest number of detected genes included testis, skin sun-exposed lower leg, and the brain cerebellar hemisphere tissue. Notably, the well-established thyroid cancer risk gene FOXE1 at 9q22.33 (P<sub>TWAS</sub>=5.68x10<sup>-81</sup>) was identified in testis tissue, and the DLG5 gene at 10q22.3 (P<sub>TWAS</sub>=5.77x10<sup>-82</sup>) which initially identified as significant in the GWAS, demonstrated a positive association
between its predicted expression in brain cerebellar hemisphere tissue and DTC. Additionally, the TERT (P<sub>TWAS</sub>=2.61x10<sup>-13</sup>) and CLPTM1L (P<sub>TWAS</sub>=5.68x10<sup>-8</sup>) genes at 5p15, associated with lower predicted expression and increased DTC risk in sun-exposed lower leg skin tissue. Using cis-eQTL from thyroid tumor tissue (significant level of 0.05/5,210 = 9.60x10<sup>-6</sup>), we identified 17 significant genes associated with DTC risk including NFIA gene (P<sub>TWAS</sub>=1.66x10<sup>-10</sup>), which was exclusively identified in thyroid tumor tissue (Supplementary Table 6).

## PWAS analysis

We analyzed 1,475 genetically predicted proteins levels from brain proteome and 1,350 from plasma blood, revealing seven significant proteins under FDR threshold of 0.05: MTHFR (1p36.22), NANS (9q22.33), KDEL C2 (11q22.3), SAMD4A (14q22.2), EIF2AK4 (15q15.1), RMDN3 (15q15.1), and JAG1 (20p12.2) (Supplementary Fig. 6 and Supplementary Table 7). Genetically predicted protein abundance level of MTHFR (OR<sub>IWV</sub>=0.51; 95% CI = 0.33–0.79; P = 2.59x10<sup>-3</sup>), KDEL C2 (OR<sub>MVMR</sub>=0.84; 95% CI = 0.76–0.92; P = 3.48x10<sup>-4</sup>), and SAMD4A (OR<sub>IWV</sub>=1.59; 95% CI = 1.33–1.90; P = 3.36x10<sup>-7</sup>) were causally associated with DTC through MR and Bayesian colocalization (Table 3 and Supplementary Table 8). Notably, the well-known locus 9q22.33 and novel locus 15q15.1 were consistently identified across multi-omics layers, from the genome, transcriptome and proteome. In contrast, SAMD4A gene at 14q22.2 was highlighted only in TWAS (thyroid tissue) and PWAS (brain proteome), without reaching significance in GWAS. Bayesian colocalization supported a shared genetic signal between both eQTL and pQTL for this gene, with PP4 exceeding 90%, further reinforcing its causality through MR analysis (Table 2–3 and Supplementary Fig. 7).

<table><caption>Table 3
Results of the seven genes identified by PWAS followed by Proteome-wide Mendelian randomization and Bayesian posterior probability of colocalization.</caption><thead><tr><th>Locus</th><th>Protein</th><th>Protein Name</th><th>PWAS</th><th>Proteome</th><th>R<sup>2</sup></th><th>Z-score</th><th><sup>a</sup>PWAS P<sub>fdr</sub></th><th>Method</th><th>nsnp</th><th>OR (95% CI)</th><th><sup>b</sup>MR P</th><th>PP4</th></tr></thead><tbody><tr><td colspan="13"><strong>Novel loci</strong></td></tr><tr><td>1p36.22</td><td>MTHFR</td><td>Methylenetetrahydrofolate reductase</td><td>Brain</td><td>0.271</td><td>-4.87</td><td>1.54x10<sup>-03</sup></td><td></td><td>IVW</td><td>4</td><td>0.51 (0.33–0.79)</td><td>2.59x10<sup>-03</sup></td><td>89.40%</td></tr><tr><td>11q22.3</td><td>KDEL C2</td><td>Protein O-glucosyltransferase 3</td><td>Plasma</td><td>0.109</td><td>-4.18</td><td>1.66x10<sup>-02</sup></td><td></td><td><sup>c</sup>MVMR</td><td>12</td><td>0.84 (0.76–0.92)</td><td>3.48x10<sup>-04</sup></td><td>96%</td></tr><tr><td>14q22.2</td><td>SAMD4A</td><td>Sterile alpha motif domain containing 4A</td><td>Brain</td><td>0.353</td><td>4.37</td><td>1.18x10<sup>-02</sup></td><td></td><td>IVW</td><td>2</td><td>1.59 (1.33–1.9)</td><td>3.36x10<sup>-07</sup></td><td>99.10%</td></tr><tr><td>15q15.1</td><td>EIF2AK4</td><td>Eukaryotic translation initiation factor 2 alpha kinase 4</td><td>Brain</td><td>0.089</td><td>3.88</td><td>4.73x10<sup>-02</sup></td><td></td><td><sup>d</sup>MVMR</td><td>4</td><td>2.60 (1.16–5.84)</td><td>2.02x10<sup>-02</sup></td><td>31.10%</td></tr><tr><td></td><td>RMDN3</td><td>Regulator of microtubule dynamics 3</td><td>Brain</td><td>0.215</td><td>3.85</td><td>4.73x10<sup>-02</sup></td><td></td><td>IVW</td><td>2</td><td>3.68 (0.67–20.1)</td><td>0.132</td><td>98.20%</td></tr><tr><td>20p12.2</td><td>JAG1</td><td>Jagged canonical Notch ligand 1</td><td>Plasma</td><td>0.004</td><td>4.17</td><td>1.66x10<sup>-02</sup></td><td></td><td>IVW</td><td>3</td><td>1.56 (1.17–2.09)</td><td>2.58x10<sup>-03</sup></td><td>12.50%</td></tr><tr><td colspan="13"><strong>Known loci</strong></td></tr><tr><td>9q22.33</td><td>NANS</td><td>N-acetylneuraminate synthase</td><td>Brain</td><td>0.344</td><td>5.97</td><td>6.64x10<sup>-06</sup></td><td></td><td><sup>e</sup>MVMR</td><td>4</td><td>3.37(1.86–6.13)</td><td>6.70x10<sup>-05</sup></td><td>0.00%</td></tr></tbody></table>

PWAS: Proteome-wide association studies; R<sup>2</sup>: R<sup>2</sup> of tissue model's correlation to gene's measured transcriptome; MR: Mendelian randomization; IVW: Inversed variance weighted-random effect; MVMR: Multivariable Mendelian randomization; OR: Odds Ratio; CI: Confidence Interval; COLOC: Bayesian colocalization; <sup>a</sup>Significance threshold of false discovery rate of 0.05; <sup>b</sup>Significance threshold (0.05/total number of gene significant after PWAS = 7.14x10<sup>-3</sup>); <sup>c</sup>adjusted for ACAT1 in plasma blood; <sup>d</sup>adjusted for SRP14 in brain; <sup>e</sup>adjusted for CORO2A in brain.

## Enrichment analysis

A gene set enrichment analysis was conducted using 50 genes associated with DTC identified through GWAS, TWAS or PWAS (Supplementary Table 9 and Supplementary Fig. 8). Gene Ontology enrichment analysis revealed significant enrichment in biological processes and molecular functions centered around epithelial-to-mesenchymal transition (EMT)<sup>25</sup>, regulation of miRNA transcription and metabolism<sup>26</sup>, and transforming growth factor beta (TGF-β) signaling. KEGG pathway analysis further highlighted significant enrichment in pathways related to adherens junctions, cell cycle regulation, and the FoxO signaling pathway, all critical for cellular proliferation and survival regulation. Moreover, pathways associated with other cancer types, such as chronic myeloid leukemia and pancreatic cancer, suggest shared oncogenic mechanisms across malignancies.

## Polygenic risk score

We constructed a PRS including previously described and novel GWAS significant SNPs and evaluated its predictive performance in the EPIC and EPITHYR study, where 25 out of 27 SNPs were available. Compared with a previously published 10-SNP PRS<sup>27</sup>, our 25-SNP PRS demonstrated an approximately 2–3% improvement in AUC. Stratification by tumor size revealed that the 25-SNP PRS provided a substantially higher predictive power for larger tumors (> 10 mm
vs ≤ 10 mm), with nearly 6% improvement over 10-SNP PRS (Fig. 3 and Supplementary Table 10), suggesting that PRS may be particularly useful for identifying individuals at higher risk of more advanced or larger carcinomas. Incorporating gene expression and protein abundance as predictors in a linear model led to a slight improvement in AUC (Supplementary Fig. 9). However, considering the computational intensity involved, the gain in predictive performance was minimal.

Discussion

A multi-omics approach integrating genome, transcriptome, and proteome uncovered novel genetic loci associated with DTC risk. Leveraging the largest GWAS to date with a substantial sample size provided robust statistical power, enabling us to identify 10 significant novel loci (1q41, 1q43, 6p21.1, 7q31.33, 8q24.22, 9p22.1, 10q22.3, 15q15.1, 19p13.2 and 22q12.1). We confirmed four previously suggested loci (1p31.3, 5p15.33, 16q23.2, and 19p12)³, along with nine established loci (1q42.2, 2q35, 3q26.2, 5q22.1, 8p12, 9q22.33, 10q24.33, 14q13.3, and15q22.33)⁴. TWAS and PWAS enable us to detect additional potential susceptibility loci at 1p36.22, 11q22.3, 14q22.2 and 20p12.2 not highlighted by GWAS.

Well-established loci

We validated nine previously reported DTC-associated loci, including 1q42.2, 2q35, 3q26.2, 5q22.1, 8p12, 9q22.33, 10q24.33, 14q13.3, and 15q22.33. Some of those loci has been analyzed in previous fine-mapping and/or functional studies³⁰⁻³³. Among these, 9q22.33 locus stood out significant across all omics layers. For instance, NANS was identified in TWAS, alongside other genes such as *PTCS2*, *FOXE1*, *TRMO*, *TRIM14*, *CORO2A*, and *HEMGN*. NANS protein plays a role in amino acid and organic acid synthesis pathways, which are essential for tumor cell metabolism and growth. Interestingly, the effect (Z-scores) of NANS was negative in both normal and tumor tissues, aligning with findings from a recent proteome-wide MR study³⁴ that associated NANS with a decreased risk of DTC based on plasma proteomics. However, our PWAS using brain proteomic data revealed an opposing association, suggesting a potential conflicting effect between plasma and brain proteomic expression levels. The exact function of NANS in thyroid cancer remains unclear and warrants further investigation.

Previously suggested loci

Four loci (1p31.3, 5p15.33, 16q23.2, and 19p12) previously suggested as potential susceptibility genomics regions ($P < 5.0 \times 10^{-6}$), are now confirmed at genome-wide significance in our analysis.

At 1p31.3, we identified rs334723 [G] (OR = 1.36; P = 3.66x10<sup>-11</sup>), within the *NFIA* locus. A previous GWAS found that rs334699[G] (r<sup>2</sup> = 0.97 with rs334723[G]) is associated with decreased TSH level<sup>35</sup>, and another previous MR analysis supports that decreased TSH levels increases DTC risk<sup>36</sup>. At 5p15.33, the lead variant rs7734992[C], (OR = 1.18; P = 1.40 x 10<sup>-16</sup>) is an intronic variant located in *TERT* gene. In the development of thyroid carcinoma, telomerase becomes activated by increasing transcription of the *TERT* gene, leading to elevated levels of *TERT* protein and subsequent activation of telomerase<sup>37</sup>. These promoter alterations are more common in aggressive subtypes such as poorly differentiated and anaplastic thyroid carcinoma<sup>38</sup>. At 16q23.2, we identified rs11645076[T], located within an intronic region. TWAS highlighted *MPHOSPH6* in thyroid-relevant tissues, and *MAFTRR* when considering all tissues. *MPHOSPH6* gene is associated with leucocyte telomere length and other cancers such as glioma or lung cancer<sup>39,40</sup>, though functional studies in thyroid cancer are limited. *MAFTRR* was reported to be over-expressed in patients with Hashimoto's Thyroiditis (HT)<sup>41</sup>. Regarding the 19p12 locus, the lead variant rs8105767[G] maps to the gene *ZNF257*, which shows colocalization between its *cis*-eQTL and DTC, but not causality by MR analysis. It encodes a transcription factor from the zinc finger family, though functional characterization of *ZNF257* in thyroid cancer remains limited. According to GTEx data, *ZNF208* gene, located near rs8105767, exhibits the highest expression in brain hypothalamus and thyroid tissues (Supplementary Fig. 10).

Novel significant loci

Among the 10 novel GWAS-significant loci, six (6p21.1, 7q31.33, 8q24.22, 10q22.3, 19p13.2, and 22q12.1) were identified through GWAS alone, three (1q41, 1q43, and 9p22.1) were further validated by TWAS, and one locus (15q15.1) was identified by all GWAS, TWAS, and PWAS.

Novel loci identified by GWAS only

Experimental studies have indicated that dysregulation of the *VEGFA* gene, located within ± 500kb of the lead variant rs1326141[A] at 6p21.1, may contribute to thyroid cancer progression⁴². The SNP rs11077 (approximately 3,000 bp downstream of rs1326141) was reported to be associated to the expression level of *XPO5* gene and with increased risk of TC in a case-control study conducted in China⁴³. 7q31.33 locus was previously associated to multiple cancers risk such as prostate, chronic lymphotic leukemia, colorectal, glioma⁴⁴,⁴⁵. The gene *POT1* (located 400 kb away from our lead variant rs2299903[G]) plays a critical role in maintaining telomere integrity. A rare germline mutation in *POT1* identified in a non-medullary thyroid cancer family impairs telomere binding and causes telomere elongation in vitro, suggesting that *POT1* dysfunction may elevate thyroid cancer risk through telomere dysregulation⁴⁶. 8q24 locus has been associated to risk for several cancers, including those of the prostate, colon, and ovary⁴⁷. A multi-ethnic case-control study confirmed a significant association with approximately a 13% increased risk of thyroid cancer⁴⁸, despite earlier inconsistent findings regarding its association with thyroid cancer. We reported a considerable number of significant SNPs at 10q22.3 located within *DLG5* and showing high LD with the lead variant 1650149[G]. The *DLG5* gene is identified as a partner gene in *RET* fusion, which commonly occur at the somatic level and are observed predominantly in approximately 2% of lung cancer and 10–20% of thyroid cancer⁴⁹. At 19p13.2, rs4804416[G] located in an intronic region of the *INSR* gene, encodes the insulin receptor. A different SNP, rs919275, in the *INSR* gene (not LD with rs4804416) was weakly associated with PTC risk⁵⁰. Interestingly, this gene was also previously identified in thyroid cancer GWAS in Korean population⁵¹, but did not achieve significance. At 22q12.1, rs186430430[C] (OR = 3.37; P = 1.40x10<sup>-22</sup>) in the intronic region
of the CHEK2 gene has been associated to an increased risk of breast, prostate and colon cancer, while showing a protective association for lung and
laryngeal cancer⁵². The variant is uncommon and relatively rare, with minor allele frequency (MAF) below 1%, and was analyzed only in deCODE genetics and
FinnGen studies, which applied a broader quality control threshold (MAF < 0.01). However, it was not found in other studies (EPITHRY, EPIC, UKBB, and
Italian), even with the same MAF threshold. The significance of the CHEK2 gene appears only in FinnGen study, likely due to the Finnish population's genetic
isolation and unique genetic variations compared to other European population⁵³.

## Novel loci identified by GWAS and TWAS

Three loci (1q41, 1q43, and 9p22.1) were identified through both GWAS and TWAS, with multiple genes located at the same loci. Determining causality among the genes located at the same loci is challenging as they may be associated with disease phenotypes through their correlation with disease-causal genes within the same LD regions. MR analyses using strong IVs and colocalization analyses revealed that expression of *TGFB2* (1q41) and *SDCCAG8* (1q43) were significantly associated with DTC risk.

TGFB2 gene, encoding transforming growth factor beta 2, plays crucial roles in various cellular processes, including cell growth, differentiation, apoptosis, and immune response. Our results align with the TCGA data analysis, which showed an inverse correlation between the expression of TGF-Beta family ligands, including TGFB2, and the thyroid differentiation score⁵⁴,⁵⁵. This suggests that higher expression of these TGF-Beta family ligands is associated with decreased thyroid differentiation, contributing to cancer progression. SDCCAG8 gene, encoding a regulator of sonic hedgehog signaling and ciliogenesis, has been identified as genetically associated with both TSH levels and thyroid diseases⁵⁶. This study highlighted that the TSH-decreasing allele variant in SDCCAG8 (rs10926981[T], r² = 0.72 with our lead SNP rs2451668) was associated with an increased risk of thyroid cancer. At 9p22.1, SAXO1 and HAUS6 genes were highlighted by the TWAS; however, MR analysis was only possible for HAUS6 and yielded no significance, while SAXO1 showed high SNP-level colocalization (PP4) but only moderate gene-level colocalization (PIP).

## Novel loci identified by GWAS, TWAS and PWAS

A novel locus at 15q15.1 has been identified across all omics layers: RAD51 in GWAS, RAD51-AS1 in TWAS (hypothalamus) and RMDN3 (brain) in PWAS. These genes are located in close proximity, within 50 kb of each other, suggesting potential regulatory interactions. RAD51 plays a critical role in homologous recombination, a key DNA repair mechanism. Overexpression of RAD51 has been observed in various cancers, including thyroid carcinoma⁵⁷. Long non-coding RNA RAD51-AS1 has been shown to be down-regulated by E2F1⁵⁸ and RAD51 protein by inhibiting its translation⁵⁹. In this context, down-regulation of RAD51-AS1 could lead to upregulation of both E2F1 and RAD51, potentially promoting thyroid carcinogenesis⁵⁷,⁶⁰. Our MR and eQTL analyses revealed an inverse association between RAD51-AS1 gene expression and a DTC risk, suggesting that RAD51-AS1 gene might mediate the regulation of E2F1 and RAD51 genes. Although the direct role of RMDN3 in thyroid cancer remains to be elucidated, its proximity to RAD51 and RAD51-AS1 at the 15q15.1 locus suggests that it may be part of a regulatory network affecting thyroid carcinogenesis.

## Novel loci identified by TWAS and PWAS

The locus 14q22.2 was significant only in TWAS and PWAS, but not in the GWAS. These findings suggest that SAMD4A may influence DTC risk primarily through regulatory effects on gene expression and protein levels, rather than through direct genetic variants. SAMD4A has been identified as a tumor suppressor in breast cancer, potentially playing a role in inhibiting tumor-induced angiogenesis. Its expression is significantly reduced in breast cancer tissues and correlates strongly with poor patient survival⁶¹, while its role in DTC remains unstudied.

## TWAS findings from diverse tissues

Leveraging all available tissues allowed to detect genes such as TERT and CLPTM1L in sun-exposed lower leg skin tissue, which might have been missed if we had focused solely on thyroid-relevant analysis. TERT promoter mutations, commonly induced by UV radiation, are associated with poorer prognosis in melanoma and are generally found in tumors located in intermittently sun-exposed areas⁶². Previous prospective studies have also shown that factors such as a higher number of nevi and greater residential UV exposure are associated with an increased risk of thyroid cancer⁶³. While the direct effect of TERT mutations on thyroid cancer remains unclear, the shared oncogenic pathways between melanoma and thyroid cancer, such as the RAS-RAF-MEK-ERK signaling pathway, suggest a potential role for these mutations in tumor progression and aggressiveness in thyroid cancer⁶⁴.

## PWAS findings

Proteins detected in plasma or brain may not accurately reflect the same biological processes occurring in thyroid-relevant tissues. PWAS results, based on proteomic data from plasma and brain, has limited statistical power due to a smaller sample size, a weaker biological connection to thyroid cancer, and the restricted number of proteins analyzed in ARIC and ROS/MAP studies. We identified seven proteins coded by genes across six loci. Three loci (1p36.22, 11q22.3, and 20p12.2) not identified in the GWAS and TWAS, were exclusively detected through PWAS.

An interesting gene is *KDELC2*, located on chromosome 11q22.3 alongside *ATM* gene, a key player in telomere maintenance and a master regulator of the DNA double stranded break (DSB) response. Both genes lie within the same LD block, therefore association with *KDELC2* may reflect LD patterns with *ATM*. Indeed, rare pathogenic coding variants in *ATM* have been associated with breast, prostate, thyroid cancer, and pancreatic cancers⁶⁵,⁶⁶. The *JAG1* gene, located at 20p12.2, was found to have significantly elevated expression levels in human thyroid cancer tissues and cell lines compared to normal thyrocytes, as demonstrated by quantitative real-time PCR analysis, suggesting the progression or development of thyroid cancer⁶⁷.

## Enrichment analysis
TGF-β signaling induces EMT and regulates miRNA expression, which modulates the TGF-β pathway and promotes the transition of epithelial cells to a more migratory mesenchymal phenotype. TGF-β plays an important role for the normal functioning of thyroid cells by preventing excessive cellular growth⁶⁸. Also, downregulation of specific miRNAs, such as miR-145, has been linked to increased tumor aggressiveness²,⁶,⁶⁹. Additional enriched processes encompass cellular stress responses, immune regulation, gene expression control, protein kinase activation, and telomere maintenance function. Interestingly, several pathways related to cardiac and pulmonary development were also implicated, suggesting possible shared molecular mechanisms between organogenesis and thyroid carcinogenesis.

## Strengths and limitations

This study presents several key strengths. It is the largest GWAS meta-analysis of DTC to date, with 7,681 cases and 963,550 controls, providing strong statistical power to identify novel risk loci. A major strength lies in the integration of multi-omics approaches, combining *cis*-eQTL, *cis*-pQTL and MR, enabling causal inference and uncovering biological mechanisms underlying DTC risk. An updated PRS improves risk prediction, particularly for larger tumors (> 10 mm), with a ~ 6% increase in detection accuracy.

The study also has limitations. First, restricting primary TWAS analysis to four thyroid-relevant tissues limited detection of well-established DTC-associated genes such as *DIRC3* and *FOXE1*, yet expanding tissue inclusion could introduce false positives by incorporating gene expression from unrelated tissues, potentially identifying non-causal genes due to partial correlation. To mitigate this, we applied MR under strict criteria, including MVMR for multiple genes within a single tissue, though it was not feasible to assess all gene-tissue combinations at a specific locus within a single MR model. Generally MR studies use a stringent clumping threshold ($r^2 < 0.001$) to avoid LD bias, but such thresholds often left only a single SNP per gene⁷⁰, particularly for QTL-based exposures, which limited statistical power and made it difficult to detect and correct for pleiotropy. As a trade-off, we used a relaxed threshold ($r^2 < 0.1$) to increase the number of IVs, which may introduce some correlation among SNPs and weaken conditional F-statistics. Despite this, we observed consistent causal estimates even with relatively low conditional F-statistics (Supplementary Table 11). Secondly, our study focused on *cis*-QTL, excluding *trans*-QTL due to interpretative complexity. Lastly, as the study population was consisted of European ancestry, the generalizability to other ethnic groups is limited. Further researches including participants of Asian and African ancestry are essential to better understand genetic variations in DTC risk across diverse populations.

In conclusion, this multi-omics approach enhances our understanding of how genetic variation, gene expression, and protein abundance contribute to DTC risk. The identification of the novel locus 14q22.2, significant only in TWAS and PWAS, demonstrates the value of this integrative strategy. Our findings emphasize the need for further functional studies to validate these genes and their biological roles in DTC.

## Declarations

### Acknowledgments

We acknowledge the use of data and biological samples from the EPIC-Ragusa cohort, principal investigator Rosario Tumino; EPIC-Asturias, principal investigator José Ramón Quirós García; EPIC-Bilthoven, principal investigator Monique Verschuren, EPIC-Utrecht, principal investigator Roel Vermeulen.

We acknowledge deCODE Genetics, and FinnGen for sharing their GWAS summary statistics. This research has been conducted using the UK Biobank Resource under application number 92392. This work is also part of the Inserm Cross-Cutting Project GOLD.

### Funding

See Hyun Park was the recipient of a PhD fellowship from Paris-Saclay University.

The coordination of EPIC is financially supported by International Agency for Research on Cancer (IARC) and also by the Department of Epidemiology and Biostatistics, School of Public Health, Imperial College London which has additional infrastructure support provided by the NIHR Imperial Biomedical Research Centre (BRC). The national cohorts from EPIC are supported by: Danish Cancer Society (Denmark); Ligue Contre le Cancer, Institut Gustave Roussy, Mutuelle Générale de l'Éducation Nationale, Institut National de la Santé et de la Recherche Médicale (INSERM) (France); German Cancer Aid, German Cancer Research Center (DKFZ), German Institute of Human Nutrition Potsdam-Rehbruecke (DIfE), Federal Ministry of Education and Research (BMBF) (Germany); Associazione Italiana per la Ricerca sul Cancro-AIRC-Italy, Compagnia di SanPaolo and National Research Council (Italy); Dutch Ministry of Public Health, Welfare and Sports (VWS), Netherlands Cancer Registry (NKR), the National Institute for Public Health and the Environment (RIVM) (The Netherlands), LK Research Funds, Dutch Prevention Funds, Dutch ZON (Zorg Onderzoek Nederland), World Cancer Research Fund (WCRF), Statistics Netherlands (The Netherlands); Health Research Fund (FIS) - Instituto de Salud Carlos III (ISCIII), Regional Governments of Andalucía, Asturias, Basque Country, Murcia and Navarra, and the Catalan Institute of Oncology - ICO (Spain); Cancer Research UK (14136 to EPIC-Norfolk; C8221/A29017 to EPIC-Oxford), Medical Research Council (1000143 to EPIC-Norfolk; MR/M012190/1 to EPIC-Oxford). (United Kingdom).

The genotyping of EPIC samples was supported by the Association pour la Recherche sur le Cancer (ARC) (#RF20180207126). The EPITHYR genome-wide association study was funded by INCA (#9533) and ARC (#PGA120150202302).

### Author Disclosure Statement

None of the authors have any disclosure to report nor competing financial interest.

### Conflict of interest

The other authors declare no conflict of interest.
## Author contributions

S.P. and T.T designed the study. S.P., Y.A., P.E.S., M.K., and T.T conducted statistical analyses. H.T., A.F., R.E., F.G., S.L., F.d.V., F.L., C.M., P.L.P., E.O., P.G., J.F.D., C.C.D., F.C., G.I., C.S., M.G., A.J., K.S.B., J.Y., S.R., G.S and TT provided the data. Principal collaborators were F.V., F.L., C.M., P.L.P., E.O., P.G., and J.F.D., for EPITHYR; C.C.D., F.C., G.I., C.S., M.G., A.J., K.S.B., J.Y., S.R., and G.S for EPIC; H.T., A.F., R.E., F.G., and S.L., for the Italian study. All authors contributed to the review and editing. All authors contributed read and approved the final version of the manuscript.

## Data availability

Individual level data are not publicly available due to GDPR-related restrictions on personal data sharing. Data from EPIC access can be requested via https://epic.iarc.fr/access/index.php. All requests are evaluated by the EPIC Steering Committee. Data from EPITHR may be requested from the corresponding author and will be evaluated by the data access committee. Data from the UK Biobank can be accessed through their website upon request (https://www.ukbiobank.ac.uk/), FinnGen GWAS data are publicly available at https://www.finngen.fi/en and summary statistics from deCODE genetics can be accessed at https://www.decode.com/summarydata/. The R scripts used for analyses can be provided upon reasonable request to the corresponding author.

## Disclaimer

Where authors are identified as personnel of the International Agency for Research on Cancer / World Health Organization, the authors alone are responsible for the views expressed in this article and they do not necessarily represent the decisions, policy or views of the International Agency for Research on Cancer / World Health Organization.

## References

1. Olson, E., Wintheiser, G., Wolfe, K. M., Droessler, J. & Silberstein, P. T. Epidemiology of Thyroid Cancer: A Review of the National Cancer Database, 2000-2013. *Cureus* **11**, e4127.

2. Hemminki, K. & Li, X. Familial risk of cancer by site and histopathology. *International Journal of Cancer* **103**, 105-109 (2003).

3. Truong, T. et al. Multiethnic genome-wide association study of differentiated thyroid cancer in the EPITHYR consortium. *International Journal of Cancer* **148**, 2935-2946 (2021).

4. Gudmundsson, J. et al. A genome-wide association study yields five novel thyroid cancer risk loci. *Nat Commun* **8**, 14517 (2017).

5. Sanderson, E. et al. Mendelian randomization. *Nat Rev Methods Primers* **2**, 6 (2022).

6. Köhler, A. et al. Genome-wide association study on differentiated thyroid cancer. *J Clin Endocrinol Metab* **98**, E1674-1681 (2013).

7. Riboli, E. et al. European Prospective Investigation into Cancer and Nutrition (EPIC): study populations and data collection. *Public Health Nutr* **5**, 1113-1124 (2002).

8. Kurki, M. I. et al. FinnGen provides genetic insights from a well-phenotyped isolated population. *Nature* **613**, 508-518 (2023).

9. Mägi, R. & Morris, A. P. GWAMA: software for genome-wide association meta-analysis. *BMC Bioinformatics* **11**, 1-6 (2010).

10. Lonsdale, J. et al. The Genotype-Tissue Expression (GTEx) project. *Nat Genet* **45**, 580-585 (2013).

11. Zhou, D. et al. A unified framework for joint-tissue transcriptome-wide association and Mendelian randomization analysis. *Nat Genet* **52**, 1239-1246 (2020).

12. Barbeira, A. N. et al. Exploring the phenotypic consequences of tissue specific gene expression variation inferred from GWAS summary statistics. *Nat Commun* **9**, 1825 (2018).

13. Vogel, C. & Marcotte, E. M. Insights into the regulation of protein abundance from proteomic and transcriptomic analyses. *Nat Rev Genet* **13**, 227-232 (2012).

14. Bennett, D. A. et al. Religious Orders Study and Rush Memory and Aging Project. *J Alzheimers Dis* **64**, S161-S189 (2018).

15. Zhang, J. et al. Plasma proteome analyses in individuals of European and African ancestry identify cis-pQTLs and models for proteome-wide association studies. *Nat Genet* **54**, 593-602 (2022).

16. Wingo, A. P. et al. Integrating human brain proteomes with genome-wide association data implicates new proteins in Alzheimer's disease pathogenesis. *Nat Genet* **53**, 143-146 (2021).

17. Gusev, A. et al. Integrative approaches for large-scale transcriptome-wide association studies. *Nat Genet* **48**, 245-252 (2016).

18. Burgess, S., Zuber, V., Valdes-Marquez, E., Sun, B. B. & Hopewell, J. C. Mendelian randomization with fine-mapped genetic data: Choosing from large numbers of correlated instrumental variables. *Genet Epidemiol* **41**, 714-725 (2017).

19. Sanderson, E. Multivariable Mendelian Randomization and Mediation. *Cold Spring Harb Perspect Med* **11**, a038984 (2021).

20. Burgess, S., Butterworth, A. & Thompson, S. G. Mendelian Randomization Analysis With Multiple Genetic Variants Using Summarized Data. *Genetic Epidemiology* **37**, 658-665 (2013).

21. Giambartolomei, C. et al. Bayesian Test for Colocalisation between Pairs of Genetic Association Studies Using Summary Statistics. *PLoS Genet* **10**, e1004383 (2014).

22. Mancuso, N. et al. Probabilistic fine-mapping of transcriptome-wide association studies. *Nat Genet* **51**, 675-682 (2019).

23. Ashburner, M. et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. *Nat Genet* **25**, 25-29 (2000).
24. Kanehisa, M. & Goto, S. KEGG: Kyoto Encyclopedia of Genes and Genomes. *Nucleic Acids Res* **28**, 27–30 (2000).

25. Shakib, H. et al. Epithelial-to-mesenchymal transition in thyroid cancer: a comprehensive review. *Endocrine* **66**, 435–455 (2019).

26. Hamidi, A. A., Taghehchian, N., Basirat, Z., Zangouei, A. S. & Moghbeli, M. MicroRNAs as the critical regulators of cell migration and invasion in thyroid cancer. *Biomarker Research* **10**, 40 (2022).

27. Liyanarachchi, S. et al. Assessing thyroid cancer risk using polygenic risk scores. *Proc Natl Acad Sci U S A* **117**, 5997–6002 (2020).

28. Comiskey, D. F. et al. Characterizing the function of EPB41L4A in the predisposition to papillary thyroid carcinoma. *Sci Rep* **10**, 19984 (2020).

29. Shen, Z., Sun, Y. & Niu, G. Variants in TPO rs2048722, PTCSC2 rs925489 and SEMA4G rs4919510 affect thyroid carcinoma susceptibility risk. *BMC Med Genomics* **16**, 19 (2023).

30. He, H. et al. The Role of NRG1 in the Predisposition to Papillary Thyroid Carcinoma. *J Clin Endocrinol Metab* **103**, 1369–1379 (2017).

31. Guibon, J. et al. Fine-mapping of two differentiated thyroid carcinoma susceptibility loci at 2q35 and 8p12 in Europeans, Melanesians and Polynesians. *Oncotarget* **12**, 493–506 (2021).

32. Tcheandjieu, C. et al. Fine-mapping of two differentiated thyroid carcinoma susceptibility loci at 9q22.33 and 14q13.3 detects novel candidate functional SNPs in Europeans from metropolitan France and Melanesians from New Caledonia. *Int J Cancer* **139**, 617–627 (2016).

33. Jendrzejewski, J. et al. Fine mapping of 14q13 reveals novel variants associated with different histological subtypes of papillary thyroid carcinoma. *Int J Cancer* **144**, 503–512 (2019).

34. Fan, Q. et al. Assessment of circulating proteins in thyroid cancer: Proteome-wide Mendelian randomization and colocalization analysis. *iScience* **27**, (2024).

35. Malinowski, J. R. et al. Genetic Variants Associated with Serum Thyroid Stimulating Hormone (TSH) Levels in European Americans and African Americans from the eMERGE Network. *PLoS One* **9**, e111301 (2014).

36. Yuan, S. et al. Causal associations of thyroid function and dysfunction with overall, breast and thyroid cancer: A two-sample Mendelian randomization study. *Int J Cancer* **147**, 1895–1903 (2020).

37. Yuan, X., Liu, T. & Xu, D. Telomerase reverse transcriptase promoter mutations in thyroid carcinomas: implications in precision oncology—a narrative review. *Annals of Translational Medicine* **8**, 1244–1244 (2020).

38. Yuan, X., Yuan, H., Zhang, N., Liu, T. & Xu, D. Thyroid carcinoma-featured telomerase activation and telomere maintenance: Biology and translational/clinical significance. *Clin Transl Med* **12**, e1111 (2022).

39. Cui, Y. et al. The effects of gene polymorphisms on glioma prognosis. *J Gene Med* **19**, 345–352 (2017).

40. Cortez Cardoso Penha, R. et al. Common genetic variations in telomere length genes and lung cancer: a Mendelian randomisation study and its novel application in lung tumour transcriptome. *Elife* **12**, e83118 (2023).

41. Peng, H. et al. Elevated Expression of the Long Noncoding RNA MAFTRR in Patients with Hashimoto's Thyroiditis. *J Immunol Res* **2021**, 3577011 (2021).

42. Stuchi, L. P. et al. VEGFA and NFE2L2 Gene Expression and Regulation by MicroRNAs in Thyroid Papillary Cancer and Colloid Goiter. *Genes* **11**, (2020).

43. Wen, J. et al. Association of microRNA-related gene XPO5 rs11077 polymorphism with susceptibility to thyroid cancer. *Medicine (Baltimore)* **96**, e6351 (2017).

44. Paiss, T. et al. Linkage of aggressive prostate cancer to chromosome 7q31-33 in German prostate cancer families. *Eur J Hum Genet* **11**, 17–22 (2003).

45. Saunders, C. N. et al. Relationship between genetically determined telomere length and glioma risk. *Neuro Oncol* **24**, 171–181 (2021).

46. Srivastava, A. et al. A Germline Mutation in the POT1 Gene Is a Candidate for Familial Non-Medullary Thyroid Cancer. *Cancers (Basel)* **12**, 1441 (2020).

47. Ahmadiyeh, N. et al. 8q24 prostate, breast, and colon cancer risk loci show tissue-specific long-range interaction with MYC. *Proc Natl Acad Sci U S A* **107**, 9742–9746 (2010).

48. Sahasrabudhe, R. et al. 8q24 rs6983267G variant is associated with increased thyroid cancer risk. *Endocr Relat Cancer* **22**, 841–849 (2015).

49. Santoro, M., Moccia, M., Federico, G. & Carlomagno, F. RET Gene Fusions in Malignancies of the Thyroid and Other Tissues. *Genes (Basel)* **11**, 424 (2020).

50. Kitahara, C. M. et al. Common obesity-related genetic variants and papillary thyroid cancer risk. *Cancer Epidemiol Biomarkers Prev* **21**, 2268–2271 (2012).

51. Son, H.-Y. et al. Genome-wide association and expression quantitative trait loci studies identify multiple susceptibility loci for thyroid cancer. *Nat Commun* **8**, 15966 (2017).

52. Cybulski, C. et al. Constitutional CHEK2 mutations are associated with a decreased risk of lung and laryngeal cancers. *Carcinogenesis* **29**, 762–765 (2008).

53. Översti, S. et al. Human mitochondrial DNA lineages in Iron-Age Fennoscandia suggest incipient admixture and eastern introduction of farming-related maternal ancestry. *Sci Rep* **9**, 16883 (2019).

54. Han, H. et al. Expression and Prognostic Value of m6A RNA Methylation-Related Genes in Thyroid Cancer. *Iran J Public Health* **52**, 1902–1916 (2023).

55. Alves, L. F. & Geraldo, M. V. MiR-495-3p regulates cell migration and invasion in papillary thyroid carcinoma. *Front Oncol* **13**, 1039654 (2023).

56. Williams, A. T. et al. Genome-wide association study of thyroid-stimulating hormone highlights new genes, pathways and associations with thyroid disease. *Nat Commun* **14**, 6713 (2023).

57. Sarwar, R. et al. Upregulation of RAD51 expression is associated with progression of thyroid carcinoma. *Exp Mol Pathol* **102**, 446–454 (2017).

58. Gazy, I. et al. TODRA, a lncRNA at the RAD51 Locus, Is Oppositely Regulated to RAD51, and Enhances RAD51-Dependent DSB (Double Strand Break) Repair. *PLoS One* **10**, e0134120 (2015).
59. Chen, C.-C. et al. Corylin increases the sensitivity of hepatocellular carcinoma cells to chemotherapy through long noncoding RNA RAD51-AS1-mediated inhibition of DNA repair. *Cell Death Dis* **9**, 1–13 (2018).

60. Onda, M. et al. Up-regulation of transcriptional factor E2F1 in papillary and anaplastic thyroid cancers. *J Hum Genet* **49**, 312–318 (2004).

61. Zhou, M. et al. RNA-binding protein SAMD4A inhibits breast tumor angiogenesis by modulating the balance of angiogenesis program. *Cancer Sci* **112**, 3835–3845 (2021).

62. Pópulo, H. et al. TERT promoter mutations in skin cancer: the effects of sun exposure and X-irradiation. *J Invest Dermatol* **134**, 2251–2257 (2014).

63. Mesrine, S. et al. Nevi, Ambient Ultraviolet Radiation, and Thyroid Cancer Risk: A French Prospective Study. *Epidemiology* **28**, 694–702 (2017).

64. Zerfaoui, M. et al. New Insights into the Link between Melanoma and Thyroid Cancer: Role of Nucleocytoplasmic Trafficking. *Cells* **10**, 367 (2021).

65. Kang, J., Deng, X.-Z., Fan, Y.-B. & Wu, B. Relationships of FOXE1 and ATM genetic polymorphisms with papillary thyroid carcinoma risk: a meta-analysis. *Tumor Biol.* **35**, 7085–7096 (2014).

66. Concannon, P. et al. Variants in the ATM gene associated with a reduced risk of contralateral breast cancer. *Cancer research* **68**, 6486 (2008).

67. Chen, J., Wang, X., Zhang, X., Yin, J. & Zheng, Y. Jagged 1 Regulates The Proliferation and Metastasis of Human MDA-T68 Thyroid Cancer Cells. *Cell J* **25**, 399–406 (2023).

68. Colletta, G., Cirafici, A. M. & Di Carlo, A. Dual effect of transforming growth factor beta on rat thyroid cells: inhibition of thyrotropin-induced proliferation and reduction of thyroid-specific differentiation markers. *Cancer Res* **49**, 3457–3462 (1989).

69. Toraih, E. A. et al. A miRNA-Based Prognostic Model to Trace Thyroid Cancer Recurrence. *Cancers (Basel)* **14**, 4128 (2022).

70. Richardson, T. G., Hemani, G., Gaunt, T. R., Relton, C. L. & Davey Smith, G. A transcriptome-wide Mendelian randomization study to uncover tissue-dependent regulatory mechanisms across the human phenome. *Nat Commun* **11**, 185 (2020).

## Figures

<figure><img src="image_3.png" /><figcaption>Figure 1</figcaption></figure>

Legend not included with this version
<figure><img src="image_4.png" /><figcaption>Figure 2</figcaption></figure>

Legend not included with this version
<figure><img src="image_5.png" /><figcaption>Figure 3</figcaption></figure>

Legend not included with this version

Supplementary Files

This is a list of supplementary files associated with this preprint. Click to download.

6.GWASTWASPWASSupplementaryTableS1.xlsx

6.GWASTWASPWASSupplementaryTableS2.xlsx

6.GWASTWASPWASSupplementaryTableS3.xlsx

6.GWASTWASPWASSupplementaryTableS4.xlsx

6.GWASTWASPWASSupplementaryTableS5.xlsx
6.GWASTWASPWASSupplementaryTableS6.xlsx

6.GWASTWASPWASSupplementaryTableS7.xlsx

6.GWASTWASPWASSupplementaryTableS8.xlsx

6.GWASTWASPWASSupplementaryTableS9.xlsx

6.GWASTWASPWASSupplementaryTableS10.xlsx

6.GWASTWASPWASSupplementaryTableS11.xlsx

5.GWASTWASPWASSupplementaryMethodFigures.pdf