Revealing the genetic complexity of hypothyroidism: integrating complementary association methods

Abstract

Hypothyroidism is a common endocrine disorder whose prevalence increases with age. The disease manifests itself when the thyroid gland fails to produce sufficient thyroid hormones. The disorder includes cases of congenital hypothyroidism (CH), but most cases exhibit hormonal feedback dysregulation and destruction of the thyroid gland by autoantibodies. In this study, we sought to identify causal genes for hypothyroidism in large populations. The study used the UK-Biobank (UKB) database, reporting on 13,687 cases of European ancestry. We used GWAS compilation from Open Targets (OT) and tuned protocols focusing on genes and coding regions, along with complementary association methods of PWAS (proteome-based) and TWAS (transcriptome-based). Comparing summary statistics from numerous GWAS revealed a limited number of variants associated with thyroid development. The proteome-wide association study method identified 77 statistically significant genes, half of which are located within the Chr6-MHC locus and are enriched with autoimmunity-related genes. While coding GWAS and PWAS highlighted the centrality of immune-related genes, OT and transcriptome-wide association study mostly identified genes involved in thyroid developmental programs. We used independent populations from Finland (FinnGen) and the Taiwan cohort to validate the PWAS results. The higher prevalence in females relative to males is substantiated as the polygenic risk score prediction of hypothyroidism relied mostly from the female group genetics. Comparing results from OT, TWAS, and PWAS revealed the complementary facets of hypothyroidism's etiology. This study underscores the significance of synthesizing gene-phenotype association methods for this common, intricate disease. We propose that the integration of established association methods enhances interpretability and clinical utility.

Keywords: FinnGen; GWAS; Hashimoto’s thyroiditis; PWAS; UK-Biobank; congenital hypothyroidism; genotyping; open targets.

FIGURE 1

Summary of independent loci identified from major GWAS results as compiled in the OTG portal. (A) The number of participants in each study and the number of hypothyroidism cases are indicated by N (all) and n (cases). There are 21 variants that are shared by all six studies (colored red). The chromosomal position is shown (bottom, light blue). (B) STRING analysis of the 21 mapped associated genes resulted in a network of 13 genes (interaction score >0.4). The nodes are colored by PPI clusters. Evidence of connectivity between the clusters is indicated by dashed lines. (C) Connectivity of the 21 associated genes (Table 1) and their major functional annotations.

FIGURE 2

GWAS results for hypothyroidism (ICD-10, E03) with ∼10.2 M variants. (A) Manhattan plot covering Chr. 1 to Chr. 22. For visualization clarity we capped the p-value at <1e-60. Red frame indicates the MHC extended locus on Chr6. (B) Zoom in of the Manhattan plot covering part of the extended region of MHC from Chr 6. The significant threshold of 5e-05 and 5e-08 are marked by the red and green horizontal lines, respectively. (C) Quantile-quantile (Q–Q) plot based on the results of GWAS using 10.2 M variants. The red line shows that there is no signal in the data, the inflation factor l = 1.089.

FIGURE 3

Partition of the significant coding GWAS variants at different thresholds. (A) Position of the variants in the Chr6 MHC locus and in other locations. We consider the MHC locus to span 6 M based in the MHC region of Chr6. (B) Partition according to the trend of variants that are protective or increase the risk for hypothyroidism.

FIGURE 4

Associated genes from PWAS results. (A) Statistically significant genes from PWAS for ICD-10 E03 with q-value <1e-07 (total 26 genes). Genes with an increased and decreased risk are colored purple/red and blue, respectively. (B) Effect size (Cohen’s d) for PWAS results for the dominant model. The genes within the dashed frames are associated with Cohen’s d >|0.06|. Positive (green font) and negative (red font) Cohen’s d values are associated with reduced and increased risk, respectively. Supplementary Material S2: Supplementary Table S5 lists all genes and their statistics.

FIGURE 5

Venn diagram for the overlapping genes according to multiple association studies protocols. (A) Venn diagram of gene-length GWAS (g-GWAS, 134 genes), coding GWAS (c-GWAS, 72 genes with p-value <5E-07) and PWAS. Each of the overlap section is shown by label the gene as part of the MHC locus (blue) or others (orange). (B) Venn diagram of PWAS (77 genes), GWAS (OT, by genetic association score >0.5 (138 genes), and transcription-based association study (TWAS, see Methods) for hypothyroidism/myxedema by UTMOST model (71 genes; Supplementary Material S2: Supplementary Table S7). The subsets of overlapping genes are color-coded according to their main functional annotations.

FIGURE 6

Network relationship and functional enrichment of PWAS results (77 genes). (A) The STRING network represents the genes connected at an interaction score >0.9. Dashed lines mark the connections between clusters. The unified function for each cluster is colored and annotated (e.g., antigen processing). (B) Enrichment analysis using the FUMA-GWAS Gene2Func protocol. In red, the fraction of genes in the gene set; blue, the adjusted p-value; orange, the overlapping genes for each term. The top 13 KEGG pathways and bottom, the GO_MF annotations. Note the enrichment of MHC genes (HLA-DPA1, HLA-DRB1, HLA-DRB5, HLA-B, HLA-DPA1, HLA-G) in KEGG and GO-MF analyses.

FIGURE 7

Genetic association with GWAS compiled by OT. (A) Ranked genes by their genetic association (by GA score, total of 715 genes). The overlap of 77 PWAS genes and 222 OT genes with GA scores >0.3 for all 715 genes. (B) A network relation of genes that are ranked by the OT global score >0.5 for the phenotype of permanent CH (total 36, 22 are connected, STRING PPI confidence score >0.7). The nodes are colored according to the match with the findings of CH causal genes from independent cohorts from India (Kollati et al., 2020) and China (Wang et al., 2020).

FIGURE 8

Gene-based association analysis by sex (A) Distribution of a polygenic risk score (PRS) among individuals with and without E03 diagnosis, marked as cases (pink) and controls (blue). PRS scores were calculated for all GWAS (A) and coding-GWAS (B) for the entire cohorts (both), females and males. (C) PRS prediction by the coefficient of determination (R ², left) AUC-ROC (right) for coding GWAS (orange) and all GWAS (blue) for the entire cohort (both) and by sex. Coding GWAS variants partitioned by sex are listed in Supplementary Material S2: Supplementary Table S9.