Comprehensive analysis of our GLIS3 ChIP-Seq data of the mouse thyroid gland [19] by HOMER indicated that binding motifs of the NKX, FOX, and PAX TF family are frequently localized near GLIS binding sites (GLISBS) (Fig. 1A). Since NKX2.1, FOXE1, and PAX8 have well-established roles in the regulation of gene expression in thyroid follicular cells [12, 13, 42, 43], we hypothesized that GLIS3 regulates transcription of a subset of target genes in these cells in coordination with these TFs.

Global analysis of GLIS3, NKX2.1, and PAX8 genomic binding in mouse thyroid glands. A–C HOMER analysis, heatmap, and ChIP-Seq read density of GLIS3 (A), NKX2.1 (B), and PAX8 (C) binding data. Heatmap of the 2 kb region is centered on each of the binding peaks identified. D Genomic context of the GLIS3, NKX2.1, and PAX8 peaks within the whole mouse genome (mm10)

To obtain support for this hypothesis, we performed NKX2.1 and PAX8 ChIP-Seq analyses with thyroid glands from mice fed a low iodine diet (LID), in which TSH blood levels are highly elevated [44]. As far as we know, this is the first time ChIP-Seq analyses with thyroid glands and endogenous PAX8 and NKX2.1 have been performed. These analyses identified 29,464 NKX2.1 and 41,044 PAX8 binding peaks. De novo motif analysis of the NKX2.1-enriched sequences identified, in addition to the NKX binding motif, consensus binding sites for members GLIS, PAX, and FOX families were among the top motifs (Fig. 1B). Motif analysis of PAX8 ChIP-Seq data identified a PAX binding sequence as the top binding motif together with binding motifs for NKX, GLIS, and FOX family members (Fig. 1C). NKX2.1 and PAX8 binding loci were most highly enriched within the gene body and intergenic regions, while 5–8% were localized within 1 kb upstream of TSS (TSS proximal) (Fig. 1D).

Comparison of GLIS3, PAX8 and NKX2.1 binding revealed substantial overlaps between GLIS3 binding loci and those of NKX2.1 and PAX8 (Fig. 2A, B). The majority of GLIS3 binding loci contained both NKX2.1 and PAX8 binding regions (referred to as G+N+P+). Subsets of GLIS3 binding regions overlapped with those of either NKX2.1 or PAX8 (G+N+P− or G+N−P+, respectively), while some of the GLIS3, NKX2.1 and PAX8 binding regions did not exhibit any overlap (G+N−P−, G−N+P−, and G−N−P+, respectively). These observations are consistent with the concept that the transcription of subsets of GLIS3 target genes are regulated in coordination with PAX8 and NKX2.1, and that some genes are regulated by only one or two of the three TFs.

GLIS3, NKX2.1 and PAX8 binding to the mouse thyroid gland genome partially overlap. A The number of collapsed binding regions that are positive for GLIS3, NKX2.1, and PAX8 signal, and the percent overlap between them are indicated. B Heatmap showing overlap between GLIS3, NKX2.1, and PAX8 binding. Heatmap of the 2 kb region centered on each of the binding regions with ChIP-seq signal normalized to 10 million reads. C Venn diagram showing the number of target genes with GLIS3, NKX2.1, and/or PAX8 binding regions. D Venn diagram showing the overlap of G+N+P+ genes with genes up- or down-regulated in the thyroid gland of Glis3KO-LID mice compared to that of WT-LID mice. E KEGG analysis of the 4502 G+N+P+ genes

Although enhancers can reside within the gene body and very distant from TSSs, we limited our analysis to genes with binding peaks within 5 kb up- or downstream from the TSS. This analysis identified binding of GLIS3, NKX2.1 and PAX8 near 5240, 8975 and 8832 genes, respectively (Fig. 2C). The majority (85.9%) of the GLIS3-bound genes shared binding both with NKX2.1 and PAX8 (G+N+P+), 9.9% with NKX2.1 only (G+N+P−) and 2.4% with PAX8 only (G+N−P+), while 1.8% bound GLIS3 only (G+N−P−) (Fig. 2C). A summary of GLIS3, NKX2.1, and PAX8 bound genes is presented in Additional file 6: Table S2A.

We previously reported that loss of GLIS3 function particularly suppresses the expression of a subset of genes that are required for TH biosynthesis, induced when mice are fed an LID or known to be induced by TSH [11, 19, 45, 46]. Therefore, we were interested in determining which of these differentially expressed genes were regulated by GLIS3 in coordination with PAX8 and/or NKX2.1. Among the 4502 G+N+P+ genes (Fig. 2C, D), 345 and 196 genes were, respectively, down- or up-regulated in Glis3KO mice fed an LID (Glis3KO-LID) compared to WT mice fed an LID (WT-LID) thyroid (Additional file 6: Table S2B, C). KEGG pathway analysis of the 345 down-regulated G+N+P+ genes identified TH biosynthesis as the top pathway (Fig. 2E). No pathway was found to be significantly associated with G+N+P+ up-regulated genes. Table 1 shows GLIS3, PAX8, and NKX2.1 binding to several gene clusters in relation to their repression in Glis3KO-LID thyroid and their induction in WT-LID thyroid. This comparison shows that GLIS3, PAX8, and NKX2.1 bound to many genes critical for TH biosynthesis, including Slc5a5, Slc26a4, Tpo, Tg, Iyd, and Duoxa2, and several genes known to be induced by TSH, such as Adm2, Sod3, Dio1, and Cdh16 [45,46,47,48,49]. The expression of several of these genes (e.g., Slc5a5, Slc26a4, Adm2, Sod3, Dio1, Cdh16) was significantly repressed in Glis3KO-LID mice compared to WT-LID and induced in the thyroid of WT-LID mice compared to WT-ND (Table 1). Other thyroid genes, including Tpo, Iyd, and Duoxa2, were induced in WT-LID thyroid, but not significantly affected by the loss of GLIS3 function, while the expression of certain GLIS3, PAX8, and NKX2.1 bound genes, including Tg, Txnrd1, and Duox2, were not significantly altered in WT-LID thyroid nor suppressed in Glis3KO-LID (Table 1).

The genome browser tracks in Fig. 3A indicate the shared locations of the binding of endogenous PAX8, NKX2.1, and GLIS3 in several genes critical for TH biosynthesis, including Slc5a5, Slc26a4, Duoxa2, Iyd, Tpo, and Tg, and Slc16a2. In several genes (e.g., Slc5a5, Slc26a4, Tpo, Tg), PAX8, NKX2.1, and GLIS3 bound within the same region of the proximal promoter. In Slc16a2 and Iyd only NKX2.1 and GLIS3 shared binding to the proximal promoter region, while NKX2.1 and PAX8 bound to Tshr, but not GLIS3. In several genes (e.g., Tpo, Tg, Slc26a4, Iyd, Slc16a2) binding of these 3 TFs overlapped in more than one region suggesting that their transcription may be controlled by multiple regulatory regions (Fig. 3A). These observations support the hypothesis that GLIS3 regulates gene transcription in coordination with several other thyroid TFs. In the case of Duoxa2, the GLIS3/PAX8/NKX2.1 binding region is in intron 1 of Duoxa2, which is within a 35 kb region on mouse chromosome 2 that also encompasses Duox2 and Duoxa1 (Fig. 3A). We cannot rule out that his enhancer region might play a role in the transcriptional regulation of all three genes.

GLIS3, NKX2.1, and PAX8 share binding loci within the regulatory regions of several TH biosynthetic genes. A Colocalization of GLIS3, NKX2.1, and/or PAX8 ChIP-seq loci in genes critical for TH biosynthesis. The NUE region in Slc5a5 (Nis) is indicated. B, C NKX2.1 and PAX8 analysis at Cdh16 (− 0.2 kb), Slc5a5 (− 2.8 kb), and Tpo (− 0.1 kb) with thyroid glands from WT and Glis3KO mice. Binding to Gapdh and Tpo (− 1.9 kb) served as negative controls

The transcriptional regulation of Slc5a5 has been extensively studied in thyroid follicular cell lines and reported to be controlled by the proximal promoter and a region, referred to as Nis upstream enhancer (NUE, located − 2.8 kb from the TSS) that has been reported to bind several TFs, including PAX8 and NKX2.1 [13, 14, 43, 50,51,52,53,54,55]. The Slc5a5 genome browser tracks show the localization of binding peaks for endogenous GLIS3, PAX8, and NKX2.1 within the NUE region, whereas no major binding was observed within the proximal promoter region (Fig. 3A). These observations are consistent with the view that NUE is a major enhancer region driving the activation Slc5a5 transcription by GLIS3 in the thyroid of mice fed an LID. As Scl5a5 is one of the genes most strongly regulated by GLIS3 (Table 1), this raised the question whether GLIS3 binding was required for the binding of PAX8 and NKX2.1 to the NUE region. To investigate this, we compared NKX2.1 and PAX8 binding to the NUE region in thyroids isolated from WT and Glis3KO mice. ChIP Q-PCR analysis demonstrated that NKX2.1 and PAX8 occupancy at the NUE region was not significantly affected by the absence of GLIS3 (Fig. 3B, C) indicating that GLIS3 is not required for PAX8 or NKX2.1 binding to the NUE regulatory region. The lack of GLIS3 also did not significantly affect PAX8 or NKX2.1 binding to the proximal promoters of Cdh16 and Tpo (region − 0.1 kb upstream from TSS). No significant binding of PAX8 or NKX2.1 was observed to Gapdh or the − 1.9 kb upstream region of Tpo, which served as negative controls (Fig. 3B, C). ChIP-Seq analysis of several epigenetic markers showed no significant differences in the level of H3K4me3 and H3K27me3 signals in GLIS3-bound genes that were differentially expressed between Glis3KO-LID and WT-LID thyroid glands (Additional file 1: Fig. S1A). Moreover, little difference in H3K4me3 and H3K27me3 signals was observed at the NUE region of Slc5a5 between WT and Glis3KO thyroid glands, although some small changes in these signals were observed at its proximal promoter (Additional file 1: Fig. S1B). Together, these observations suggest that in the absence of GLIS3 genomic regions, such as NUE, remain accessible and in an active (open) state.

Genome browser tracks in Fig. 4A show the overlap between the binding of endogenous PAX8, NKX2.1, and GLIS3 to regulatory regions of several other thyroid genes reported to be induced by TSH, induced in WT-LID, and repressed in Glis3KO-LID thyroids, including Adm2, Sod3, Dio1, and Cdh16 [45,46,47,48].

Genome browser tracks of several genes showing overlap of GLIS3, NKX2.1, and PAX8 binding loci in mouse thyroid gland. A Genes known to be induced by TSH. B Collagen and chemokine genes

GLIS3, NKX2.1, and PAX8 binding to extracellular matrix (ECM) and inflammatory genes.

Next, we examined whether any differentially expressed ECM and inflammatory genes are co-regulated by GLIS3, PAX8, and/or NKX2.1. In contrast to TH biosynthetic genes, relatively few ECM and inflammatory genes (e.g., Col18a1, Ccl9) showed binding of all 3 TFs, while a few genes (e.g., Col4a1, Col16a1) bound PAX8 and NKX2.1, but not GLIS3 (Table 1). Of course, we cannot rule out that these TFs might regulate some of these genes by binding distant enhancers. The genome browser tracks in Fig. 4B show overlaps between GLIS3, PAX8, and NKX2.1 binding regions in Col18a1, Col5a3, Col15a1, Col4a2, and Ccl9. The transcription of other differentially expressed ECM and inflammatory genes, such as Ccl7, that do not show GLIS3, PAX8, and NKX2.1 binding, are likely regulated by other TFs.

Our ChIP-Seq analysis further identified binding peaks of all three TFs within the same genomic region(s) of Glis3, Pax8, and Nkx2.1, as well as two other thyroid TF genes, FoxE1 and Hhex (Additional file 2: Fig. S2) suggesting transcriptional regulation of each other consistent with previous observations [3, 4, 56].

Since HOMER also identified FOX binding motifs near GLISBS (Fig. 1A), we were interested in examining co-regulation of TH biosynthetic genes by GLIS3 and FOXE1. Since several attempts to perform FOXE1 ChIP-Seq analysis in mouse thyroid glands were unsuccessful, we performed FOXE1, NKX2.1, and GLIS3-HA ChIP-Seq analysis in rat thyrocyte PCCl3 cells and compared these data with that of PAX8 in PCCl3 cells reported previously [43]. Nuclear expression of GLIS3-HA and endogenous FOXE1 in PCCl3 was confirmed by immunofluorescence staining (Additional file 3: Fig. S3). Consistent with our de novo motif analysis of the GLIS3 ChIP-Seq data from the mouse thyroid gland, analysis of the GLIS3-enriched sequences identified GLISBS as the top motif as well as consensus binding motifs for NKX, FOX, and PAX family members (Fig. 5A). Similarly, de novo motif analysis of the NKX2.1, PAX8 and FOXE1-enriched sequences identified, in addition to their own consensus binding motif, binding motifs of the three other TFs (Fig. 5B–D). The genomic contexts of the GLIS3, FOXE1, PAX8, and NKX2.1 binding peaks are shown in Fig. 5E. Comparison of GLIS3, PAX8, NKX2.1, and FOXE1 binding revealed substantial overlaps between the GLIS3 binding loci and the other three TFs (Fig. 6A, B). Most FOXE1-binding loci (74.8%) were in the proximity of GLIS3-binding loci. The different clusters of GLIS3, NKX2.1, PAX8 and/or FOXE1 bound genes (based on regions within 5 kb from the TSS) are presented in Additional file 7: Table S3. The data showed that GLIS3, NKX2.1, PAX8, and FOXE1 (G+N+P+F+) were bound near each other in 3118 genes, including several TH biosynthetic and TSH-responsive genes [e.g., Scl5a5 (NUE region), Duoxa2, Duox2, Cdh16, as well as Pax8, Nkx2.1, Hhex, and Foxe1], but not several other G+N+P+ thyroid genes (e.g., Slc16a2, Tshr, and Dio1). The genome browser tracks display the co-localization of GLIS3, NKX2.1, PAX8, and FOXE1 binding peaks in Scl5a5, Cdh16, Duoxa2, and Adm2 (Fig. 6C) and Glis3, Nkx2.1, Pax8, and Foxe1 (Additional file 4: Fig. S4) in PCCl3 cells. FOXE1 binding did not overlap with GLIS3, PAX8, and NKX2.1 binding at the proximal promoter of Adm2 but did show an overlap within a downstream intergenic region. The GLIS3, NKX2.1, and PAX8 binding patterns in PCCl3 were very similar to those as observed in mouse thyroid (Fig. 3A; Additional file 2: Fig. S2). However, binding of GLIS3, NKX2.1, and PAX8 in PCCl3 cells did not always match those observed in the mouse thyroid gland, including their binding to Tpo, Glis3, and Iyd (Table 2). This might be due to epigenomic differences between mouse thyroid follicular cells in vivo and immortalized rat thyrocyte PCCl3 cells, for which the gene expression profile is likely different from thyroid follicular cells in vivo.

Global analysis of GLIS3, NKX2.1, PAX8, and FOXE1 genomic binding in rat thyroid follicular PCCl3 cells. A–D HOMER analysis, heatmap, and ChIP-Seq read density of GLIS3 (A), PAX8 (B), FOXE1 (C), and (D) NKX2.1 binding data. E Genomic context of the GLIS3, PAX8, FOXE1, and NKX2.1 peaks within the whole rat genome (rn6)

GLIS3, NKX2.1, PAX8, and FOXE1 binding to the PCCl3 genome partially overlap. A Percent NKX2.1, PAX8, and FOXE1 binding loci overlapping with those of GLIS3. B Heatmaps of the 2 kb region centered on each of the binding regions with ChIP-seq signal normalized to 10 million reads shows overlaps between GLIS3, NKX2.1, PAX8, and FOXE1 binding. C Genome browser tracks of several thyroid genes showing overlap of GLIS3, NKX2.1, PAX8, and FOXE1 binding loci in PCCl3 cells