Genetic determinants of thyroid function in children

Abstract Objective Genome-wide association studies in adults have identified 42 loci associated with thyroid stimulating hormone (TSH) and 21 loci associated with free thyroxine (FT4) concentrations. While biologically plausible, age-dependent effects have not been assessed. We aimed to study the association of previously identified genetic determinants of TSH and FT4 with TSH and FT4 concentrations in newborns and (pre)school children. Methods We selected participants from three population-based prospective cohorts with data on genetic variants and thyroid function: Generation R (N = 2169 children, mean age 6 years; N = 2388 neonates, the Netherlands), the Avon Longitudinal Study of Parents and Children (ALSPAC; N = 3382, age 7.5 years, United Kingdom), and the Brisbane Longitudinal Twin Study (BLTS; N = 1680, age 12.1 years, Australia). The association of single nucleotide polymorphisms (SNPs) with TSH and FT4 concentrations was studied with multivariable linear regression models. Weighted polygenic risk scores (PRSs) were defined to combine SNP effects. Results In childhood, 30/60 SNPs were associated with TSH and 11/31 SNPs with FT4 after multiple testing correction. The effect sizes for AADAT, GLIS3, TM4SF4, and VEGFA were notably larger than in adults. The TSH PRS explained 5.3%-8.4% of the variability in TSH concentrations; the FT4 PRS explained 1.5%-4.2% of the variability in FT4 concentrations. Five TSH SNPs and no FT4 SNPs were associated with thyroid function in neonates. Conclusions The effects of many known thyroid function SNPs are already apparent in childhood and some might be notably larger in children as compared to adults. These findings provide new knowledge about genetic regulation of thyroid function in early life.

childhood has been associated with deceleration in growth and skeletal maturation and delayed puberty, whereas hyperthyroidism in childhood has been associated with growth acceleration, advanced bone age, and delayed puberty. 2,3 Several non-genetic determinants of childhood thyroid function have been identified, such as child age, sex, ethnicity, and anthropometry. 4 The genetic heritability of thyroid function parameters has been estimated at 65% for thyroid-stimulating hormone (TSH) and 39%-80% for free thyroxine (FT4). 5,6 In a large genome wide association study (GWAS) in 2018, multiple novel genetic variants of TSH and FT4 concentrations were identified that explain 9.4% of variability in serum TSH and 4.8% of serum FT4 concentrations in adults. 7 However, it is currently unknown whether genetic determinants of TSH and FT4 concentrations in adults have a similar, or any, effect during earlier stages of life.
Complex and dynamic changes occur in thyroid function across the lifespan. 8 Most likely, childhood thyroid function is less influenced by autoimmunity, medication usage, comorbidities or aging-related factors that occur throughout life, as compared to adult thyroid function. Therefore, the genetic component of variation in thyroid hormone concentrations may be greater in children as compared to adults. While it seems apparent that the effects of individual genetic variants might be of greater importance if already apparent during childhood, different effects of thyroid hormone on growth and development may warrant age-specific effects of genetic thyroid system determinants. However, there remains a relevant knowledge gap as studies on genetic determinants of childhood thyroid function to date were small, did not report the effects of single genetic variants and have not yet been able to replicate recently discovered loci. 7,[9][10][11] Therefore, we investigated the association of adult thyroid function-related genetic variants with TSH and FT4 concentrations during childhood in three different population-based cohorts. We hypothesized that only a small number of genetic loci identified in adults would also be associated with childhood serum TSH and FT4 concentrations in the same direction and with the same magnitude of effect.

Methods
This study was embedded in three prospective birth cohorts: Generation R (The Netherlands), the Avon Longitudinal Study of Parents and Children (ALSPAC, United Kingdom), and the Brisbane Longitudinal Twin Study (BLTS, Australia).

Study design and participants
In Generation R, all pregnant women with an expected delivery date between April 2002 and January 2006 and living in Rotterdam were invited to participate. 12 Single nucleotide polymorphism (SNP) data were available for 5732 children in the first wave of data collection and for 1794 children in the second wave. Of the children of European ancestry (n = 4315), 2466 children had TSH or FT4 measurements available at the age of 6 years and 2726 children had TSH or FT4 measurements at birth. The Medical Ethics Committee of the Erasmus Medical Centre, Rotterdam, approved the Generation R study.
In ALSPAC, eligible women were those living in the former Avon area in Southwest England, United Kingdom, with an expected delivery date between April 1991 and December 1992. SNP data were available for 7975 unrelated European children. Of these, 3596 children had TSH or FT4 measurements at the age of 7 years. Ethical approval for the study was obtained from the ALSPAC Ethics and Law Committee and the Local Research Ethics Committees.
BLTS was conducted at the Queensland Institute of Medical Research (QIMR) Berghofer, Australia and began in 1992 with approximately 100 pairs of 12-year-old twins (alongside their non-twin siblings and parents) with additional participants recruited annually thereafter. 13 In BLTS, SNP data were available for 2832 children. Of these, 1746 children had TSH or FT4 measurements available at the age of 12 years. Ethical approval for the study was obtained from the Human Research Ethics Committee of the QIMR Berghofer Medical Research Institute.
Exclusion criteria were pre-existing thyroid disease, thyroid medication, and TSH concentrations outside the cohortspecific 2.5-97.5th centile range. In addition, one child was randomly excluded from each participating sibling pair in Generation R in cord blood and child analyses. Written informed consent was obtained from all participants and/or the children's parents or guardians in all three cohorts. The study was performed in accordance with the principles of the Declaration of Helsinki.

Genetic variants and risk scores
We extracted data on 60 TSH and 31 FT4 SNPs that were previously identified in adults 7 from readily available SNP arrays in each cohort. These SNPs were identified in the most recent GWAS that studied both TSH and FT4 SNPs and were independently associated with TSH or FT4 in adults. The 28 TSH SNPs that were identified hereafter 14 were accessible in Generation R and BLTS only and analysed in a post-hoc analysis.
In Generation R, SNP data were generated in two separate rounds of genotyping and imputation with possible batch effects and are therefore analysed in two separate subsets in this study. ALSPAC and BLTS were also analysed as two separate cohorts.
Unweighted polygenic risk scores (PRSs) were calculated by summation of allele dosages for each SNP. Weighted PRSs were calculated by multiplying the dosage data of each SNP by the previously reported effect estimate in adults. 7 The weighted and unweighted PRSs were rescaled to range from 0 to 10.

Thyroid function measurements
In all three cohorts, all participants were invited as part of the follow-up to provide blood samples, and no exclusion criteria were applied beforehand.
In Generation R, serum samples in children were obtained at the age of 6 years, and cord blood samples were obtained directly after birth. Plain tubes were centrifuged, and serum was stored at −80 °C. TSH and FT4 were measured using an electrochemiluminescence immunoassay on the Cobas e601 immunoanalyzer (Roche Diagnostics). The intra-and interassay coefficients of variation were 1.1%-3.0% for TSH at a range of 0.04-0.4 mU L −1 and 1.6%-5.0% for FT4 at a range of 1.6-24.1 pmol L −1 .
In BLTS, serum TSH, FT4, and FT3 were measured at age 12 years using a chemiluminescent immunoassay on an Abbott ARCHITECT analyser (Abbott Diagnostic). 15

Statistical analyses
To make concentrations comparable across cohorts and to approach normality of the residuals, TSH, FT4, and FT3 values were analysed after inverse normal transformation, resulting in SD-scores for the thyroid hormone concentrations. We used multivariable linear regression models to study the association of the single TSH and FT4 SNPs as well as the two respective PRSs with child and neonatal TSH and FT4. Analyses were carried out separately for children and neonates. The explained variance was calculated as the difference in explained variance between the basic model (with covariates) and the full model (including covariates and additionally the PRS).
Assumptions of linear regression models, including linearity, homoscedasticity, and normal distributions of the model residuals, were met for all models in the three cohorts. All models were adjusted for age at blood sampling, child sex, and the first four principal components of genetic ancestry in Caucasians to account for population stratification. In BLTS, we used a linear mixed effect model with a nested random effects structure with separate intercepts for each subject within zygosity group within family to account for relatedness. A symmetric correlation structure between observations at the lowest level of nesting was assumed. Results from the three cohorts were pooled and analysed in a fixed-effects metaanalysis using a two-step approach. To take heterogeneity into account, we subsequently ran models with random effects for SNPs with moderate heterogeneity (I 2 > 30). 16 However, we expected the effects sizes to be similar across cohorts. Furthermore, the statistical power of a random-effects model is always lower than a fixed-effects model and further decreases with higher heterogeneity. 17 Therefore, we first used models with fixed effects only and subsequently added random effects in a sensitivity analysis. The genetic data were generated in two runs in Generation R and were therefore analysed as two different subsets (wave 1 and wave 2).
Several additional analyses were performed. First, the analyses of FT4 concentrations were repeated in children with FT4 concentrations in the cohort-specific 2.5-97th centile range. Second, we studied the association of the FT4 SNPs and PRS with FT3 concentrations in children from ALSPAC and BLTS. Third, the TSH and FT4 PRSs were studied as determinants of FT4 and TSH concentrations, respectively. Fourth, we constructed different PRSs including using SNPs associated with thyroid function in childhood, applying estimates from one cohort to another cohort, and excluding SNPs which showed heterogeneity or large deviation from adult estimates. Fifth, we performed a univariate metaregression with mean age per cohort as predictor and the effect size or the explained variance of FT4 as the dependent variable, using the sample sizes of the cohorts as weights.
To adjust for multiple testing, a false discovery rate (FDR) correction for the number of individual SNPs per outcome was applied (ie, 60 for TSH and 31 for FT4) using the Benjamini and Hochberg method with a q-value threshold of 0.05. 18 All statistical analyses were performed with R statistical software version 4.2.1 and the R metafor package for the meta-analysis.
Further details on exclusion criteria and measurements are in the Supplemental Methods.

Results
The final study population comprised 7231 children in the childhood thyroid function analyses ( Figure 1). Descriptive statistics are shown in Table S1. Median (95% range) age at serum thyroid measurements was 6.0 (5.7-7.5) years in Generation R, 7.5 (7.3-8.9) years in ALSPAC, and 12.1 (12.0-12.5) years in BLTS. Median TSH concentrations were 2.4 mU L −1 in Generation R, 2.1 mU L −1 in ALSPAC, and 1.5 mU L −1 in BLTS. Median FT4 concentrations were 16.5 pmol L −1 in Generation R, 15.6 pmol L −1 in ALSPAC, and 12.6 pmol L −1 in BLTS. The minor allele frequencies were similar across cohorts (Table S2). The final study population for newborn thyroid function analyses consisted of 2388 neonates.

Individual SNPs
Thirty out of 60 SNPs were associated with childhood TSH concentrations and 11 out of 31 SNPs with childhood FT4 concentrations after correction for multiple testing (Tables 1  and 2, Figures 2 and 3). The direction of effect was similar to adults for 55 out of 60 TSH SNPs and 28 out of 31 FT4 SNPs. None of the SNPs showed an inverse direction of effect in all cohorts as compared to adults (Tables S3 and S4).
For 29 out of 60 TSH SNPs and 9 out of 31 FT4 SNPs, the effects estimates were within ± 30% deviation (arbitrary cutoff to provide an indication) from the adult effect estimates (Tables 1 and 2). Twenty-three (77%) SNPs associated with childhood TSH and 8 (73%) associated with childhood FT4 exceeded the adult effect estimate. The FT4-associated SNP in the AADAT locus and three TSH-associated SNPs in the GLIS3, TM4SF4, and VEGFA locus had effect estimates that were notably larger than those in adults.
When SNPs with an I 2 value >30 were analysed using a random effects model, 6 out of 30 TSH SNPs, and 1 out of 11 FT4 SNPs that were statistically significant in the fixed effects model were not associated with TSH or FT4 anymore (Tables 1  and 2).
In newborns, 5 TSH SNPs were associated with TSH concentrations, and these SNPs were also associated with TSH in childhood (Table S5 and Table 1). None of the FT4 SNPs identified in adults were associated with FT4 concentrations in newborns (Table S6).
The weighted PRS including TSH SNPs was associated with TSH concentrations in newborns, whereas the FT4 PRS was not associated with FT4 concentrations in newborns (Table S7).
The effect estimates for unweighted PRSs were comparable to the weighted PRSs, albeit with lower explained variabilities for the variation in TSH and FT4 (Tables S8 and S9).

Additional analyses
Four FT4 SNPs were associated with FT3 concentrations after multiple testing correction, but three of those showed substantial heterogeneity with I-squared values >80 and none were associated with FT3 in models with random intercepts (Table 5). However, 8 SNPs were associated with FT3 in ALSPAC only (Table S10). The FT4 PRS was associated with FT3 concentrations (β = −.04, SE = 0.01, P≤.001) and the explained FT3 variability was 0.2% in ALSPAC and 0.3% in BLTS (Table 4).
Results were similar after excluding FT4 concentrations outside the normal range (Tables S11 and S12).
Small differences in explained variance were apparent when the PRS was constructed using different approaches (Table S13). The explained variance increased in all cohorts when only including SNPs that were associated with TSH or FT4 in childhood, ranging 6.5%-7.4% instead of 5.3%-6.7% for TSH and 2.0%-5.4% instead of 1.5%-4.2% for FT4. When only SNPs within ± 30% deviation from adults or SNPs without heterogeneity were included, the explained variance decreased. The clearest decline was found in ALSPAC and Generation R for FT4, in which the explained variances changed from 2.9% to 0.8% and from 4.2% to 1.5%, respectively.
Three out of 28 SNPs that were discovered after the previously mentioned GWAS were associated with TSH concentrations in childhood (Tables S14 and S15).

Discussion
In this meta-analysis of three population-based prospective cohorts, 30 out of 60 TSH SNPs and 11 out of 31 FT4 SNPs identified in adults were also associated with childhood TSH and FT4. The weighted polygenic risk score (PRS) including TSH SNPs explained 5.3%-8.4% of childhood TSH variability, whereas the FT4 PRS explained 1.5%-4.2% of childhood FT4 variability. Furthermore, we identified five SNPs associated with TSH at birth. Our results seem to suggest that the effect sizes of some SNPs are notably larger in children as compared to adults.
In the current study, multiple SNPs that were associated with thyroid function in adults were also associated with TSH or FT4 during childhood. The explained variability in FT4 concentrations as assessed by the PRS differed between the cohorts, with what seemed to follow a trend of a lower explained variability with a higher median age of the cohort and this finding is supported by a negative association of age with the explained variance. In line with this, thyroid function is likely to be less explained by genetics and more variable within a population as participants get older, accumulate exposure to environmental factors and acquire thyroid autoimmunity.     In our study, however, the explained variability of TSH and FT4 as assessed by the PRS including the same SNPs as in adults was lower than in adults (ie, 9.4% for TSH and 4.8% for FT4 in adults). A possible explanation is that the SNPs that we studied were identified in adults of which some might be important for auto-immunity, which is likely less involved in childhood thyroid function. Moreover, there might be age-specific effects of genetic thyroid system determinants that are only apparent in childhood and therefore were not identified in adults and subsequently not investigated in our study. Therefore, future studies should include a GWAS on childhood thyroid function to identify childhood-specific genetic determinants of thyroid function. In addition, future studies should investigate gene-age interactions and the role of epigenetic regulation. This could help to increase the explained variability of TSH and FT4. Importantly, there could be other factors than age contributing to the difference in explained variability of FT4 concentrations between the cohorts, such as iodine status. Iodine intake in the Netherlands (Generation R) is adequate, whereas the United Kingdom (ALSPAC) is moderately iodine-deficient and iodine deficiency re-emerged in Australia in the 1990s when BLTS participants were recruited. [19][20][21] The contribution of iodine status to the variation in FT4 concentrations as compared to the genetic contribution might be higher in countries with a median insufficient iodine status and a wider range of iodine levels, as compared to countries with an adequate iodine intake. Indeed, the highest explained variability in FT4 concentrations was observed in Generation R, whereas the explained variability was lower in ALSPAC and BLTS.
Since childhood estimates from previous studies were not available, we were limited by using adult estimates for generating PRSs. Calculating PRSs with different methods did not yield large differences in most analyses. The decline in explained variance of FT4 when only SNPs without heterogeneity were included in the PRS in Generation R and ALSPAC (from 4.2% and 2.9% to 1.5% and 0.8%, respectively) can be explained by the fact that three out of five excluded SNPs had large effect estimates in adults, Generation R and/or ALSPAC and thus contributed to a relative large extent to the explained variance. The effect estimates of these SNPs in the DIO1, AADAT, and SLC25A52 locus were smaller in BLTS, resulting in heterogeneity in the meta-analysis.
We show that a higher TSH PRS was associated with lower FT4 concentrations and a higher FT4 PRS was associated with higher TSH concentrations in our meta-analysis in childhood. This is not in line with the classic TSH stimulation model including negative feedback of T4 on the level of the hypothalamus and pituitary. Therefore, our results suggest that the identified TSH SNPs mostly act through changes in the TSH setpoint or TSH (receptor) sensitivity, while the identified FT4 SNPs are likely to act through intracellular mechanisms including sulphation, deiodination, or glucuronidation. This is in line with the limited overlap between TSH and FT4 genetic determinants observed in adults, as well as the lack of associations between FT4 SNPs and hypo-or hyperthyroidism in adults. 7 These findings underline the importance of understanding the underlying mechanisms of individual SNPs in order to interpret their associations with clinical outcomes.
While genetic studies are less prone to confounding bias as compared to observational studies in general, inferring causality can still be hampered by possible pleiotropic effects, ie, genetic variants can influence determinants of thyroid function such as BMI and subsequently affect thyroid function. Indeed, some genetic variants of thyroid function have been associated with BMI, height, and weight circumference, [22][23][24][25][26] but none of these genetic variants were associated with childhood thyroid function in the current study. In addition, the pleiotropic genetic variants might influence BMI through thyroid function instead of vice versa.
In this study, the recently discovered genetic variant of the enzyme AADAT was associated with FT4 concentrations with an effect estimate almost twice the size of that in adults. 7 This could be caused by differential expression of AADAT with age, as several genes show age-dependent expression in target tissues of thyroid hormone including the brain, skin and adipose tissue. 27 Therefore, the association between genetic variation of AADAT with FT4 concentrations might be more apparent at child age through higher expression of AADAT. To our knowledge, however, it is not known if AADAT expression also differs with age. In-vitro studies showed AADAT-dependent conversion of T4 and T3 to their metabolites. The AADAT variant was also associated with FT3 in our study, but was not significant anymore when including random intercepts to account for heterogeneity, which is likely explained by the different directions of effect. The AADAT variant was positively associated with FT3 in our meta-analysis and in ALSPAC, yet amongst the older children in BLTS, as well as in the adult GWAS there was a negative association of the AADAT variant with FT3. 7 These results of an age-dependent magnitude and direction of effect might suggest that the effect of AADAT on thyroid function is age specific. However, this potential age-dependent direction of effect is based on limited data and should be considered as a hypothesis that requires further replication.
Our results seem to suggest that the effect sizes of three TSH loci, GLIS3, TM4SF4, and VEGFA are also notably larger in children than in adults. This difference could be related to its underlying biology related to early life rather than aging as GLIS3 is associated with thyroid development and mutations in this gene are associated with neonatal diabetes and congenital hypothyroidism. In addition, GLIS3 is also involved in the development of the eye, liver, kidney, and pancreatic beta cells. 28 Similarly, TM4SF4 is a member of a protein family involved in the regulation of cell development, activation, growth and motility. 28 Thus, the large effect of GLIS3 and TM4SF4 on childhood TSH concentrations might be explained by their involvement in organ and cell development, which are likely more evident during early life. VEGFA promotes conversion of T4 to T3 which consequently suppresses TSH concentrations. 29 The large effect size of this locus can be explained by its age-dependent expression, as animal studies have shown that VEGF expression decreases with age. 30,31 The loci with the largest effect sizes in childhood relative to adults have been associated with an elevated blood glucose level (AADAT) and diabetes (GLIS3). This might suggest that any relation between thyroid function and diabetes in childhood could be caused by pleiotropic genetic effects, whereas a suboptimal thyroid function in adulthood could result in long-term metabolic effects. In contrast, the loci with the smallest effect sizes in childhood relative to adults have been associated with weight (FOXE1) and height (CAPZB). This suggests that any relation between anthropometry and thyroid function in adulthood might be partially explained by pleiotropic genetic effects, whereas anthropometry might have rather direct effects on thyroid function in childhood. To our knowledge, this is the largest study on genetic determinants of childhood thyroid function. Previous studies included either small subsets of children, 9 assessed the combined effects of a much smaller amount of currently known genetic determinants of thyroid function 7,10 or performed a GWAS approach in adolescents in a relatively small subset of our study sample. 11 We were able to perform this study using a multi-cohort approach with a large sample size available for analyses. This enabled us to optimize the generalizability of the results of our study by replicating adult genetic determinants and to crossreplicate our own findings. An important limitation of our study is that the design did not allow for identification of new, childhood thyroid function SNPs. We did not perform a meta-GWAS as we focused on known SNPs in adults and an even larger sample size would be preferable for this approach. Importantly, the difference in power between our study and the GWAS in adults might partially explain why we did not replicate more thyroid function SNPs. The fact that only a few of the most recently discovered TSH SNPs were replicated in children could also be a power issue, since the largest cohort was not included in this analysis and the sample size of the GWAS in which these extra TSH SNPs were identified constituted more than one and one half times the size of the previous GWAS. Future studies should include a GWAS in childhood, as new treatment targets and strategies such as the incorporation of genetics in risk stratification for treatment decisions might evolve by improving the knowledge on normal regulation of childhood thyroid function.
In conclusion, our study provides new data on genetic regulation of childhood thyroid function. Our results are evidence that the effects of many known genetic variants are already apparent in childhood. Further, we cautiously interpret our results as suggesting that the extent of the effect of certain genes may be age-specific. These findings advance the understanding of child thyroid function and aid in untangling the effects of maternal and child thyroid function on offspring outcomes.

Contributors
Tessa Mulder performed the data analyses and was involved in writing the report. Purdey Campbell performed the analyses in the BLTS cohort and was involved in writing the report. Peter Taylor, Scott Wilson, Marco Medici, Colin Dayan, Vincent Jaddoe, John Walsh, and Nicholas Martin contributed to writing of the report. Robin Peeters and Henning Tiemeier contributed to data analyses, data interpretation, and writing of the report. Tim Korevaar supervised data analyses, contributed to data interpretation and writing of the report, and directed the project. Conflicts of interest: All authors (T.M., P.C., P.T., S.W., M.M., C.D., V.J., J.W., N.M., R.P., H.T., and T.K.) have no declarations of interest.

Data availability
The datasets generated during and/or analysed during the current study are not publicly available but may be accessed through the corresponding author on reasonable request.