Effects of Thyroid Status on Regional Brain Volumes

A Diagnostic and Genetic Imaging Study in UK Biobank

Tom Chambers; Richard Anney; Peter N. Taylor; Alexander Teumer; Robin P. Peeters; Marco Medici; Xavier Caseras; D. Aled Rees


J Clin Endocrinol Metab. 2021;106(3):688-696. 

In This Article



The UK Biobank is a population-based, volunteer cohort of 500 000 participants (40–69 years of age at recruitment) who have provided extensive medical, lifestyle, and genetic data.[9] Our study used the brain MRI imaging-derived phenotype (IDP) data of 22 000 participants who had been scanned by UK Biobank (2014–2018) and had their data released by the time of study initiation. Using International Statistical Classification of Diseases and Related Health Problems revision-10 (ICD-10) primary and secondary diagnostic codes from participants' hospital inpatient records (from 1996 to present) and from self-reports of medical conditions (provided at UK Biobank centre), we removed individuals with severe neurological and psychiatric conditions (Supplementary note 1[10]). Using the thyroid-relevant hospital inpatient ICD-10 code recordings (E00-E07), we derived a categorical variable of having had anytime recorded diagnoses for hypothyroidism, hyperthyroidism, or other thyroid disorders (Supplementary note 1[10]). For genetic analyses, we further removed those subjects with ICD-10 or self-reported thyroid disorder. Additional imaging and genetic exclusion criteria and quality control are outlined below. Ethical approval for UK Biobank was granted by the North West Multi-Centre Ethics Committee. Data for this study were obtained under application number #17044.

MRI Data

A full description of the imaging acquisition, quality control and imaging-derived phenotype (IDP) generation from UK Biobank can be found elsewhere.[11] Briefly, UK Biobank participants attended 1 of 2 imaging sites using the same scanner design (3-Tesla Siemens Skyra scanner; 32-channel head coil) and underwent an extensive imaging protocol, which included a T1-weighted structural scan (3D Magnetization Prepared Rapid Acquisition Gradient Echo with 1mm[3] isotropic resolution), whose data we used here. We used the released IDP produced by UK Biobank using processing tools from FMRIB (Functional Magnetic Resonance Imaging of the Brain) Software Library[12] (FSL), including FAST[13] (FMRIB's Automated Segmentation Tool) registration of cerebellar lobules atlas and FIRST[14] (FMRIB's Integrated Registration and Segmentation Tool) registration of subcortical regions atlas. Of the cerebellar lobule IDPs, we excluded Crus I vermis due to its small size and likely low signal-to-noise ratio. We grouped these regions of interest by hemisphere, and all cerebellar lobules into a single total cerebellar volume measure, analyzing any lobule-specific effects as secondary analyses. We chose this approach because no cerebellar-specific registration process was applied, likely reducing the signal-to-noise ratio in the cerebellum, because cerebellar lobules are highly correlated and because the lobule boundaries show only small association with functional boundaries, meaning that detection of strong lobule-specific effects was deemed less likely. Subjects with any covariate imaging value (see below) deemed as outliers (>5 × median absolute deviation from the overall median) were removed from further analysis and, of these, those with remaining outlier values were excluded from each respective region analysis (Supplementary Table 1[10]).

Genetic Data

A full description of UK Biobank's data collection, quality control, and imputation process can be found elsewhere (http://www.ukbiobank.ac.uk/scientists-3/genetic-data/). Locally, we further applied additional quality control of the imaging-sample raw genotypes using the genotypeqc function (https://github.com/ricanney/stata). All Stata functions described leverage PLINK[15] (v1.90b5.4; www.cog-genomics.org/plink/1.9/). Briefly, all markers were harmonized to genome build hg19 and common nomenclature was applied based on the Haplotype Reference Consortium (HRC) r1.1. We excluded markers based on individual missingness (>2%), low minor allele count (<5), deviations from Hardy-Weinberg equilibrium (P < 10–10), and from expected minor allele frequency (MAF; defined as 4 SD from the reported 1000 Genomes Project phase-3 GBR MAF). Participants were excluded based on their overall marker missingness (>2%) and marker heterozygosity (>4 × SD from sample mean). Using the bim2ancestry_keep function (https://github.com/ricanney/stata), we limited to GBR ancestry (>4 × SD from 1000 Genomes Project phase-3 GBR sample mean for the first 3 principal components) and using bim2unrelated we removed one of each pair of close relatives (>0.0442 kinship coefficient ie, third-degree relatives). A total of 7 726 488 markers were included in downstream analysis.

Polygenic score training data were obtained from the most recent genome-wide association study (GWAS) meta-analysis on circulating thyroid hormone levels, which tested up to 72 168 participants of European ancestry and provided summary statistics for thyroid-stimulating hormone (TSH) and free thyroxine (fT4) levels in subjects with TSH levels within the cohort reference range, as well as for elevated (hypothyroidism) and reduced (hyperthyroidism) TSH compared with these ranges.[16] Each thyroid disorder GWAS was based on 7 858 695 and 7 980 324 genetic markers, respectively. We utilized the summary statistics across genders, rather than sex-specific summary statistics. Summary statistics were harmonized to hg19 build and HRC reference nomenclature using summaryqc (https://github.com/ricanney/stata).

Polygenic scores were created using the PLINK summaryqc2score wrapper function (https://github.com/ricanney/stata). Regions of known high linkage disequilibrium[17] were excluded prior to calculating polygenic scores. Scores were subsequently calculated from linkage disequilibrium–independent (r[2] < 0.2) markers for hyper- and hypothyroidism, and TSH and free-T4 level summary statistics for single nucleotide polymorphism (SNP) inclusion P value thresholds (P T) of P T < 0.5, <0.1, <0.05, <0.01, <1e−3, <1e−4, <1e−5, <1e−6, <1e−7, and <1e−8. An individual's polygenic score is a linear weighted sum of their genotype values weighted by their effect sizes for each thyroid trait.

Statistical Methods

For our analyses we used univariate multiple linear regression in R(3.6.0) (https://www.R-project.org/). Our first primary analyses investigated the effect of hypothyroid and hyperthyroid diagnoses on total cerebellar and subcortical volumes. The control group in this analysis were those with no thyroid-related diagnoses. For each model, we controlled for potential confounding by other imaging variables of head size (volumetric scaling applied for registration, analogous to the inverse of head size; 25000); imaging center attended (54-2.0); X-, Y-, and residual Z-head position in the scanner (25756, 25757, 25758, with Z-position having starting table-Z-position regressed out; 25759); along with demographic covariates of age (21003-2.0, quadratic); sex;[18] and their interaction. Histograms were used to confirm normally distributed residuals. We tested for any mediation of body mass index (BMI) on overt diagnoses differences using the mediation package[19] with 1000 bootstrapping simulations, while also controlling for the above covariates. Using the medsens function,[19] we tested for violations of the sequential ignorability assumption that there are no unmeasured confounders of the mediator-outcome pathway. In the subset with no overt or self-reported thyroid disorders, we repeated the above analysis using polygenic scores for thyroid-related traits. For these genetic predictors, in addition to the above covariates, we also controlled for BMI (log transformed; 21001–2.0) as well as the first 10 genetic principal components to correct for residual population structure. Secondary analyses investigated any lobule-specific effects. Standardized β coefficients reflecting changes in units of SD, corresponding 95% CI, and the variance explained in the outcome uniquely attributable to the phenotype of interest (ΔR2; calculated by subtracting each model R2 with the phenotype against those without) are provided. The Benjamini-Hochberg method[20] was used to control the false discovery rate (FDR < 0.05) for the number of models assessed in each section (thyroid diagnoses and polygenic scores).

We also aimed to investigate evidence for regional pleiotropy between ADHD and hypo- and hyperthyroidism using the GWAS-pairwise method[21] utilizing recent ADHD GWAS meta-analysis summary statistics of 20 183 ADHD cases and 35 191 controls.[22] Following application of the same quality control protocols outlined above, 5 907 045 genetic markers were carried forward for analyses. We performed 2-sample meta-analysis using METAL,[23] with the second pair in the GWAS performed with the reported β and with the sign of the β flipped so as to not avoid omitting those markers with opposing effects. Three models were examined; model-1 where the association is driven by thyroid trait GWAS (gwas-1), model-2 where it is driven by ADHD GWAS (gwas-2), and model-3 where it is driven by signal at both gwas-1 and gwas-2. For any marker, model-3 was accepted if P model-3 < 5e-8 and independently P gwas-1 and P gwas-2 were both observed at P < 1e-5.