A Meta-Analysis of Thyroid-Related Traits Reveals Novel Loci and Gender-Specific Differences in the Regulation of Thyroid Function

Thyroid hormone is essential for normal metabolism and development, and overt abnormalities in thyroid function lead to common endocrine disorders affecting approximately 10% of individuals over their life span. In addition, even mild alterations in thyroid function are associated with weight changes, atrial fibrillation, osteoporosis, and psychiatric disorders. To identify novel variants underlying thyroid function, we performed a large meta-analysis of genome-wide association studies for serum levels of the highly heritable thyroid function markers TSH and FT4, in up to 26,420 and 17,520 euthyroid subjects, respectively. Here we report 26 independent associations, including several novel loci for TSH (PDE10A, VEGFA, IGFBP5, NFIA, SOX9, PRDM11, FGF7, INSR, ABO, MIR1179, NRG1, MBIP, ITPK1, SASH1, GLIS3) and FT4 (LHX3, FOXE1, AADAT, NETO1/FBXO15, LPCAT2/CAPNS2). Notably, only limited overlap was detected between TSH and FT4 associated signals, in spite of the feedback regulation of their circulating levels by the hypothalamic-pituitary-thyroid axis. Five of the reported loci (PDE8B, PDE10A, MAF/LOC440389, NETO1/FBXO15, and LPCAT2/CAPNS2) show strong gender-specific differences, which offer clues for the known sexual dimorphism in thyroid function and related pathologies. Importantly, the TSH-associated loci contribute not only to variation within the normal range, but also to TSH values outside the reference range, suggesting that they may be involved in thyroid dysfunction. Overall, our findings explain, respectively, 5.64% and 2.30% of total TSH and FT4 trait variance, and they improve the current knowledge of the regulation of hypothalamic-pituitary-thyroid axis function and the consequences of genetic variation for hypo- or hyperthyroidism.


Introduction
Through the production of thyroid hormone (TH), the thyroid is essential for normal development, growth and metabolism of virtually all human tissues. Its critical role in heart, brain, bone, and general metabolism is illustrated by the clinical manifestations of thyroid disease, which affects up to 10% of the population. Low thyroid function (i.e., hypothyroidism) can lead to weight gain, high cholesterol, cognitive dysfunction, depression, and cold intolerance, whereas hyperthyroidism may result in weight loss, tachycardia, atrial fibrillation, and osteoporosis. Mild variation in thyroid function, both subclinical and within the normal range, is associated with these TH-related clinical outcomes as well [1][2][3][4].
The thyroid gland secretes predominantly the pro-hormone thyroxine (T4), which is converted into the active form triiodothyronine (T3) in peripheral tissues. The production of TH by the thyroid gland is regulated by the hypothalamus-pituitary-thyroid (HPT) axis, via a so-called negative feedback loop. Briefly, low levels of serum TH in hypothyroidism result in an increased release of thyroid stimulating hormone (TSH) by the pituitary, under the influence of hypothalamic thyrotropin releasing hormone (TRH) [5]. TSH, a key regulator of thyroid function, stimulates the synthesis and secretion of TH by the thyroid. When circulating TH levels are high, as in hyperthyroidism, TRH and TSH synthesis and secretion are inhibited.
In healthy (euthyroid) individuals, TSH and free T4 (FT4) levels vary over a narrower range than the broad inter-individual variation seen in the general population, suggesting that each person has a unique HPT axis set-point that lies within the population reference range [6]. Besides environmental factors such as diet, smoking and medication, little is known about the factors that influence this inter-individual variation in TSH and FT4 levels [7][8][9]. The heritability of TSH and FT4 has been estimated from twin and family studies at about 65% and 40%, respectively [10][11][12]. However, the underlying genetic variants are not fully established, and the contribution of those discovered so far to the overall variance is modest. Single nucleotide polymorphisms (SNPs) in the phosphodiesterase type 8B (PDE8B), upstream of the capping protein (actin filament) muscle Z-line, b (CAPZB) and, more recently, of the nuclear receptor subfamily 3, group C, member 2 (NR3C2) and of v-maf musculoaponeurotic fibrosarcoma oncogene homolog (MAF/ LOC440389) genes have been implicated in TSH variation by genome-wide association studies (GWAS) [13][14][15], whereas SNPs in the iodothyronine deiodinase DIO1 have been associated with circulating levels of TH by candidate gene analysis [16][17][18].
To identify additional common variants associated with thyroid function, we performed a meta-analysis of genome-wide association data in 26,420 euthyroid individuals phenotyped for serum TSH and 17,520 for FT4 levels, respectively. In addition, we also assessed gender-specific effects and correlation with subclinical thyroid dysfunction.

Results
To identify common genetic variants associated with serum TSH and FT4 levels, we carried out a meta-analysis of genome-wide association results from 18 studies for TSH and 15 studies for FT4 levels, which assessed the additive effect of ,2.5 million genotyped and HapMap-imputed SNPs in relation to those traits in individuals of European ancestry (for cohort description see Table 1 and Table S1). In order to avoid bias due to the presence of thyroid pathologies, prior to analysis we excluded all individuals with TSH values outside the normal range (TSH,0.4 mIU/L and TSH.4.0 mIU/L) and those taking thyroid medication for known thyroid pathologies whenever the relevant information was available. Our meta-analysis was thereby carried out in up to 26,420 and 17,520 euthyroid subjects, respectively for TSH and FT4. Additional exclusion criteria used by individual cohorts are detailed in Table S1.
Using the standard genome-wide threshold of 5610 28 , we observed significant associations for SNPs at 23 loci, of which 19 were associated with TSH, and 4 with FT4 ( Figure S1). The results are presented in Table 2 and Figure 1, Figure 2, Figure 3, Figure 4, Figure 5. In Table S2 single cohort results for each GW significant SNP are reported.
At each locus, a single variant was sufficient to explain entirely the observed association, except for the VEGFA locus, which contained an independent signal located 150 kb downstream of the gene, detected by conditional analyses ( Figure 1F and Table 2).
Of all 24 independent markers, significant evidence for heterogeneity (P,0.002, corresponding to a Bonferroni threshold of 0.5/24) was only observed at ABO (P = 1.22610 24 ). Iodine nutrition, which may profoundly affect thyroid function, is quite different in some of the cohorts under study (i.e., Europe vs North America). To test whether the observed heterogeneity could be attributable to different iodine intake, we combined cohorts from South Europe (an iodine-deficient region) and compared effect sizes with those observed in a meta-analysis of North American samples (an iodine-replete region). Interestingly, the effect size of the top marker at ABO was three times larger in Europeans vs North American, and this difference remained significant after Bonferroni correction (P = 7.0.9610 24 ) (Table S3). However, the relation of the ABO SNP, a tag for the blood group O, to iodine intake remains to be determined.

Gender-specific analyses
Given the reported clinical differences in thyroid function in males and females [21][22][23], we searched for gender-specific loci by whole-genome sex-specific meta-analysis, analyzing males and females separately in each cohort. Some of the loci detected in the main meta-analysis were seen at genome-wide significance level only in females (NR3C2, VEGFA, NRG1 and SASH1) or in males (MAF/LOC440389, FGF7, SOX9, IGFBP5) with either the same top SNP or one surrogate, but effect sizes at their variants were significantly gender-specific only at PDE8B, PDE10A and MAF/ LOC440389, considering a false discovery rate of 5% [24]. In addition, effects at MAF/LOC440389 were significantly different also at the more stringent Bonferroni threshold of 1.9610 23 ( = 0.05/26), and close to significance at PDE8B and PDE10A (Table 3). At these latter loci, the TSH-elevating alleles showed a stronger impact on trait variability in males compared to females ( Figure 6). In addition, the gender specific meta-analysis for FT4, revealed a novel female-specific locus on chromosome 18q22, and a novel male-specific locus on chromosome 16q12.2, that had not been detected in the main meta-analysis (Table 3, Figure 6 and Figure S2). The female-specific signal (rs7240777, P = 3.49610 28 ) maps in a ''gene desert'' region, with the nearest genes NETO1 (neuropilin (NRP) and tolloid (TLL)-like 1), located, about 550 kb upstream and FBXO15 (F-box only protein 15) 500 kb downstream ( Figure 5D). The male-specific association is located in intron 11 of the LPCAT2 (lysophosphatidylcholine acyltransferase 2) gene, and near CAPNS2 (calpain, small subunit 2) (rs6499766, P = 4.63610 28 ), a gene which may play a role in spermatogenesis [25]. The FT4elevating alleles in the NETO1/FBXO15 and LPCAT2/CAPNS2 were fully gender-specific, i.e. there was no effect in males and in females, respectively (P.0.01).
Overall, the 20 TSH and the 6 FT4 associations account, respectively, for 5.64% and 2.30% of total trait variance.

Common loci regulating TSH and FT4 levels
To explore overlap between TSH-and FT4-associated loci and their involvement in the HPT-negative feedback loop, we assessed the associations of the top TSH-associated SNPs on FT4 levels, and vice versa. For the SNPs in or near PDE8B, MAF/LOC440389, VEGFA, IGFBP5, NFIA, MIR1179, MBIP and GLIS3 the TSHelevating allele appeared to be associated with decreasing FT4 levels (P,0.05, Table S4). However, after application of Bonferroni correction (threshold for FT4 association of TSH SNPs, P = 2.5610 23 ), none of these reciprocal associations remained significant.
By contrast, a positive relationship was seen for one of the FT4 associated loci, since the variant at the LHX3 locus was significantly associated with higher levels of both FT4 and TSH (P = 5.25610 23 , with Bonferroni threshold 0.05/6 = 0.008).
As the presence of reciprocal associations between TSH and FT4 regulating SNPs would be expected from physiology, we tested the power of our study to detect such a relationship. Power calculation for the top SNP at PDE8B, which has the largest effect on TSH levels, revealed that our meta-analysis only has 9% power to detect an association of FT4 at a Bonferroni P = 2.5610 23 . We also carried out a bivariate analysis in the SardiNIA study using poly software to estimate specific contributions [26]. This analysis showed that most of the observed negative feedback correlation is due to environmental factors (environmental correlation = 20.130, genetic correlation = 20.065).

Association of loci with hypothyroidism and hyperthyroidism
To assess possible clinical implications, we investigated whether the variants identified in individuals without overt thyroid pathologies (i.e., with TSH levels within the normal range and not taking thyroid medication) were also associated in individuals with abnormal TSH values (i.e., outside the reference range), who were not included in the initial meta-analysis as potentially affected by thyroid pathology. Towards this, we first assessed the global impact of TSH-and FT4-associated SNPs on the risk of increased or decreased TSH levels by comparing weighted genotype risk score (GRS) quartiles in the individuals with abnormal TSH values that were discarded for the GWAS analyses. For the TSHassociated SNPs, the odds of increased TSH levels were 6.65 times greater in individuals with a GRS in the top quartile compared to individuals in the bottom quartile (P = 3.43610 220 ) (

Author Summary
Levels of thyroid hormones are tightly regulated by TSH produced in the pituitary, and even mild alterations in their concentrations are strong indicators of thyroid pathologies, which are very common worldwide. To identify common genetic variants associated with the highly heritable markers of thyroid function, TSH and FT4, we conducted a meta-analysis of genome-wide association studies in 26,420 and 17,520 individuals, respectively, of European ancestry with normal thyroid function. Our analysis identified 26 independent genetic variants regulating these traits, several of which are new, and confirmed previously detected polymorphisms affecting TSH (within the PDE8B gene and near CAPZB, MAF/LOC440389, and NR3C2) and FT4 (within DIO1) levels. Gender-specific differences in the genetic effects of several variants for TSH and FT4 levels were identified at several loci, which offer clues to understand the known sexual dimorphism in thyroid function and pathology. Of particular clinical interest, we show that TSH-associated loci contribute not only to normal variation, but also to TSH values outside reference range, suggesting that they may be involved in thyroid dysfunction. Overall, our findings add to the developing landscape of the regulation of thyroid homeostasis and the consequences of genetic variation for thyroid related diseases.
FT4-associated SNPs we found no significant associations for any of the tested comparisons (data not shown).
We also assessed the 20 independent TSH SNPs individually in relation to the risk of abnormal TSH levels by case-control metaanalysis in subjects with high (cases) versus low (controls) TSH values. This analysis showed that variants at PDE8B, CAPZB, FGF7, PDE10A, NFIA and ITPK1 loci are significantly associated (Bonferroni threshold P = 2.5610 23 ) with abnormal TSH levels ( Table 4, bottom panel). PDE8B, CAPZB and FGF7 were also strongly associated with the risk of decreased TSH levels in an analysis of individuals with low (cases) versus normal range TSH (controls). In addition, variants at VEGFA were also significantly associated in this comparison. Finally, when individuals with high TSH values were analyzed versus controls, the NR3C2 locus appeared significantly associated in addition to PDE8B and CAPZB.

Association of TSH lead SNPs in pregnant women
Normal thyroid function is particularly important during pregnancy and elevated TSH levels are implicated in a number of adverse outcomes for both mother and offspring. We therefore assessed whether the TSH lead SNPs were also associated with elevated TSH during pregnancy, when increased TH production is necessary. We tested 9 of the 20 lead TSH variants (or their proxies, see Text S1) in a cohort of 974 healthy pregnant women at 28 weeks gestation [27] and found, as expected, that mean TSH levels were correlated with the number of TSH-elevating alleles (P = 3.0610 212 , Table S5). Effect size estimates in pregnant women were not significantly different when compared to those of women in the main gender-specific meta-analysis (heterogeneity P value.0.05), suggesting that the effects of the TSH-elevating alleles are no greater during pregnancy (data not shown). However, there was evidence of association between the number of TSH-raising alleles and subclinical hypothyroidism in pregnan-cy, both in the whole sample (OR per weighted allele: 1.18 [95%CI: 1.01, 1.37], P = 0.04) and in TPO antibody-negative women (1.29 [95%CI: 1.08, 1.55], P = 0.006) (Table S6).

Discussion
We report 26 independent SNPs associated with thyroid function tests in euthyroid subjects, 21 of which represent novel signals (16 for TSH and 5 for FT4). Overall they explain 5.64% and 2.30% of the variation in TSH and FT4 levels, respectively.
We observed that carriers of multiple TSH-elevating alleles have increased risk of abnormal TSH levels, and also found association between the number of TSH-elevating alleles and subclinical hypothyroidism in pregnancy. These results are potentially clinically relevant, because abnormal TSH values are the most sensitive diagnostic markers for both overt and subclinical thyroid disease [4]. The variants identified in the current study, or those in LD with them, may thus contribute to the pathogenesis of thyroid disease. Of note, we found eight loci significantly associated with abnormal TSH levels (PDE8B, PDE10, CAPZB, VEGFA, NR3C2, FGF7, NFIA and ITPK1), of which two were specifically associated with either abnormally low (VEGFA) or elevated (NR3C2) TSH values, suggesting differential mechanisms for the contribution of these variants to hyper-and hypothyroidism, respectively. Interestingly, the mineralocorticoid receptor NR3C2 gene has recently been found to be up-regulated in adultonset hypothyroidism [28], and PDE8B and CAPZB have been suggestively associated with hypothyroidism by GWAS [29]. Alternatively, it may be that carriers of these alleles are healthy individuals who may be misdiagnosed as having thyroid disease because their genetically determined TSH concentrations fall outside the population-based reference range. More research is required to determine which of these interpretations is correct, and the relevance of these variants as markers for thyroid dysfunction or thyroid-related clinical endpoints. The evidence for gender-specific differences at several TSH and FT4 regulatory loci is intriguing. They included variants at PDE8B, PDE10A, and MAF/LOC440389, which showed significantly stronger genetic effects with pituitary-thyroid function in males, and variants at NETO1/FBX015 and LPCAT2/CAPNS2 which seems to have an effect only in females and males, respectively. Sex differences in the regulation of thyroid function have generally been linked to the influence of sex hormones and autoimmune thyroid disease, resulting in a higher prevalence of thyroid dysfunction in women, without clear understanding of underlying molecular mechanisms [21][22][23]. Our study suggests that differential genes and mechanisms are potentially implicated in the regulation of thyroid function in men and women. Given the impact of thyroid function on several disease outcomes as well as male and female fertility and reproduction, clarifying the underlying associations may provide additional insight for future interventions.
Although it is well known that TSH and FT4 levels are tightly regulated through a negative feedback loop involving the HPT axis, we detected significant overlap between TSH and FT4 signals only at the LHX3 locus, which was primarily associated in our study with FT4. The LHX3 allele is associated with an increase of both TSH and FT4, which is consistent with the essential role of this transcription factor in pituitary development. Inactivating mutations in LHX3 cause the combined pituitary hormone deficiency-3 syndrome [CPHD3 (MIM#221750)] [30,31], characterized by low TSH and FT4 levels. The positive association of the LHX3 variant with both TSH and FT4 suggests an effect of this allele at the level of the HPT-axis, resulting in an increased exposure to thyroid hormone throughout life. In contrast, although several of the TSH-elevating alleles appeared to be associated with decreasing FT4 levels, none of these reciprocal associations remained significant after Bonferroni correction. Lack of loci associated in a reciprocal manner with both TSH and FT4 is somewhat puzzling, as their presence would be expected from physiology. However, these findings are consistent with initial reports by Shields et al. [27] and more The table shows the association results for SNPs that reached genome-wide level (p,5610 208 ) in the main meta-analysis. SNPs at LPCAT2/CAPNS2 and NETO1/FBXO15 reached the GW threshold in the gender-specific meta analysis (see Table 3), and here the p-value in the main meta-analysis is reported. For each SNP, the best candidate gene is showed, as well as its genomic position in build 36, the effect allele (A1) and the other allele (A2), its combined frequency across studies and its standard error, the effect size and its standard error, the p-value for association, the number of samples analyzed, and the p-values for heterogeneity of effects across the cohorts meta-analyzed. Effect sizes are standardized, so they represent the estimated phenotypic change, per each copy of the effect allele, in standard deviation units. doi:10.1371/journal.pgen.1003266.t002 recent findings by Gudmundsson et al. [32]. A power analysis showed that our study -in spite of being one of the larger conducted so far on these traits -is underpowered to detect an inverse relationship between TSH and FT4 variants, considering a Spearman rank correlation of 20.130 between these traits [12]. As a consequence, contrasting studies on smaller sample sizes may also lack power and cannot be considered robust when testing this relationship [33]. In addition, we estimated that most of the observed negative feedback correlation is due to environmental factors; so it is unlikely that negative feedback is controlled by a genetic locus with large effect. This observation can rationalize the lack of reciprocal, significant associations detected for both TSH and FT4 in this and other studies, and further supports the crucial role of the HPT-axis in maintaining normal levels of thyroid hormone. At present the relationship between the associated variants and specific mechanisms involved in regulating TSH and FT4 levels has not been established, but we have identified strong candidates at the majority of the loci by literature-mining approaches, as detailed below and in Table 5.
Most of the 16 novel loci implicated in the regulation of TSH are highly represented in the thyroid with the exception of PRDM11, expressed in brain, ABO, in blood, and MIR1179. PDE10A encodes a cAMP-stimulated phosphodiesterase, which was previously only suggestively associated with TSH levels and hypothyroidism [13,34], although the tested variants were weakly correlated with our top signal (r 2 = 0.55 with rs2983521 and r 2 = 0.15 with rs9347083). The presence of linkage at this gene in families reaching accepted clinical criteria of thyroid dysfunction reinforces the observation that variants in this gene may contribute to clinical thyroid disorders [34]. PDE10A, together with PDE8B and CAPZB, emerged in our study as the strongest currently known genetic determinants of this trait. Both PDE8B and PDE10A are implicated in cAMP degradation in response to TSH stimulation of thyrocytes. In addition, the activity of both PDE10A and CAPZB appear modulated by cAMP [35,36]. These  Table 5), but further studies will be necessary to clarify the exact biological mechanisms and the specific genes involved at each locus. The association of TSH levels with IGFBP5, INSR and NR3C2 is, however, an indication of a specific role of the growth hormone/insulin-like growth factor (GH/IGF) pathway in thyroid function. Remarkably, expression of IGFBP5 is tightly regulated by cAMP, again underlying the pivotal role of this second messenger in determining net TSH levels [37]. For FT4, the DIO1, FOXE1 and LHX3 identified loci have strong biological support as potential effectors. While both DIO1 and FOXE1 were previously associated with FT4 levels and hypothyroidism by candidate gene analysis and functional studies [17][18][19][38][39][40][41], association at LHX3 is novel and is consistent with the essential role of this transcription factor in pituitary development (see above) [30,31,42,43]. Consistent with the role of pituitary in growth, this locus has also recently been associated with height in Japanese [44]. The associations of AADAT, NETO1/ FXBO15 and LPCT2/CAPNS2 with FT4 levels are currently less clear. It may be relevant that AADAT catalyzes the synthesis of kynurenic acid (KYNA) from kynurenine (KYN), a pathway that has been associated with the induction in brain of proinflammatory cytokines that are known to activate the hypothalamopituitary-adrenal (HPA) axis, in turn affecting the HPT axis and thyroid function, including FT4 levels [45][46][47][48][49].
Additional pathway analyses by MAGENTA [50], GRAIL [51], and IPA (Ingenuity Systems, www.ingenuity.com) to look for functional enrichment of the genes mapping to the regions associated with TSH, FT4 or both, yielded no novel interactions. However, IPA highlighted an over-representation of genes implicated in developmental processes (11/26, P = 6.27610 26 -8.85610 23 ) and cancer (16/26 loci, P = 2.44610 26 -9.30610 23 ). This is consistent with the notion that a normally developed thyroid gland is essential for both proper function and thyroid hormone synthesis, and that defects in any of the essential steps in thyroid development or thyroid hormone synthesis may result in morphologic abnormalities, impaired hormonogenesis and growth dysregulation. It is also interesting to note that 11 of the 20 TSH signals and 3 of the 6 FT4 signals are connected in a single protein network, underlying the biological interrelationship between genes regulating these traits ( Figure S3).
While our manuscript was in preparation, a GWAS of comparable sample size was published on levels of TSH in the general Icelandic population, which confirmed 15 of our reported loci (E. Porcu et al., 2011, ESHG, abstract), and inferred a role for three TSH-lowering variants in thyroid cancer [32]. Four additional TSH loci identified by Gudmundsson and colleagues were also associated in our sample-set of euthyroid individuals with p,0.05 and consistent direction of effects (VAV3, NKX2-3, TPO and FOXA2). Finally, 2 loci (SIVA1, ELK3) could not be tested because the corresponding SNPs or any surrogate (r2.0.5) were not available in our data set (Table S7). Our study shows that most of the loci described in Icelanders are reproducible in other populations of European origin; differences in sample size, phenotype definition (i.e., selection of euthyroid subjects vs general population) and in the genetic map used to detect associations most likely explain non-overlapping genome-wide significant signals. Among them, the reported signals at SOX9, ABO, SASH1, GLIS3 and MIR1179 will need to be confirmed in other studies; but one of them -GLIS3is a prime candidate, because it is involved in congenital hypothyroidism [52]. Interestingly, despite the use of variants detected through whole-genome sequencing in Icelanders, the top signals at seven overlapping loci (PDE8B, PDE10A, CAPZB, MAF/LOC440389, VEGFA, NR3C2, IGFBP5) were either coincident or in high LD (r2.0.9) with those detected in our HapMap-based meta-analysis. Thus, such variants are likely to be the causative ones.
In conclusion, our study reports the first GWAS meta-analysis ever carried out on FT4 levels, adds to the existing knowledge novel TSH-and FT4-associated loci and reveals genetic factors that differentially affect thyroid function in males and females. [32], FGF7 [56]). Furthermore, the TSH-associated variants were found to contribute to TSH levels outside the reference range. Overall, our findings add to the developing landscape of the regulation of hypothalamic-pituitary-thyroid axis function and the consequences of genetic variation for hypo-or hyperthyroidism.

Ethics statement
All human research was approved by the relevant institutional review boards, and conducted according to the Declaration of Helsinki.

Cohort details
Cohort description, genotyping and statistical methods for individual study cohorts are reported in Text S1 and Table S1.

Statistical analyses
We carried out a meta-analysis including up to 26,523, individuals from 18 cohorts for TSH and up to 17,520 individuals from 15 cohorts for FT4 (see Table 1). FT4 measures were not available for all 21,955 individuals with TSH levels of the 15 participating cohorts. We combined evidence of associations from single GWAS using an inverse variance meta-analysis, where weights are proportional to the squared standard error of the beta estimates, as implemented in METAL [57]. Prior to GWAS, each study excluded individuals with known thyroid pathologies, taking thyroid medication, who underwent thyroid surgery, and with outof-range TSH values (,0.4 mIU/L and .4 mIU/L), and an inverse normal transformation was applied to each trait (Table  S1). Age, age-squared, and gender were fitted as covariates, as well as principal components axes or additional variables, as required (Table S1). Family-based correction was applied if necessary (see Table S1). Uniform quality control filters were applied before meta-analysis, including MAF ,0.01, call rate ,0.9, HWE P,1610 26 for genotyped SNPs and low imputation quality (defined as r 2 ,0.3 or info ,0.4 if MACH [58] or IMPUTE [59,60] were used, respectively) for imputed SNPs.
Genomic control was applied to individual studies if lambda was .1.0. The overall meta-analysis showed no significant evidence for inflated statistics (lambda for TSH, FT4 and were 1.05 and 1.03 respectively). To evaluate for heterogeneity in effect sizes across populations, we used a chi-square test for heterogeneity, The table shows the association results in males and females separately for all independent SNPs associated with TSH and FT4 in the main meta-analysis (Table 2), as well as for the marker found to be associated only in females in the gender-specific meta-analysis. The last two columns report the p-value (Het P) and the false discovery rates (FDRs) for differences of effect sizes. SNPs with significantly different effect sizes at 5% FDRs and/or Bonferroni threshold (p = implemented in METAL [57]. The same test was used to evalute heterogeneity related to iodine intake, by comparing effect sizes obtained in a meta-analysis of studies assessing individuals from South Europe (InChianti, MICROS, Val Borbera, SardiNIA, totaling up to 7,488 subjects) with those estimated in a metaanalysis of studies assessing individuals from North America (BLSA, CHS, FHS, OOA, totaling up to 5,407 subjects). Finally, the main meta-analysis was carried out independently by two analysts who obtained identical results.

Conditional analysis
To identify independent signals, each study performed GWA analyses for both TSH and FT4 by adding the lead SNPs found in the primary analysis (19 for TSH, and 4 for FT4, see Table 2) as additional covariates to the basic model, and removing those from the test data set. When lead SNPs were not available, the best proxies (r 2 .0.8) were included. We then performed a meta-analysis on the conditional GWAS results, using the same method and filters as described above. We used the standard genome-wide significance cutoff (P,5610 28 ) to declare a significant secondary association.

Gender-specific analysis
To identify sex-specific effects, each study performed GWA analyses for each gender separately, using the same covariates and transformation as in the basic model (with the exception of gender covariate). We then performed a meta-analysis on association results using the same method and filters described for the primary analysis. To evaluate sex-specific differences we tested heterogeneity between effect sizes as described above. False-discovery rates (FDRs) on the 26 associated SNPs were calculated with R's p.adjust() procedure via the method of Benjamini and Hochberg [24].

Variance explained
The variance explained by the strongest associated SNPs was calculated, for each trait and in each cohort, as the difference of R 2 adjusted observed in the full and the basic models, where the full model contains all the independent SNPs in addition to the covariates. The estimates from each cohort were combined using a weighted average, with weights proportional to the cohort sample size.

Extreme phenotype analysis
To evaluate the impact of the detected variants with clinically relevant TSH levels, we compared the allele frequencies observed in different categories of individuals in a case-control approach. Specifically, we compared individuals in the upper and lower TSH tails (individuals with TSH .4 mIU/L and TSH ,0.4 mIU/L, respectively, whom were excluded for the GWAS analyses), as well as individuals in each tail with those in the normal TSH range. In the first case, individuals in the lower tail were considered controls and those in the upper tail cases. In the other two cases, we defined individuals in the normal range as controls and individuals on the two tails cases. To avoid sources of bias, individuals taking thyroid medication and/or with thyroid surgery were excluded. Only unrelated individuals were selected from the family-based cohort SardiNIA, while GEE correction was applied to the TwinsUK dataset. Results from single cohorts were then meta-analyzed. We first assessed the global impact of the 20 TSH-and 6 FT4associated variants by defining a genotype-risk score (GRS) for each individual as the weighted sum of TSH-and FT4-elevating alleles, with weights proportional to the effect estimated in the meta-analysis. For each comparison, we then calculated quartiles from the global distribution (cases+controls) of the genotype score and used quartile 1 as the baseline reference to compare the   Encodes a dual specificity phosphodiesterase abundant in the thyroid, which can hydrolyze both cAMP and cGMP to the corresponding nucleoside 59 monophosphate, but has higher affinity for cAMP, and is more efficient with cAMP as substrate. This gene was previously suggestively associated with TSH levels and hypothyroidism and linkage has been observed over this gene in families with individuals reaching the clinical criteria for sub-clinical and clinical thyroid disorders [13,34]. The top SNP is only weakly correlated with previously reported variants (r2 = 0.449 with rs2983521 and r2 = 0.184 with rs9347083), but is a perfect proxy of a SNP recently reported in association with TSH [32].
rs9472138; rs11755845 6p12 VEGFA intergenic TSH/FT4 Encodes a growth factor implicated in angiogenesis, which acts as an important regulator of both benign and malignant processes in the thyroid [62]. Angiogenesis is particularly critical for thyroid function, as the local microvasculature exerts an essential role in the continuous supply of iodine, the key element of thyroid hormone synthesis. In response to a reduction in intracellular iodine concentration, thyrocytes rapidly release angiogenic signals, including an increase in VEGFA expression and secretion [63,64].
Notably, thyroid hormone stimulation in rat brain has been shown to induce VEGFA upregulation [65], which is consistent with the nominal association of the VEGFA locus with FT4 levels observed in this study. The two independent signals, rs9472138 and rs11755845, associated with TSH levels map 40 kb downstream of VEGFA, which is the nearest gene in the region. SNPs in this region were recently reported in association with TSH levels (r2 = 0.874 and r2 = 0.947) [32].
rs13015993 2q33-36 IGFBP5 intergenic TSH It belongs to a protein family that interacts with insulin-like growth factors (IGFs) and plays a major role in regulating cell proliferation, differentiation, apoptosis and transformation. IGFBP5 is significantly over expressed in thyroid papillary carcinoma [54,66]. Studies in vitro showed that TSH, through cAMP, inhibits IGFBP5 transcription, whereas enhanced production of IGFBP5 is correlated with inhibition of thyroid function [67]; [68]. In addition, IFGBP5 has been found up-regulated in response to thyroid hormone in bone, where it interacts with the growth hormone/insulin-like growth factor (GH/IGF) system to contribute to bone formation [69], suggesting that thyroid hormone may potentate the effect of IGF-1 at the receptor level. The top SNP maps about 60 kb upstream of IGFBP5. A proxy of this SNP has been recently found associated with TSH and FT4 levels (rs737308, r2 = 0.927) [32]. rs9915657 17q23 SOX9 39UTR TSH Encodes a transcription factor involved in chondrocyte differentiation and male sex determination, although other specific functions are known. The TA domain of SOX9, which is expressed both in the pituitary and in the thyroid, has been reported to interact with a component of the thyroid hormone receptor complex (TRAP230) [70]. How this interaction could affect TSH levels is at present unclear. The top SNP maps 5 kb downstream of SOX9, which is the nearest gene in the region. rs334699 1p31.3-p31.2 NFIA intron 3 TSH Encodes a member of the NF1 (nuclear factor 1) family of transcription factors. NFI proteins have been implicated in regulating developmental processes by their specific expression pattern during embryonic development and by analysis of NFI-deficient mice [71]. In addition, they play crucial roles in the transcription of many cellular genes. Members of this family, including NFIA, have been shown to interact with thyroid transcription factor 1 (TTF1) [72], a transcription factor essential for thyroid-specific gene expression [73]. The top SNP maps in intron 3 of the gene. A proxy of this SNP has been recently found associated with TSH levels (rs334725, r2 = 1) [32].
Encodes a member of the fibroblast growth factor (FGF) family. FGF family members are involved in a variety of biological processes, including embryonic development, cell growth, morphogenesis, tissue repair, tumor growth and invasion. FGF signals play a role in the development of the thyroid gland and mice deficient for corresponding receptors show thyroid agenesis [74]. The top SNP maps in intron 2 of FGF7. A proxy of this SNP has been recently associated by GWAS with thyroid volume and goiter (rs4338740, r2 = 0.874) [36]. rs17723470 11p11 PRDM11 intron 2 TSH Encodes a member of the family of PR-domain genes involved in human cancers [75]. The function of PRDM11 and its correlation with TSH levels is unclear; however the association with low TSH values in the extreme phenotype analysis supports a role of this gene in this trait. The top SNP maps in intron 2 of the gene and is moderately correlated with a SNP associated with TSH levels (rs7128207, r2 = 0.55) [32].   [80]. SNPs in this locus have been recently reported as associated with TSH levels and thyroid cancer risk [32]. The top associated SNP maps at about 190 kb downstream of MBIP, which is the nearest gene in the region. Recently, a long, intergenic, noncoding RNA gene (lincRNA) named Papillary Thyroid Carcinoma Susceptibility Candidate 3 (PTCSC3) has been mapped 3.2 kb downstream of rs944289 (r2 = 0.708 with our top SNP), whose expression is strictly thyroid specific and acts as a PTC tumor suppressor gene [81]. rs9497965 6q24.3 SASH1 intergenic TSH Encodes a member of the SLY-family of signal adapter proteins and is a candidate tumor suppressor in breast and colon cancer. However, the biological function of SASH1 and its involvement in malignant transformation remain largely unknown. Of note, SASH1 has been identified as a downstream target of the insulin/IGF1/PI 3-kinase signaling pathway [82], which appears implicated in TSH levels in the current study. The associated SNP is located about 130 kb upstream of SASH1, which is the nearest gene in the region. rs1571583 9p24.2 GLIS3 intron 2 TSH Encodes a nuclear protein with five C2H2-type zinc finger domains, which is a member of the GLI-similar zinc finger protein family. GLIS3 functions as both a repressor and activator of transcription and is specifically involved in the development of pancreatic beta cells, the thyroid, eye, liver and kidney. Mutations in this gene have been associated with neonatal diabetes and congenital hypothyroidism (NDH) [52]. The top SNP at this locus maps in intron 2 of the gene. rs7860634 9q34.3 LHX3 intron 6 FT4/TSH Encodes a transcription factor with an essential role in pituitary development [30]. Mutations in LHX3 cause the combined pituitary hormone deficiency-3 (CPHD3) syndrome (OMIM #221750), characterized by a complete deficit in growth hormone, prolactin, gonadotropin, and TSH, a rigid cervical spine leading to limited head rotation, as well as an extended spectrum with variable sensorineural hearing loss and ACTH deficiency [31,42,43], which is consistent with its association also with TSH levels observed in this study. This locus has been also recently associated with height in Japanese [44] The top SNP maps in intron 6 of the gene.
number of cases and controls in the other quartiles. In addition, for TSH-associated variants we conducted single SNP comparisons. GRS quartile and single SNP analyses were performed by each study separately. Cohort specific results were then meta-analyzed for both the GRS score and single SNP results only if they had at least 50 cases and 50 controls. Specifically, cohorts included were: CHS, Lifelines, PROSPER, RS, SardiNIA and TwinsUK.

Bivariate analysis
Bivariate analysis was carried out with the software poly [26] in the SardiNIA cohort using the same individuals included in the GWAS and considering the same covariates and transformation for TSH and FT4 levels.

Web resources
The URLs for data presented herein are as follows:  Figure S1 Manhattan plots from meta-analysis results of serum TSH (panel A) and FT4 (panel B) levels. SNPs are plotted on the x axis according to their position (build 36) on each chromosome against association with TSH (A) and FT4(B) on the y axis (shown as -log10 P value) in. The loci highlighted in green are those that reached genome-wide significance (P,5610 28 ). In each panel, quantile-quantile plots obtained with all SNPs (red dots) and after removal of SNPs within associated regions (blue dots) are also shown. The gray area corresponds to the 90% confidence region from a null distribution of P values (generated from 100 simulations). (TIF) Figure S2 Panel A. Manhattan plots from meta-analysis results of serum TSH and FT4 levels for men and women separately are shown as indicated. SNPs are plotted on the x axis according to their position (build 36) on each chromosome; association with TSH and FT4 is indicated on the y axis (as -log10 P value). Signals reaching  [40,41]. In addition, FOXE1 is a susceptibility locus for thyroid cancer [19,[83][84][85]. The top SNP maps 25 kb upstream of FOXE1 and is also associated in our study with FT4 levels. SNPs in the FOXE1 locus were previously associated with FT4 levels in a recent candidate gene analysis (rs1443434, r2 = 0.776) [18], as well as by GWAS with both low serum TSH and T4 levels (rs965513, r2 = 1) [19], and hypothyroidism (rs7850258, r2 = 0.625) [20]. rs716822 4q33 AADAT intron 4 FT4 AADAT catalyzes the synthesis from kynurenine (KYN) of kynurenic acid (KYNA), which is implicated in the pathophysiology of several diseases of the central nervous system involving inflammation-induced brain injury [46][47][48][49].

Supporting Information
The KYN pathway has been associated with the induction of proinflammatory cytokines in the brain, which are known to activate the hypothalamo-pituitaryadrenal (HPA)-axis, involved in stress response and affecting the HPT-axis and thyroid function, including FT4 levels [45]. The top SNP maps in intron 4 of the gene.
The top SNP maps in a gene desert region, with NETO1 located about 550 kb upstream and FBXO15 about 500 kb downstream. None of these genes has a clear role in thyroid function. NETO1 is expressed in brain and encodes a predicted transmembrane protein containing two extracellular CUB domains followed by a low-density lipoprotein class A (LDLa) domain. A similar gene in mice plays a critical role in spatial learning and memory by regulating the function of synaptic N-methyl-D-aspartic acid receptor complexes in the hippocampus. FBOX15 encodes a member of the F-box protein family characterized by an approximately 40-amino acid F-box motif. SCF complexes, formed by SKP1 (MIM#601434), cullin (see CUL1; MIM#603134), and F-box proteins, act as protein-ubiquitin ligases. F-box proteins interact with SKP1 through the F box, and they interact with ubiquitination targets through other protein interaction domains.
LPCAT2 encodes a member of the lysophospholipid acyltransferase family. The enzyme may function to catalyze both the biosynthesis of platelet-activating factor and of glycerophospholipid precursors from arachidonyl-CoA and lysophosphatidylcholine. The encoded protein may function in membrane biogenesis and production of platelet-activating factor in inflammatory cells.
The associated SNP maps in intron 11 of the gene, and is also near CAPNS2, which is contained in the LPCAT2 gene. None of these genes has a clear role in thyroid function. Of note, calpain subunit genes (CAPNS1 and CAPNS2) play important functions in mammalian reproduction [25].
The genome-wide statistical significance in the gender specific analysis are shown in green. Panel B. Quantile-quantile plots are shown for all SNPs (red dots) and after removal of SNPs within associated regions (blue dots). (TIF) Figure S3 Ingenuity pathway analysis (IPA) results for candidate genes in the TSH and FT4 associated loci. A single protein network connects most of the identified loci. (PDF)  For each meta-analysis, we report the frequency of the effect allele (FreqA1), the effect size and the corresponding standard error (Effect, StdErr), the combined association pvalue (P), the pvalue for heterogeneity between studies (Het P), and the number of samples analyzed (N). The last two columns report the pvalue for differences of the effect size estimated in the two groups, and the total number of samples analyzed.
(DOC)     [32] with TSH and FT4 levels available in our data-set. When the same marker was not available, we reported a proxy (r 2 .0.8) and the relative r 2 . StdErr, standard error. Loci reaching genome-wide significance in our data set with either the same SNP or a proxy are highlighted in bold. (DOC) Text S1 Supplemental acknowledgments and funding information, cohort description, phenotyping, genotyping, analysis methods, supplemental references. (DOC)