Edinburgh Research Explorer Genome-wide analyses identify a role for SLC17A4 and AADAT in thyroid hormone regulation

Thyroid dysfunction is an important public health problem, which affects 10% of the general population and increases the risk of cardiovascular morbidity and mortality. Many aspects of thyroid hormone regulation have only partly been elucidated, including its transport, metabolism, and genetic determinants. Here we report a large meta-analysis of genome-wide association studies for thyroid function and dysfunction, testing 8 million genetic variants in up to 72,167 individuals. One-hundred-and-nine independent genetic variants are associated with these traits. A genetic risk score, calculated to assess their combined effects on clinical end points, shows signi ﬁ cant associations with increased risk of both overt (Graves ’ disease) and subclinical thyroid disease, as well as clinical complications. By functional follow-up on selected signals, we identify a novel thyroid hormone transporter (SLC17A4) and a metabolizing enzyme (AADAT). Together, these results provide new knowledge about thyroid hormone physiology and disease, opening new possibilities for therapeutic targets. Theo J. Visser, Marco Medici. Deceased: Theo J. Visser. A full list of consortium members appears below.

T hyroid dysfunction is a common clinical condition, affecting~10% of the general adult population 1 . Adequate thyroid hormone levels are essential for normal growth and differentiation, regulation of energy metabolism, and physiological function of virtually all human tissues. Thyroxine (T4) is the prohormone produced by the thyroid, which is largely converted into the active hormone 3,3′,5-triiodothyronine (T3) in peripheral tissues. Circulating T4 levels are regulated by the hypothalamus-pituitary-thyroid (HPT) axis, in which pituitary thyroid-stimulating hormone (TSH) stimulates T4 production. In turn, T4 and T3 negatively regulate TSH synthesis via a negative feedback loop.
To exert their actions, T4 and T3 cross the membranes of target cells via specific transporters. Once intracellular, they are metabolized, including the conversion of T4 to T3, followed by binding of T3 to its nuclear receptor to regulate transcription of target genes. Both T4 and T3 transport and metabolism are therefore key determinants of thyroid hormone action.
In daily clinical practice, thyroid function is assessed by measuring circulating TSH and free T4 (FT4) levels, with increased TSH indicating hypothyroidism and decreased TSH indicating hyperthyroidism. FT4 levels are decreased in overt hypothyroidism, increased in overt hyperthyroidism and in the reference range in subclinical hypo and hyperthyroidism. In the last decade, it has become clear that not only overt but also subclinical hypo and hyperthyroidism are associated with several pathological conditions, such as atrial fibrillation, coronary heart disease, stroke, depression, as well as cardiovascular and overall mortality [2][3][4][5][6][7] . More recently, studies have shown that even variation in thyroid function within the normal range is associated with many of these complications 4,[8][9][10] . Despite the physiological significance of thyroid hormones, as well as the prevalence and clinical importance of thyroid dysfunction, many key players in the regulation of thyroid hormone bioavailability and action, including its transport and metabolism, still need to be elucidated.
Genome-wide association studies (GWAS) performed so far have revealed genetic variants in about 30 loci robustly associated with thyroid function [11][12][13] . However, these variants explain only <9% of the heritability of TSH and FT4 variation 14 , while in total, it has been estimated at 65 and 39-80% for TSH and FT4, respectively 15,16 , suggesting that many loci still await discovery.
Here, we report the results of a large meta-analysis of GWAS for circulating TSH and FT4 levels, as well as for hypo and hyperthyroidism, followed by independent replication and functional studies. Results are complemented with genetic risk score (GRS) analyses, gene expression, co-localization analyses, and associations with various clinical phenotypes (Supplementary Figure 1) to discover new pathways underlying thyroid function and disease. We identify 109 significant independent genetic associations with these traits. The GRS shows a significant association with increased risk of both Graves' disease and subclinical thyroid disease, as well as clinical complications. Finally, we identify a novel thyroid hormone transporter and a metabolizing enzyme. Together, these results enhance our knowledge about thyroid hormone physiology and disease.

Results
New loci affecting thyroid hormone levels. Our GWAS metaanalyses and replication in up to 72,167 subjects of European ancestry with TSH levels within the reference range (Supplementary Data 1) discovered 19 novel loci for circulating TSH levels and 16 novel loci for circulating FT4 levels (Tables 1 and 2, Supplementary Figures 2-5), leading to a total of 42 and 21 known and novel associated loci for these two traits. As illustrated in Fig. 1, TSH and FT4 capture distinct and complementary genetic underpinnings of thyroid function. Some of the novel loci include genes that have been previously implicated in thyroid development (GLIS3), thyroid hormone action and transport (NCOR1, TTR, SLCO1B1), thyroid hormone metabolism (DIO2, DIO3OS), and thyroid cancer (e.g., HES1, SPATA13, DIRC3, ID4) by various candidate gene studies of monogenic diseases and animal models. Multiple independent variants were found for PDE8B, DIO1, DIO2, TSHR, and CAPZB.
Across the 42 TSH and 21 FT4 loci, allelic heterogeneity (i.e., independently associated single-nucleotide polymorphisms (SNPs) at the same locus) was detected at 11 and 7 loci, respectively, by using linkage structure information and summary statistics-based conditional analyses (Supplementary Tables 1 and  2). All significant associations together accounted for 33% and 21% of the genetic variance of TSH and FT4, respectively, explained by all common and low-frequency variants with a minor allele frequency (MAF) >1%.
Since TSH and FT4 regulation are inversely correlated through the HPT axis, we investigated the association of the TSHassociated loci with FT4 levels, and vice versa. As shown in Fig. 1 and in Supplementary Tables 1 and 2, we observed overlapping associations (Bonferroni-corrected threshold p < 8.2 × 10 −4 ) at various loci (TSH: FGF7, PDE8B, DET1, ITPK1, VEGFA, GLIS3, NFIA, and MBIP, and FT4: FOXE1 and GLIS3), although only GLIS3 showed genome-wide significance for both traits. All alleles associated with higher TSH were associated with lower FT4, with the exception of MBIP and FOXE1.
Hypo and hyperthyroidism are more prevalent in women than in men. However, sex-stratified GWAS meta-analyses for TSH and FT4 did not show any significant gene-by-sex interaction in our samples (Supplementary Figure 6, Supplementary Tables 1  and 2).
Given the high degree of functional homology between the mouse and human genome, we selected from The International Mouse Phenotyping Consortium database 17 genes that when manipulated in mice cause abnormal thyroid physiology (i.e., hormone levels, n = 26) or morphology (n = 51), and assessed whether their human homologs contained SNPs significantly (p < 2.5 × 10 −5 for physiology and p < 1.9 × 10 −5 for morphology, see Methods) associated in our FT4 and TSH GWAS (Supplementary Data 2). Of these candidate genes, SNPs in CGA (rs6924373) and TPO (rs9678281) contained significant associations that did not reach genome-wide significance in our GWAS for TSH. These associations were tested for replication in 9011 independent samples and achieved genome-wide significance for TSH (p < 5 × 10 −8 ) in the combined dataset (Supplementary Figure 7A). Overall, these results highlight the potential of nested candidate gene approaches in GWAS summary results and emphasize the functional conservation of genes regulating thyroid function between mice and humans.
Relation to hypo and hyperthyroidism. Genetic variants that determine variation in circulating TSH and FT4 levels within the reference range (i.e., the individual HPT-axis setpoint) are expected to differ from variants that underlie thyroid dysfunction (hypo or hyperthyroidism). To clarify this, we also conducted a case-control GWAS meta-analysis of increased TSH levels (i.e., hypothyroidism), including cases with TSH levels above the cohort-specific reference range (n = 3340) and controls with TSH levels within the reference range (n = 49,983). The decreased TSH level (i.e., hyperthyroidism) GWAS meta-analysis included cases with TSH below the reference range (n = 1840 cases) and the same controls as in the increased TSH GWAS. The distribution of sex and age groups of these subjects is provided in Supplementary  Table 3. Since in both GWAS analyses, cases were defined on the ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-06356-1 2 basis of a TSH level above or below the reference range, these groups included subjects with overt but also mild subclinical forms of hypothyroidism and hyperthyroidism, respectively.
We detected seven loci for hypothyroidism and eight loci for hyperthyroidism (Supplementary Figure 8, Supplementary  Table 4). At some of the loci, the variant was significantly associated with both hypo and hyperthyroidism, with effects in opposing directions. For example, a variant at PDE10A (rs2983514) was associated with both higher risk of hypothyroidism and lower risk of hyperthyroidism. Some of the hypothyroidism loci had already previously been implicated in hypothyroidism through GWAS, including TPO, FOXE1, VAV3, and a variant in ATXN2 (rs597808) in high linkage disequilibrium (LD) with the R262W polymorphism in SH2B3 18 . However, we did not detect variants in a number of well-known autoimmune thyroid disease genes (e.g., CTLA4, HLA class I and II). This may be due to the fact that, in these population-based cohorts, patients receiving medication for autoimmune thyroiditis were excluded. Thus, thyroid autoimmunity caused by auto-antibodies may have a different set of predisposing variants. All variants associated with hyperthyroidism have not been previously found in association with hyperthyroidism, except for FOXE1 19 . However, all of these variants were in high LD with variants associated with TSH or FT4 levels within the reference range in the current or previous GWAS 11,13 . The same holds true for variants associated with hypothyroidism, suggesting that the effects of many genetic variants on thyroid function extend beyond the physiological range, thus affecting the risk of thyroid dysfunction.
As complementary analyses to investigate whether the TSH loci are also related to autoimmune thyroid diseases, we tested all variants or their proxies for association with thyroid peroxidase antibody (TPOAb) positivity of a former GWAS 20 as an early   Table 5). SPATA13 and VAV3 have previously been linked to self-reported diagnosed hypothyroidism 18,21 , while no studies have investigated their potential autoimmune origin. The observed associations of these gene variants with variation in TSH levels within the normal range could therefore be due to a mild early stage of thyroid autoimmunity, instead of reflecting physiological differences in the HPT-axis setpoint. For Graves' disease, only the psoriasis 22 PSORS1C1 locus showed a significant association, consistent with shared genetic determinants between these two autoimmune diseases 23 .
Detailed results of the mouse candidate analysis, results of pathway analyses, and look-ups for pleiotropy of the TSH, FT4, hypo and hyperthyroidism loci are described in Supplementary Note 1-3.
Gene expression analyses. To obtain insights into gene expression patterns and potential effector transcripts at the identified loci, we assessed whether the 94 independent index variants from the TSH, FT4, hypo and hyperthyroidism GWAS were correlated with transcript levels of nearby (cis-) or distant (trans-) genes. The results of 22 published expression quantitative trait loci (eQTL) studies were interrogated (Methods), assessing the relation between the genetic variants and gene expression patterns in a total of 127 different tissues and cell types. First, we evaluated the presence of eQTLs in at least one tissue or cell type: 38 variants showed eQTL effects ( Fig. 2a and Supplementary Data 3). While many variants were associated with transcript expression in only one or few (≤8) tissues, two variants located on chromosome 17 showed ubiquitous associations with gene expression: the FT4-associated variant rs11078333 at NCOR1 locus and the TSH-associated variant rs199461 at the NSF locus. The FT4-increasing allele at rs11078333 was associated with higher expression levels of NCOR1 in blood and brain, but also affected the expression of ADORA2B (increased) and ZSWIM7 and TTC19 (decreased) in many other tissues (including thyroid for TTC19). NCOR1 is an essential nuclear co-repressor that is recruited by thyroid hormone receptors in the absence of thyroid hormone to mediate transcriptional repression. At the NSF locus, the TSH increasing allele at rs199461 increased expression of KANSL1 and LRRC37A2 and decreased expression of WNT3 in several tissues, including thyroid. Consistent with known thyroid physiology, the majority of TSH-associated variants acted as eQTLs in thyroid tissue (Fig. 2a), with 45% of the variants being thyroid-specific eQTLs. In contrast, none of the nine FT4 eQTLs acted exclusively on the thyroid but were also associated with transcript expression changes in multiple known thyroid hormone effector organs, including liver, muscle, and adipose tissue.
Second, we used a summary-based Mendelian randomization (SMR) method coupled with testing for heterogeneity of effects (HEIDI) to assess co-localization, i.e., to investigate whether the overlap between eQTLs and GWAS hits could be attributable to the same underlying causative variant 24 . In thyroid tissue, we found evidence for co-localization with differential gene expression at 13 different GWAS loci: PD8EB, PRDM11, MBIP (with RP11-116N8.1 expression), NKX2-3, NSF (with WNT3 expression), IGF2BP2, FOXA2, SLC25A37, and C9orf92 for TSH; AADAT, NEK6 (with both NEK6 and PSMB7 expressions) for FT4; and TPO, PDE8B, and PDE10A for hypothyroidism (Supplementary Data 4). At these loci, our findings implicate the causal gene among the many genes present in the locus. For example, the FT4-associated variant rs6854291-influenced transcript levels of AADAT in thyroid while there were no effects on 50   . We also observed that independent variants at the same locus were associated with gene expression in different tissues. For example, while the index variant rs6885099 at PDE8B co-localized with changes in PDE8B expression in thyroid, the independent variant rs1119208 was associated with PDE8B expression in pancreas (Supplementary Data 4). Notably, also the variants at AADAT and SLC17A4 colocalized with gene expression in pancreas (Fig. 2c), which is of interest given the close interrelations between thyroid hormone signaling, insulin regulation, and glucose metabolism 25,26 .  In vitro studies. Thyroid hormone action in target tissues is importantly determined by the amount of T3 available for receptor binding inside the cell. Therefore, the transport of thyroid hormone across the cell membrane and its metabolism inside the cell represent crucial regulatory layers in thyroid hormone signaling. Although several key players in thyroid hormone signaling have been described over the last decades, including deiodinases and several thyroid hormone transporters, many others remain to be identified. Based on their associations with circulating FT4 levels and the co-localization studies, SLC17A4 and AADAT were further studied in vitro to explore a direct role in thyroid hormone signaling. SLC17A4 is an organic anion transporter that is particularly expressed in the liver, kidney, and gastrointestinal tract 27 . We transiently over-expressed human SLC17A4 (hSLC17A4) in COS-1 cells and observed increased cellular T3 (Fig. 3a) and T4 (Fig. 3b) accumulation compared to empty-vector transfected control cells. These effects were even stronger upon cotransfection with the intracellular thyroid hormone-binding protein mu-crystallin (CRYM) (Fig. 3a, b) and were similar in magnitude to those obtained by the monocarboxylate transporter (MCT) 8 ( Supplementary Figure 9), the most specific thyroid hormone transporter identified to date. Saturation experiments in the absence of CRYM showed a dose-dependent decrease in the uptake of T3 (Fig. 3c) and T4 (Fig. 3d). The estimated IC 50 values for T3 (0.35 ± 0.13 μM, n = 4) and T4 (0.06 ± 0.01 μM, n = 4) transport by SLC17A4 are considerably lower than those of MCT8 (T3: 20.61 ± 1.26 μM and T4: 23.22 ± 1.22 μM, n = 3, Supplementary Figure 9), and other currently known thyroid hormone transporters [28][29][30][31][32][33] , and indicate a high substrate affinity.
Together, these findings strongly indicate that SLC17A4 encodes a high-affinity T3 and T4 transporter.
AADAT encodes a mitochondrial aminotransferase with broad substrate specificity, which acts on kynurenic acid and α-aminoadipate, important intermediates in tryptophan, and lysine metabolism 34,35 . The association of circulating FT4 with the AADAT locus suggested that AADAT may also be involved in thyroid hormone metabolism. In that case, it could facilitate the oxidative deamination of the alanine side-chain of thyroid hormone, yielding a pyruvic acid moiety 36 . Therefore, lysates of AADAT over-expressing COS-1 cells were incubated with T4 and T3 in the presence of the co-factor pyridoxal phosphate and the co-substrate α-ketoglutaric acid, and the reaction mixtures were analyzed by ultra-performance liquid chromatography (UPLC). The results demonstrated effective timeand AADAT concentration-dependent conversion of T4 and T3 to their pyruvic acid metabolites TK3 and TK4 (Fig. 4), with saturation occurring at substrate concentrations between 10 and 100 μM. Importantly, this is well below the reported Km values of AADAT for α-aminoadipate (0.9 mM) and kynurenine (4.7 mM) 34 .
Given the observed effects of SLC17A4 and AADAT on T3 and T4 transport and metabolism, we additionally tested the associations of the identified genetic variants in SLC17A4 and AADAT with T3 levels and the T3/T4 ratio (Supplementary Table 6). SLC17A4-rs9356988 was associated with the T3/T4 ratio, while AADAT-rs6854291 was associated with both the T3/ T4 ratio and T3 levels.
Genetic TSH and FT4 risk score associations. To assess the cumulative clinical impact of our GWAS findings, we calculated a weighted GRS for TSH and FT4 levels, which included all independent TSH-and FT4-associated variants, respectively.
Next, these GRSs were tested for association with the risk of hypothyroidism and hyperthyroidism in up to 21,287 individuals. Figure 5 shows substantial differences in the risk of thyroid dysfunction across the range of GRS scores. Individuals with a TSH-based GRS in the highest quartile compared to individuals with a GRS in the lowest quartile had an odds ratio of 2.53 (p = 6.8 × 10 −32 ) for hypothyroidism and 0.19 (p = 9.8 × 10 −31 ) for hyperthyroidism, respectively. Conversely, the FT4-based GRS did not show any significant associations with either hypo or hyperthyroidism (Supplementary Table 7), which is consistent with the limited overlap observed between TSH and FT4 loci. A GRS using all TSH-associated variants showed a significant association with Graves' disease (p = 2.9 × 10 −5 ) that remained significant after excluding the PSORS1C1 variant (p = 2.5 × 10 −4 ), indicating a polygenic contribution of TSH-associated variants detected in the general population to Graves' disease.
As normal thyroid function is essential for the physiological function of virtually all human tissues, we tested if the TSH and FT4 GRSs were associated with a broader range of phenotypes by using available GWAS results of these phenotypes. These results are shown in Supplementary Table 8 with effects provided per increase in standard deviation of either TSH or FT4. A higher TSH GRS was associated with both a lower risk of Graves' disease (odds ratio (OR) = 0.64, p = 2.0 × 10 −5 ) and goiter (OR = 0.30, p = 3.9 × 10 −27 ), and lower thyroid volume (Δvol = −23%, p = 1.3 × 10 −37 ), whereas a higher FT4 GRS was associated with a higher risk of goiter (OR = 1.52, p = 7.9 × 10 −3 ) and higher thyroid volume (Δvol = 9%, p = 3.8 × 10 −3 ). In addition, a higher TSH GRS was associated with a lower risk of schizophrenia (OR  Table 9). These associations match clinical and epidemiological observations.

Discussion
With 8 million genetic variants tested in up to 72,167 individuals, we present results of the largest GWAS on thyroid function and dysfunction performed so far. We identified 109 significantly and independently associated genetic variants, doubling the number of loci known to regulate thyroid function, which explain a substantial part of the variation in these traits. Importantly, we detected associations between these variants and thyroid diseases as well as various clinical end points, and functionally Almost all previously identified TSH and FT4-associated SNPs were also genome-wide significantly associated with the respective trait in our analyses: the 20 TSH and 4 FT4 associations of the sex-combined GWAS of Porcu et al. 11 , the two additional TSH SNPs as well as the FT4 association revealed in Taylor et al. 12 , and 17 of the 21 genome-wide significantly TSHassociated SNPs identified by Gudmundsson et al. 13 The remaining loci of the latter study were discovered in our study via the mouse candidate analysis and were also associated with hypothyroidism (TPO), associated with FT4, hypo and hyperthyroidism (FOXE1), or had a p < 1 × 10 −6 (FOXE1, ELK3, SIVA1). Only two FT4-associated loci, LPCAT2/CAPNS2 and NETO1/FBXO15, identified in the sex-stratified analysis of Porcu et al. did not replicate in our sex-specific GWAS (p ≥ 0.01). While all but 2 of the 18 cohorts (the Old Order Amish and the Baltimore longitudinal study on Aging) of Porcu et al. as well as four of the seven cohorts (TwinsUK GWAS, SardiNIA, Val Borbera and BHS) of Taylor et al. were also included in our GWAS, we more than doubled the sample size for both traits, thereby significantly increasing our power to detect new loci.
Consistent with HPT-axis physiology, most TSH-associated variants acted on gene expression levels in thyroid tissue, while the FT4-associated variants had more widespread effects on multiple known thyroid hormone target tissues. In addition, four of the newly identified variants were either associated with the risk of TPOAb positivity or Graves' disease, suggesting an underlying autoimmune-related pathophysiology. All of these findings confirm that GWAS in the general population provides a valuable method to identify genes implicated in thyroid physiology and/or thyroid disease. Moreover, these insights can be successfully translated into experimental evidence, as illustrated by our in vitro studies on SLC17A4 and AADAT.
To investigate the combined effect of the thyroid hormone associated risk variants, we calculated a GRS. The GRS of TSHassociated variants was significantly associated with the risk of hypothyroidism and hyperthyroidism. For some FT4-associated hits, the lack of association with TSH levels can be explained on physiological grounds. For example, the identified SNP in DIO1 decreases the enzymatic activity of the protein, leading to less T4 to T3 conversion, resulting in higher T4 levels, but lower T3 levels, resulting in no net effect on feedback to the pituitary and therefore no effect on TSH levels. Similar hypotheses could be postulated for other loci involved in thyroid hormone metabolism, such as DIO3OS and AADAT. To assess whether these genetically estimated TSH levels are clinically relevant or merely reflect physiological inter-individual differences in TSH levels (i.e., HPT-axis setpoint), we tested the GRS against various clinical end points. These analyses showed significant associations with thyroid diseases, and also with altered lipid levels (total and LDL cholesterol) and height, which are both known to be affected by hypo and hyperthyroidism. For example, short stature is one of the key characteristics of patients suffering from congenital hypothyroidism. Interestingly, associations were also found with kidney function and schizophrenia, for which the causal relationships are less apparent. Thyroid hormone has been shown to influence kidney development and filtration function 37,38 . Likewise, rodent and human studies have shown that both hypo and hyperthyroidism lead to disrupted prenatal glial cell development, which is an important step in the development of schizophrenia 39 , while various psychiatric diseases including schizophrenia are also thought to influence thyroid function via central effects on the HPT axis 40 . Future studies are needed to clarify the mechanisms underlying these associations, as it is not possible to solve these potential inverse causal relationships solely with GWAS results. Irrespective of the direction of the effects, our results suggest that the presence of kidney dysfunction and psychiatric symptoms in patients with thyroid disease deserve attention. Given the substantial increase in number of TSH and FT4-associated variants, which explain a substantial part of the variation in these traits, future studies should start exploring the use of these markers to predict the individual HPT-axis setpoint. This predicted setpoint could be used to guide treatment of thyroid diseases, which is important as despite normalized TSH and FT4 levels, a substantial part of treated patients still have persistent hypo or hyperthyroid complaints, leading to a lower quality of life 41 . This could be due to the fact that the TSH and FT4 levels are normalized to within the population-based reference ranges, but still deviate from the patient's individual  setpoint. For this purpose, a GWAS on the TSH/FT4 ratio could prove to be a more sensitive method to identify more variants, which specifically affect the HPT-axis setpoint. When interpreting the results of our GWAS studies on increased and decreased TSH levels, it is important to realize that these studies were performed in population-based cohorts, and not in dedicated thyroid disease patient cohorts. Individuals on thyroid medication or a history of thyroid surgery were excluded, resulting in a relative overrepresentation of individuals with subclinical forms of thyroid dysfunction. The identified variants are therefore expected to be a mix of variants, which have been previously associated with hypo or hyperthyroidism (e.g., TPO, FOXE1, and ATXN2) and variants that lead to a TSH level, which is slightly above or below the population-based reference ranges. These latter effects can either reflect true mild thyroid dysfunction with increased risk of clinical consequences or merely reflect a deviation from the individual HPT-axis setpoint with no clinical consequences. While our GRS analyses suggest that carrying multiple risk alleles leads to an increased risk of overt thyroid dysfunction and related clinical consequences, the exact contribution of each individual variant needs to be clarified in future studies.
Thyroid hormone is importantly metabolized through enzymatic deiodination by DIO1-3, but also undergoes alternative metabolic reactions, including conjugation with sulfate or glucuronic acid and modification of the alanine side-chain. The latter includes the conversion of T3 and T4 to their respective pyruvic acid metabolites TK3 and TK4, which requires the oxidative deamination of their alanine side-chain 36 . TK3 and TK4 have been detected in urine and bile of rat injected with radio-labeled T3 and T4 42 . Although these and other studies suggested an important role for the liver and kidney in the formation of these pyruvic acid metabolites, the involved enzyme(s) had not been identified. Our functional analyses demonstrated that AADAT effectively catalyzes the transamination of T4 and in particular T3 to TK4 and TK3, respectively. Moreover, AADAT is highly expressed in the liver, gastrointestinal tract, and kidney in human 43 . Taken together, AADAT activity may thus be critical for the rate of thyroid hormone metabolism, which likely underlies the association of AADAT with circulating FT4. Although the specific impact of the associated variant on AADAT expression has not been assessed yet, our eQTL co-localization studies indicate that the index SNP decreases AADAT transcript levels in the thyroid, and this in turn leads to increased circulating FT4 levels.
The functional analyses further demonstrate that human SLC17A4 is able to transport both T4 and T3. The protein belongs to the solute carrier 17 family, whose members transport various organic anions, such as p-aminohippuric acid. Genetic variation in the SLC17A4 locus has been associated with the progression of elevated serum urate levels to gout 27,44 . According to the GTEx data resource and previously reported studies 27 , SLC17A4 is predominantly expressed in human small intestinal and colonic epithelial cells, pancreas, liver, and kidney cortex, which could imply a role for this transporter in the metabolic clearance and entero-hepatic cycle of thyroid hormone. Future studies should investigate the pharmacokinetic properties of SLC17A4, its relative contributions to thyroid hormone transport in various individual tissues, as well as the effects of the identified SLC17A4 variants on thyroid hormone transport.
The findings from our functional studies do not only provide new insights into thyroid hormone physiology, but may also have important clinical implications. Hypothyroidism is treated with levothyroxine (LT4), which is inexpensive and administered orally. In recent decades, various factors have been identified which help to determine LT4 dose, such as weight, gastrointestinal diseases, and interfering drugs 45,46 . Despite this knowledge, ineffective LT4 supplementation is still a major clinical problem, as 30-50% of patients are either under-or over-treated and therefore remain at risk for the symptoms and complications associated with thyroid dysfunction 45,46 . Therefore, the identification of SLC17A4 as a thyroid hormone transporter and AADAT as a thyroid hormone metabolizing enzyme provides new insights into thyroid hormone physiology and opens up a potential avenue for novel therapeutic targets or optimization of existing ones to improve the care of patients suffering from thyroid dysfunction. All genetic findings in our study were limited to common or lowfrequency SNPs, whereas rare SNPs or structural variants may also contribute to the yet unexplained variance of thyroid function. Large whole-exome or -genome sequencing studies are required to reveal these rare variant associations 47 . Furthermore, additional GWAS with increased sample size will help to reveal the yet undiscovered associations of common and low-frequency SNPs. Our ThyroidOmics Consortium (http://www.thyroidomics.com) provides a well-established infrastructure to address these knowledge gaps in future projects.

Methods
Included studies. Discovery meta-analyses included data from 22 independent cohorts with 54,288 subjects for the TSH analyses, and from 19 cohorts with 49,269 subjects for FT4, 53,423 subjects (3440 cases) for hypothyroidism, and 51,823 subjects (1840 cases) for hyperthyroidism (Supplementary Data 1). Selected SNPs from the TSH or FT4 analyses were carried forward for replication with in silico GWAS data from 5 cohorts (9053 subjects) and de novo genotyping in additional 5 cohorts (13,330 subjects). All subjects gave informed consent and studies were approved by the cohort-specific ethics committees.
We used the results of the GWAS of TPOAb positivity that included 18,297 subjects 20 for a look-up of all the 53 TSH-associated loci or their HapMapII proxies (r² > 0.8 in a 1 Mb window) that were available in that dataset to assess their relation to autoimmune hypothyroidism. A complementary look-up was performed for the 52 SNPs that were available in a GWAS on Graves' disease diagnosed by clinical examinations, circulating thyroid hormone and TSH concentrations, serum levels of antibodies against thyroglobulin, thyroid microsomes, and TSH receptors, ultrasonography, [99m] TCO 4 − (technetium-99m pertechnetate) (or [ 123 I] (radioactive iodine)) uptake and thyroid scintigraphy using the data of the BioBank Japan Project (BBJ) including 1747 patients and 6420 controls (Supplementary Data 1).
Trait definition. In each study, only subjects with TSH levels within the cohortspecific reference range were included for the TSH and FT4 analyses. TSH and FT4 were analyzed as continuous variables after inverse normal transformation. Increased TSH was defined by a TSH level above the upper limit of the cohortspecific TSH reference range, while decreased TSH was defined by a level below the lower limit of the reference range. For both increased and decreased TSH analyses, the comparison group consisted of subject with a TSH level within the cohortspecific reference range. Exclusion criteria for all analyses were non-European ancestry, use of thyroid medication (defined as ATC (Anatomical Therapeutic Chemical) code H03), or previous thyroid surgery.
GWAS in individual studies. In each study of the discovery GWAS, genotyping was performed on genome-wide arrays. Genome-wide data were imputed to the 1000 Genomes, phase 1 version 3 (March 2012) ALL populations reference panel, including the X chromosome. Quality control before imputation was applied in each study separately. Details on study-specific genotyping and imputation information are provided in the Supplementary Data 1.
In the individual study GWAS, the association of the SNPs was analyzed using linear regression for TSH and FT4, and logistic regression for decreased and increased TSH. The genotype-phenotype association was conducted using an additive genetic model on SNP dosages, thus taking genotype uncertainties of imputed SNPs into account. The analyses for TSH and FT4 were initially sexstratified and meta-analyzed as a second step. The analyses were adjusted for age, age-squared (to account for non-linear effects), and relevant study-specific covariates such as principal components for population stratification, study center, and family structure (e.g., by inclusion of the kinship matrix as a random effect), if applicable. The family-based cohorts GARP, SardiNIA, ValBorbera, MICROS, TwinsUK, LLS, and FHS conducted additional analyses on the men and womencombined sample, with additional adjustment for sex, to properly account for their family relatedness.
Statistical methods for meta-analysis. Result files from individual studies included in this analysis underwent extensive quality control before meta-analysis: file format checks as well as plausibility and distributions of association results including effects, standard errors, allele frequencies, and imputation quality of the SNPs were obtained by using the gwasqc() function of the GWAtoolbox package v2.2.4 48 . Additionally, the known associations of rs6885099 in PDE8B with TSH and rs2235544 in DIO1 with FT4 were checked for consistent effect direction and size in each study. All cohort-specific genomic control values (λ GC ) ranged from 0.94 to 1.14 (median 1.00) for the continuous trait and from 0.68 to 1.04 (median 0.91) for the dichotomous trait GWAS.
All meta-analyses were carried out in duplicate by three independent analysts. We conducted a fixed-effect meta-analysis applying inverse-variance weighting as implemented in Metal 49 . SNPs with MAF ≤0.005 or an imputation quality score ≤0.4 were excluded prior to the meta-analyses resulting in a median of 9,653,808 SNPs per cohort (IQR: 9,302,604-10,705,092). During the meta-analysis, the studyspecific results were corrected by their specific λ GC if >1. Results were checked for possible errors like use of incorrect unit, trait transformation, or association model by plotting the association p-values of the analyses against those from a z-scorebased meta-analysis for verifying overall concordance. SNPs that were present in <75% of the total sample size contributing to the respective meta-analysis (separately for autosomal and X-chromosomal SNPs) or with a MAF ≤0.01 (hypo and hyperthyroidism MAF ≤0.05 because of the low number of cases in the analysis) were excluded from subsequent analyses. Finally, data for up to 8,048,941 genotyped or imputed autosomal and X-chromosomal SNPs were available after the discovery stage meta-analysis of TSH, FT4, and up to 5,965,951 SNPs after hypo and hyperthyroidism.
Quantile-quantile plots of the meta-analysis results are provided in Supplementary Figures 10 and 11. To assess whether there was an inflation of pvalues in the meta-analysis results attributed to reasons other than polygenicity, we performed LD score regression 50 . The LD score-corrected λ GC values of the metaanalysis results ranged from 1.00 to 1.04, supporting the absence of unaccounted population stratification. Genome-wide significance was defined as a p-value of <5 × 10 −8 , corresponding to a Bonferroni correction of one million independent tests. Unless stated otherwise, all reported p-values are two-sided. The I 2 statistic was used to evaluate between-study heterogeneity 51 .
Gene-by-sex interaction on the circulating TSH and FT4 levels were obtained for each SNP by comparing the discovery meta-analysis results from men (TSH: n = 24,618; FT4 n = 22,315) and women using a t-test. Test statistics were calculated using the formula t = (β men − β women )/sqrt(SE men ² + SE women ²), assuming independent effect sizes between men and women.
To evaluate the presence of independent SNPs within each locus, SNPs were clustered based on their correlation with the SNP showing the lowest p-value at that locus (index SNP) using the software PLINK 52 and the genotypes of the combined individuals of the 1000Genomes phase1v3 all ethnicities reference panel (linkage disequilibrium pruning using r 2 ≤ 0.01 within windows of ±1 Mb). The loci were named according to the nearest gene of the index SNP. Genomic positions correspond to build 37 (GRCh37).
Replication analysis. The genome-wide significant index SNPs of newly identified loci from the sex-combined TSH (n = 22) and FT4 meta-analyses (n = 19) were taken forward to the replication stage (Supplementary Table 10). When SNPs were not available in the in silico replication datasets, a proxy SNP in LD with r² > 0.8 was selected.
Of the ten studies that contributed to replication, five studies used 1000Genomes imputed dosages, three studies performed de novo genotyping, and two studies were genotyped on both the Illumina ExomeChip and CardiometaboChip. No SNP or proxy for the X-chromosomal locus was available in any replication dataset. The results from the discovery meta-analysis and the results of replication studies were meta-analyzed to obtain the overall p-values of the selected SNPs. SNPs with p-values below genome-wide significance in this combined analysis and with concordant effect directions in both stages were considered as replicated 53 .
Integration of information from genetically manipulated mice. We tested whether information about thyroid function or disease from genetically manipulated mice could facilitate the detection of additional human thyroid loci that did not reach genome-wide significance in the GWAS (nested candidate gene approach). To this end, all genes that when manipulated cause abnormal thyroid physiology (MP:0002876; 26 genes) or abnormal thyroid gland morphology (MP:0000681; 51 genes) were selected from the comprehensive Mouse Genome Informatics resource in October of 2016 (http://www.informatics.jax.org/mp/). Next, genes were translated to their human homologs, followed by the calculation of the number of independent SNPs in these genes with MAF > 0.01 in the 1000 Genomes EUR populations (PLINK option-indep-pairwise 50 5 0.2) to obtain multiple testing corrected significance thresholds (p < 2.5 × 10 −5 for physiology and p < 1.9 × 10 −5 for morphology). The genes in the respective mouse lists were then queried for the presence of SNPs that showed association with TSH and/or FT4 below the significance threshold. To test whether the number of genes with significant associations was higher than expected by chance, results were compared to those from 2000 iterations of random gene lists of equal length. A p-value for enrichment was computed from a complementary cumulative binomial distribution as described in detail previously 54 . Lastly, novel loci that were not identified at genome-wide significance in the GWAS of TSH or FT4 were tested for replication in up to 9011 and 4532 independent samples for TSH and FT4, respectively.
Successful replication was defined as direction-consistent association and genome-wide significance in a meta-analysis of the discovery and replication samples.
eQTL look-up. To assess the possible effect of our lead signals on transcriptional activity, we queried expression QTL (eQTL) results from 22 publicly available studies (specific reference listed in Supplementary Data 3). These studies were carried out from 2007 to February 2017 on 127 different tissues and cell types, and used either micro-array or sequencing-based assessment of gene expression. For each study, we derived the list of top eQTLs by LD clumping, and searched top eQTLs in high LD (r 2 > 0.8 in 1000Genomes EUR samples) with the 94 thyroid function-associated index or independent SNPs (TSH, FT4, and ATXN2 and TPO as additional GWAS loci from hypothyroidism) (Supplementary Tables 1, 2,  and 4).
To evaluate the evidence of co-localization between the index GWAS and eQTL SNPs, we used the SMR method 24 , coupled with the test for heterogeneity of effects (HEIDI) 24 . The first tests whether the effect on expression seen at a SNP or at its proxies correlates with the signal observed in the GWAS (SMR test), while the second evaluates if the eQTL and GWAS associations can be attributable to the same causative variant (HEIDI test). Because direction of effects has to be taken into account, we focused this analysis only on GTEx data for which full summary results were available. For SMR, we considered the experiment-wise p-value of 2 × 10 −4 (corresponding to a Bonferroni correction for 242 gene-thyroid trait-tissue combinations assessed). Specifically, we tested all genes with an eQTL p-value <1 × 10 −7 and for which the top eQTL showed genome-wide significant association with any thyroid hormone traits, regardless of LD between the top eQTL and the thyroid hormone index SNP. For the HEIDI test, we used the suggested p-value >0.05 cutoff to declare co-localization, and further required that at least five SNPs were available for the test 24 .
Expression constructs and cloning. The cDNA of MCT8 and CRYM was cloned into pcDNA3 and pSG5 expression vectors, respectively, using standard cloning techniques 56,57 . A pCMV6-Entry_SLC17A4 expression vector containing a Cterminal Myc and Flag tag was obtained from OriGene Technologies (Rockville, USA). A pbluescript AADAT cDNA construct was obtained from Thermo Scientific (Bleiswijk, NL) and subcloned into pcDNA3 with addition of a C-terminal Flag-tag using standard cloning techniques. Any variants were substituted in agreement with the NM_001286683.1 reference sequence using Quikchange sitedirected mutagenesis according to manufacturer's protocol (Stratagene, Amsterdam, The Netherlands). All primers are available upon request. Correctness of all expression constructs was confirmed by sequencing of the inserts.
For T3 and T4 uptake studies, COS-1 cells were seeded in 24-well dishes (1 × 10 5 cells per well) and transiently transfected at 70% confluence with 250 ng empty vector (EV), SLC17A4, or MCT8, with or without the addition of 100 ng CRYM. CRYM is a cytoplasmic high-affinity thyroid hormone-binding protein, which prevents efflux of internalized thyroid hormones. For thyroid hormone metabolism assays, COS-1 cells were seeded on 10 cm dishes (5 × 10 5 cells per well) and transiently transfected with 2000 ng EV or AADAT at 70% confluence. Xtreme-Gene 9 was used as a transfection reagent according to manufacturer's protocol (Roche Diagnostics, Almere, NL).
Genetic risk score analysis. Two separate GWAS effect size-weighted GRS were generated to evaluate the combined effect of the TSH-and FT4increasing alleles, respectively, on the risk of hypo and hyperthyroidism using individual level data from four of our largest GWAS studies: ARIC, CHS, Rotterdam Study, and SHIP (hypothyroidism cases: n = 1613; hyperthyroidism cases: n = 662; controls: n = 19,674). The GRS were based on the 61 and 31 replicated GWAS SNPs for TSH and FT4, respectively (Supplementary Tables 1  and 2), normalized to a range of 0 to 100, associated in each cohort separately using a logistic regression adjusted for sex and age, and combined afterwards by a fixed-effect inverse-variance meta-analysis using R 60 . The probability of disease was calculated using the formula 1/(1 + exp(−(β 0 + β 1 *x))), where β 0 and β 1 correspond to the intercept and GRS-related effect in the regression model, respectively.
To test the combined effect of the replicated TSH and FT4 SNPs on various traits (Supplementary Tables 8 and 9), a GRS-based association on meta-analyses results was performed as described in reference 61 using the function grs.summary() of the R-package gtx. If a specific SNP was not available in the look-up GWAS dataset, a proxy SNP in LD with r² > 0.8 was included where possible. The association effects correspond to a change in the look-up GWAS trait (or natural logarithm of the odds ratio in the case of a binary trait) per standard deviation unit of TSH and FT4, respectively. In case the trait was logarithm-transformed in the GWAS, (e beta −1) × 100 corresponds to a 1% change of this trait.
Pathway analyses. We performed Data-Driven Expression Prioritized Integration for Complex Traits (DEPICT) 62 analyses to prioritize the most likely causal genes at the associated loci and identify enriched pathways and tissues. We used DEPICT version 1 rel194 and included variants with GC-corrected p-value <1 × 10 −5 from discovery GWAS as input. The following parameters were applied in the DEPICT analyses: 50 repetitions used to compute the false discovery rate, 500 permutations used to adjust for biases such as gene length, and 500 null GWAS used to run repetitions and permutations.
Genes for network analysis were selected using the associated genes from DEPICT using the lead SNPs from the analyses of the discovery GWAS. The Ingenuity Pathway Analysis Software Tool (IPA; Ingenuity® Systems, CA, USA) Network was used in order to perform pathway analysis (core analysis). Molecules and/or relationships considered were the ones available in the IPA Knowledge Base for mammals. Confidence filters considered only relationships, direct and indirect, where the confidence is experimentally observed or high (predicted). Networks were generated with a maximum size of 35 genes and allowing up to 25 networks per analysis. We did not restrict for tissue and cell lines or mutations. The networks are constructed using the IPA algorithm with the Ingenuity Knowledge Base as a reference set generating a score as well as a p-value. IPA computes a score for each network according to the fit of that network to the user-defined set of focus genes. The score indicates the likelihood of the focus genes in a network being found together due to random chance. The significance of p-value is calculated using the right-tailed fisher exact test.
Look-up of pleiotropy. The look-up of additional traits associated with the replicated GWAS findings were performed using the PhenoScanner 63 database with the default search options, whereas the independent SNPs of the replicated TSH, FT4, hypo, and hyperthyroidism-associated loci or their proxies (r² > 0.8) obtained by SNiPa 64 were used as input. Look-up results with genome-wide significant p-value (5 × 10 −8 ) were reported.