https://orcid.org/0000-0002-0940-4517
Abstract
Congenital hypothyroidism due to thyroid dysgenesis (TD), presented as thyroid aplasia, hypoplasia or ectopia, is one of the most prevalent rare diseases with an isolated organ malformation. The pathogenesis of TD is largely unknown, although a genetic predisposition has been suggested. We performed a genome-wide association study (GWAS) with 142 Japanese TD cases and 8380 controls and found a significant locus at 2q33.3 (top single nucleotide polymorphism, rs9789446: P = 4.4 × 10−12), which was replicated in a German patient cohort (P = 0.0056). A subgroup analysis showed that rs9789446 confers a risk for thyroid aplasia (per allele odds ratio = 3.17) and ectopia (3.12) but not for hypoplasia. Comprehensive epigenomic characterization of the 72-kb disease-associated region revealed that it was enriched for active enhancer signatures in human thyroid. Analysis of chromosome conformation capture data showed long-range chromatin interactions of this region with promoters of two genes, FZD5 and CCNYL1, mediating Wnt signaling. Moreover, rs9789446 was found to be a thyroid-specific quantitative trait locus, adding further evidence for a cis-regulatory function of this region in thyroid tissue. Specifically, because the risk rs9789446 allele is associated with increased thyroidal expression of FDZ5 and CCNYL1 and given the recent demonstration of perturbed early thyroid development following overactivation of Wnt signaling in zebrafish embryos, an enhanced Wnt signaling in risk allele carriers provides a biologically plausible TD mechanism. In conclusion, our work found the first risk locus for TD, exemplifying that in rare diseases with relatively low biological complexity, GWAS may provide mechanistic insights even with a small sample size.
Introduction
Congenital hypothyroidism occurs in 1 in 3000 newborns and is thus one of the most prevalent rare diseases (which are defined as fewer than 1:2000 in the European Union and 1 in 2500 in Japan) (1). Untreated congenital hypothyroidism can lead to irreversible growth restriction and intellectual disability. However today, early treatment with levothyroxine based on newborn screening programs enables a near normal outcome for the hitherto severely retarded children (2). Congenital hypothyroidism can be classified into two phenotypic forms: (i) thyroid dysgenesis (TD), in which the thyroid gland is absent (aplasia), small (hypoplasia) or abnormally located (ectopia), and (ii) thyroid dyshormonogenesis with several inborn errors of thyroid hormone biosynthesis (3). TD is the predominant form affecting 80–90% of congenital hypothyroidism cases (3).
Identification of the TD-associated risk locus at 2q33.3. (A) Manhattan plot of the χ2 test P values of GWAS. A broken line indicates the genome-wide threshold for significance (5 × 10−8). Genotyped SNPs and imputed SNPs are plotted according to the chromosomal locations with the –log10(P values) from the GWAS. (B) The upper panel shows a zoomed image of the region in chromosome 2 (hg19; chr2:208600001-209200000), showing a significant association in GWAS. Genotyped SNPs (circles) along with imputed SNPs (boxes) are shown. The lower panels display the linkage disequilibrium blocks in the 1000 Genomes project data (East Asian and European).
In contrast to other rare diseases that are diagnosed in newborn screening programs (e.g. cystic fibrosis and phenylketonuria), the pathogenesis of TD remains largely unknown with non-Mendelian features such as sporadic occurrence (4), discordance in monozygotic twins (5) and female predominance (4,6). TD is essentially considered as an isolated organ malformation, because 95% of cases have no extra-thyroidal anomalies (7). Although many efforts have been made to clarify the molecular mechanisms of TD including exome sequencing (8) and methylome analysis (9), only rare syndromic forms associated with coding variants in genes such as PAX8 (10), NKX2-1 (11), FOXE1 (12), GLIS3 (13), NTN1 (14) and JAG1 (15) have been observed in patient cohorts. Despite the general sporadic occurrence of TD, in a French nationwide survey, familial case of TD was observed in 2.3% of patients, which was 15-fold higher than expected by chance alone (4). In addition, it is known that prevalence of TD differs by the ethnic background (16). Therefore, it can be assumed that TD, like other more frequent human malformations like non-syndromic orofacial clefts (17,18), is multifactorial in cause, with a genetic contribution.
Here, we conducted a genome-wide association study (GWAS) of TD to discover variants that might confer the genetic risk. Our GWAS identified a single genome-wide significant signal in 2q33.3, and follow-up analyses revealed that the genetic risk may be linked to dysregulated Wnt pathway activation, providing unique insights into the pathogenesis of TD.
Results
GWAS of TD
A total of 142 unrelated Japanese patients with TD (13 with thyroid aplasia, 44 with hypoplasia and 85 with ectopia, based on ultrasound or scintigraphy) were enrolled. Seven patients had extra-thyroidal manifestations, including congenital heart defect (N = 3), cleft palate (N = 2), tuberous sclerosis (N = 1) and multiple hemangiomas (N = 1). After genotype data acquisition, imputation and quality control (QC) processes, we assessed association using the genotype distribution of 140 cases and 8380 ethnicity-matched controls. By analyses of 482 564 genotyped single nucleotide polymorphisms (SNPs) and 2 962 984 imputed SNPs, we identified a single genomic region in 2q33.3 that reached genome-wide threshold for statistical significance (Fig. 1 and Supplementary Material, Fig. S1). The 72-kb genomic region (hg19 chr2: 208838738-208 910 964), consisting of 80 consecutive significant SNPs, was located in the intron and upstream sequences of PLEKHM3 (top genotyped SNP, rs9789446: P = 4.4 × 10−12) and corresponded to a linkage disequilibrium block (Fig. 1).
Replication
As a replication study, we genotyped the top SNP rs9789446 in 80 unrelated German TD patients including 41 individuals with thyroid aplasia, 8 with thyroid ectopia and 31 with thyroid hypoplasia. Statistical comparison of the genotype distribution of German patients and controls (N = 4339) confirmed that the rs9789446 genotype was deviated significantly (P = 0.0056), with a higher frequency of AA genotype in the patient group (Table 1).
rs9789446 genotypes in subjects and controls
| Stage (Site) | Cohort | rs9789446 | rs9789446 genotype (%) | χ2 test P values | ||
|---|---|---|---|---|---|---|
| A allele frequency | GG | GA | AA | |||
| Discovery (Japan) | Case | 0.629 | 16.4 | 41.4 | 42.2 | 4.4 × 10−12 |
| Control | 0.449 | 30.5 | 49.2 | 20.3 | ||
| Replication (Germany) | Case | 0.850 | 1.3 | 27.5 | 71.3 | 0.0056 |
| Control | 0.737 | 6.8 | 39.0 | 54.1 | ||
| Stage (Site) | Cohort | rs9789446 | rs9789446 genotype (%) | χ2 test P values | ||
|---|---|---|---|---|---|---|
| A allele frequency | GG | GA | AA | |||
| Discovery (Japan) | Case | 0.629 | 16.4 | 41.4 | 42.2 | 4.4 × 10−12 |
| Control | 0.449 | 30.5 | 49.2 | 20.3 | ||
| Replication (Germany) | Case | 0.850 | 1.3 | 27.5 | 71.3 | 0.0056 |
| Control | 0.737 | 6.8 | 39.0 | 54.1 | ||
rs9789446 genotypes in subjects and controls
| Stage (Site) | Cohort | rs9789446 | rs9789446 genotype (%) | χ2 test P values | ||
|---|---|---|---|---|---|---|
| A allele frequency | GG | GA | AA | |||
| Discovery (Japan) | Case | 0.629 | 16.4 | 41.4 | 42.2 | 4.4 × 10−12 |
| Control | 0.449 | 30.5 | 49.2 | 20.3 | ||
| Replication (Germany) | Case | 0.850 | 1.3 | 27.5 | 71.3 | 0.0056 |
| Control | 0.737 | 6.8 | 39.0 | 54.1 | ||
| Stage (Site) | Cohort | rs9789446 | rs9789446 genotype (%) | χ2 test P values | ||
|---|---|---|---|---|---|---|
| A allele frequency | GG | GA | AA | |||
| Discovery (Japan) | Case | 0.629 | 16.4 | 41.4 | 42.2 | 4.4 × 10−12 |
| Control | 0.449 | 30.5 | 49.2 | 20.3 | ||
| Replication (Germany) | Case | 0.850 | 1.3 | 27.5 | 71.3 | 0.0056 |
| Control | 0.737 | 6.8 | 39.0 | 54.1 | ||
Subgroup analysis by thyroid morphology
Using rs9789446 genotype data of discovery and replication datasets, we performed a subgroup analysis based on thyroid morphology (aplasia, hypoplasia or ectopia). Per-allele odds ratio (OR) based on additive model of logistic regression showed similar results in the two datasets: per-allele ORs were significantly increased for thyroid aplasia (3.17, 95% confidence interval 1.84–5.45, P < 0.0001) and thyroid ectopia (3.12, 95% confidence interval 2.24–4.34, P < 0.0001), but not increased for thyroid hypoplasia (1.23, 95% confidence interval 0.87–1.74, P = 0.25) (Fig. 2).
Subgroup analysis based on thyroid morphology. The forest plot shows associations of rs9789446 with TD, expressed as per-allele OR, across three subgroups defined by thyroid morphology (thyroid aplasia, hypoplasia and ectopia). The number of German subjects with ectopia was so small (N = 6) that calculation of OR was not applicable (NA). Boxes represent point estimates of effect. Lines represent 95% confidence interval.
Characterization of the disease-associated region in 2q33.3
The disease-associated region in 2q33.3 does not contain genes previously linked to thyroid differentiation and development (19). To interrogate biological functions of the genomic region, we leveraged publicly available datasets. First, we characterized epigenome signatures of the region in human thyroid by intersecting published single-cell open chromatin data of thyroid follicular cells (20) with mRNA expression and histone modification profiles of bulk human thyroid tissue samples (21). The 72-kb TD-associated region displays chromatin features of an active enhancer state including stretches of open chromatin with enriched H3K27ac and H3K4me1 chromatin immunoprecipitation signals (Fig. 3). Secondly, we analyzed chromatin interaction profiles obtained from chromosome conformation capture (Hi-C) studies of 27 human cell and tissue types (22) to explore possible long-range interactions with cis-regulatory elements of the TD-associated region. Our analyses of a broad panel of human tissues revealed abundant examples of long-range interactions of the TD-associated region with promoters of FZD5 and CCNYL1, which are located approximately 200-kb centromeric from the region (Fig. 3). In fact, these two genes were robustly expressed in human embryonic thyroid (Fig. 3 and Supplementary Material, Fig. S2) (23). Thirdly, we queried the Genotype-Tissue Expression (GTEx) V6 dataset (24) and found that 158 expression quantitative trait loci (eQTLs) were mapped to the 72-kb region. Notably, all these eQTLs were associated with expression levels of FZD5 and CCNYL1 (Fig. 4). In addition, among 44 tissues analyzed in the GTEx V6 dataset, regulatory associations were confined to adult thyroid tissue but absent for other tissues (Supplementary Material, Fig. S3). Specifically, the risk rs9789446 allele (A genotype) was associated with increased expression of FZD5 and CCNYL1 in thyroid tissue (Fig. 4).
Summary of gene regulation analyses in the 2q33.3 locus. The upper panel shows an overview of the 72-kb TD-associated region. Data of (i) histone modifications (H3K27ac and H3K4me1), (ii) the assay for transposase-accessible chromatin with high-throughput sequencing (ATAC-seq), which are associated with enhancer activity, (iii) capture Hi-C, visualizing long-range interactions involving promoters of FZD5 (blue) and CCNYL1 (red) and (iv) RNA expression are visualized.
Summary of eQTL analysis in the 2q33.3 locus. (A) The upper panel shows an overview of the 72-kb TD-associated region and adjacent regions. In the lower panel, statistically significant eQTLs (observed in 8 organs) registered in the GTEx V6 database are indicated by bars. No eQTL was observed in this region for 36 other organs, and these organs were omitted from this figure. Note that eQTLs located in the 72-kb locus were detected only in thyroid. (B) The allelic effect of rs9789446 on thyroid expression levels of FZD5 and CCNYL1 are shown by boxplots within violin plots. The data were obtained from the GTEx V6 dataset. Violin plot shows the density plot of the data on each side; the lower and upper borders of the box correspond to the first and third quartiles, respectively; the central line indicates the median.
Collectively, our analyses of epigenome, chromatin interaction, transcriptome and eQTLs strongly show that the TD-associated region contains cis-regulatory sequences contributing to the control of FZD5 and CCNYL1 expression in thyroid tissue. Although the distance between the TD-associated region and FZD5 exceeds 200 kb, promoter capture Hi-C data (22) imply a three-dimensional chromatin organization facilitating spatial juxtaposition of these regions and long-range regulatory chromatin interactions.
Discussion
In the present study, we conducted the first GWAS for TD and identified a single disease-associated region in 2q33.3. The locus is unique and has not been listed in GWAS catalogue (https://www.ebi.ac.uk/gwas/). We showed that the rs9789446 allele confers similar risk to thyroid aplasia and ectopia, but not to hypoplasia. This implies that thyroid aplasia and ectopia are in a continuous spectrum of defects with a shared pathogenetic perturbation in early organogenesis. Thyroid hypoplasia seems to be different, either in terms of a later timepoint of pathogenesis probably due to more diverse factors affecting thyroid cell proliferation and survival or in terms of a less stringent phenotypic separation from defects in genes of thyroid hormone synthesis with an anatomically normal gland.
Epidemiologically, TD is more prevalent in Caucasians, less prevalent in blacks, and Asians are in between (16). The frequency of the risk rs9789446 allele in the population is in a corresponding order (European 0.7394; East Asian 0.5611; African/African-American 0.3668, according to gnomAD v2.1.1) (25), which may explain at least in part the ethnic difference of TD.
A particularly interesting finding of our study is that the risk rs9789446 allele was associated with enhanced expression of two genes (FZD5 and CCNYL1) involved in Wnt pathway signaling. FDZ5 encodes for Frizzled 5, a Wnt receptor known to mediate canonical β-catenin-dependent as well as non-canonical (planar cell polarity) Wnt signaling (26,27). CCNYL1 encodes for cyclin-Y-like protein 1, which can enhance Wnt signaling through phosphorylation of the Wnt receptor LRP6 (low-density lipoprotein receptor-related protein 6) (28–30). Knowledge about the role of Wnt signaling during thyroid development is still fragmentary. In a transcriptomic analysis using three ectopic thyroid tissues removed as treatment for dysphagia, increased expression of genes constituting the Wnt pathway was observed (31). More recently, in vivo zebrafish studies provided evidence for disturbed thyroid anlage specification and thyroid primordium formation in response to experimentally enhanced Wnt signaling during early embryonic development (32,33). Under the assumption that tightly controlled level of Wnt signaling is a prerequisite of normal early thyroid development, we argue here that a dysregulation of Wnt receptor expression levels (FZD5) along with altered LRP6 competence (CCNYL1) due to the TD-predisposing risk allele rs9789446 provides a biologically plausible mechanism model to frame future studies on the pathogenetic mechanisms of TD.
In summary, we performed the first GWAS of TD and identified a single strong signal at 2q33.3 with per-allele OR exceeding 3.0. Despite the small sample size (N = 142), we were able to identify disease-associated genetic signals, reveal potential regulatory effects on gene expression and gain unique insights into a new disease mechanism for TD related to perturbed Wnt signaling. Our study exemplifies that even in rare diseases with a small sample size like TD GWAS may provide significant findings, given that the biological effect is strong enough, as it was shown already for more common malformations like non-syndromic orofacial clefts (14,15).
Materials and Methods
Subjects
This study was conducted as a collaborative effort to investigate the genetic basis of TD. A total of 142 genomic DNA samples (51 males and 91 females) extracted from peripheral leukocytes of Japanese patients with TD, who were followed up in 21 institutions in Japan, were collected at Department of Pediatrics, Keio University School of Medicine in the period from January 2006 to March 2018. Clinical assessment and recruitment for the study were done at each institution. Written informed consent for participating this study was obtained from the study subjects and/or their parents. The study protocol was approved by ethics committees of Keio University School of Medicine and National Center for Child Health and Development. Diagnosis of TD was based on the images of thyroid ultrasonography and/or 123I scintigraphy. Thyroid hypoplasia was defined by ultrasonographic findings: presence of a small-sized thyroid gland (thyroid width z score less than −2.0) (34) in a normal position. In patients with undetectable thyroid gland in normal position, differential diagnosis of aplasia and ectopia was primarily based on 123I scintigraphy findings. In patients who did not receive 123I scintigraphy, we classified the morphology subgroups based on serum thyroglobulin level (aplasia, <10 ng/ml; ectopia, ≥10 ng/ml). For all subjects, next generation sequencing-based comprehensive genetic screening was performed, as previously described (35), to confirm that none had disease-causing variants in FOXE1, NKX2-1, PAX8 or TSHR.
Genotype data acquisition and imputation
Genomic DNA was extracted from the 142 subjects using a standard procedure. In addition to the 142 samples, we analyzed 47 in-house QC samples, derived from individuals without TD, to find assay-specific errors. The 189 samples (TD 142; QC 47) were genotyped for 647 694 SNPs on autosomes and X chromosome using Japonica Array V2 (Toshiba, Tokyo, Japan). Genotype calling was conducted using Axiom Analysis Suite 1.1.1.66 (Thermo Fisher Scientific, Waltham, MA, U.S.A.). Two TD samples and one QC sample had low call rate (<0.97) and were omitted from the following analyses. The quality of information of each SNP was examined by the four criteria: (i) call rate >0.99 (14 067 SNPs excluded), (ii) minor allele frequency >0.05 (127 961 SNPs excluded), (iii) no significant deviation (P > 0.01) from Hardy–Weinberg equilibrium (5978 SNPs excluded) and (iv) no significant difference (P > 0.01) in genotype distribution between QC samples and 8.3KJPN (36) (17 124 SNPs excluded). The remaining 482 564 genotyped SNPs were used for GWAS. Untyped genotypes were imputed using IMPUTE2 and the 3.5KJPN data (37). Autosome SNPs were phased using SHAPEIT as follows: --burn 10, --prune 10 and --main 25. Phased genotypes were imputed with IMPUTE2 using the 3.5KJPN panel as follows: -Ne 2000, -k hap 1000, -burnin 15 and -iter 50. The quality of information of each imputed SNP was examined by the five criteria: (i) imputation quality (posterior genotype probability of >0.8), (ii) call rate >0.97, (iii) minor allele frequency >0.05, (iv) no significant deviation (P > 0.01) from Hardy–Weinberg equilibrium and (v) no significant difference (P > 0.01) in genotype distribution between QC samples and 8.3KJPN. A total of 2 962 984 imputed SNPs were subject to GWAS. As the control data, we retrieved genotype (diplotype) distribution data of 8.3KJPN, which are derived from 8380 healthy Japanese individuals (36), from the Tohoku Medical Megabank Organization. 8.3KJPN sequence data are based on whole genome sequencing, and these data have been confirmed to be concordant with Japonica Array V2 (37).
Replication
Samples used for the replication study were derived from 80 German patients with TD, including 41 with thyroid aplasia, 31 with hypoplasia and 8 with ectopia as determined by ultrasound, that were followed in one clinical center (Charite-Universitätsmedicine-Berlin, Germany). DNA was extracted from peripheral leukocytes of the study subjects. The SNP rs9789446 was genotyped with PCR-based Sanger sequencing. The study protocol was approved by ethics committee of Charité Universitätsmedizin Berlin (EA2/131/11). As the control sample for the German population, data of the population-based Heinz Nixdorf Recall Study (38) were used. The genotypes of rs9789446 were generated from genome-wide genotyping data (GSAMD24v1-0 array; Illumina Inc., San Diego, CA, U.S.A.) of 4349 subjects by imputation (info score = 0.996).
Statistical analyses
For GWAS, we calculated P values by χ2 test to examine single-marker genetic associations with TD. SNPs with P < 5 × 10−8 were considered of genome-wide significance. The value is based on a Bonferroni correction for multiple testing, assuming that 1 × 106 is greater than the number of linkage disequilibrium blocks (39). Corrected P = 0.05/(1 × 106 independent association tests) = 5 × 10−8.
For replication study of the genetic association, we also used χ2 test to calculate a P value. For the replication, P < 0.05 was set to a significance threshold. P values were calculated with Microsoft Excel 2013. For subgroup analysis, we evaluated the risk allele (rs9789446 A allele) with additive model logistic regression, where genotypes were regarded as a predictor variable (0 for GG, 1 for AG and 2 for AA) and per-allele OR was calculated. To perform joint analysis of Japanese and German datasets, ethnicity (Japanese or German) was added as a predictor variable. We used the glm function of R (version 4.1.0, https://www.R-project.org/) to calculate ORs with 95% confidence intervals. The data used for calculation of OR are available as Supplementary Material.
Bioinformatic analyses of the risk locus
A dataset of histone modifications (H3K27ac and H3K4me1) in adult thyroid was obtained from Canadian Epigenetics, Environment and Health Research Consortium Network (http://www.epigenomes.ca/) (21). Data of the assay for transposase-accessible chromatin sequencing were obtained from Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/; accession number GSE165659) (20). Promoter capture Hi-C data were obtained from 3D-Genome Interaction Viewer and database (http://www.3div.kr/capture_hic) (22). For evaluation of eQTLs, GTEx V6 dataset (24) was obtained from GTEx Portal (https://www.gtexportal.org/home/datasets). Transcriptomics data of embryonic tissues were retrieved from Expression atlas, University of Manchester (http://www.humandevelopmentalbiology. manchester.ac.uk/) (23).
Acknowledgements
The authors thank Naoko Amano, Makoto Anzo, Yuri Etani, Junpei Hamada, Takashi Hamajima, Yukihiro Hasegawa, Naomi Hatabu, Naoaki Hori, Shinobu Ida, Mikako Inokuch, Tsutomu Kamimaki, Yasusada Kawada, Nobutake Matsuo, Akira Ohtake, Makoto Ono, Goro Sasaki, Seiji Sato, Hironori Shibata, Eri Suzuki, Masaki Takagi and Noriyuki Takubo for providing genetic samples; Sabine Jyrch and Rita Oeltjen for technical assistance; Stefan Herms for genotype imputation to generate control data. A part of the data used in this study was provided by the Center for Genome Platform Project in the Tohoku Medical Megabank Organization.
Conflict of Interest statement. None declared.
Funding
JSPS KAKENHI (15K09630); The Japanese Society for Pediatric Endocrinology Future Development Grant supported by Novo Nordisk Pharma Ltd; The German Research Council (DFG) grant ‘beyond the exome’ (Forschergruppe 2841).
