Thyroid hormone-regulated gene expression in juvenile mouse liver: identification of thyroid response elements using microarray profiling and in silico analyses

The thyroid participates in the regulation of basic physiological processes by producing thyroid hormones (THs), which include thyroxine (T4) and triiodothyronine (T3). THs exert their effects on growth, development and metabolism of practically every cell and organ [1]. Their primary effect is the transcriptional regulation of target genes. This occurs when THs interact with TH receptors (TRs). Similar to other nuclear receptors, TR contains a DNA-binding domain that is capable of interacting with specific DNA sequences known as thyroid response elements (TREs). Typically, TREs are composed of two or more hexamer half-site sequences arranged in tandem array. TRs have the ability to bind to various imperfect TRE half-sites. The number of half-sites, the spacing between half-sites and their orientation are all features that can vary between TREs [2]. TRs can bind to TREs as monomers or homodimers, although it is thought that TRs interact as heterodimers with the retinoid × receptor (RXR) the majority of the time [3, 4].

Disruption of TH physiology during critically sensitive periods in development can lead to adverse outcomes. Studies have shown that severe TH insufficiencies lead to somatic and brain growth retardation [5]. Less severe and/or transient hypothyroidism has also been shown to cause adverse structural and functional effects [6, 7]. In humans, subclinical hypothyroidism during fetal development has been shown to result in reduced cognitive function [8].

Exposures to a broad array of substances, including both natural and synthetic chemicals, have been shown to alter TH physiology (see [9, 10] for review). Environmental contaminants such as perchlorates [11], polychlorinated biphenyls [12], bisphenol A [13], polybrominated diphenyl ethers [14], triclosan [15] and nitrate [16] have all been shown to have a negative effect on TH function, ranging from reducing circulating THs to altering TH-signalling.

Much effort has been expended to identify direct transcriptional targets of TH [17–19]. Identification of the thyroid regulated transcriptome will provide insights into the molecular impact of THs in directing tissue function/development. We have previously characterized TH-regulated global gene expression in the livers of juvenile mice [18]. This study evaluated transcriptional effects in animals in which TH production had been severely depressed through a relatively long-term treatment with a high concentration of an antithyroid drug (6-propyl thiouracil (PTU)). This model proved inefficient for the identification of genes directly regulated by TH as the lengthy period of treatment with antithyroid substance prevented a clear interpretation of the mode of action for altered genes. It was unclear if altered gene expression was caused by direct TH action, or due to downstream effects in the liver resulting from an altered developmental trajectory.

There is overwhelming evidence that environmental contaminants can act as endocrine disrupters with possible negative consequences for human health. However, major needs in the field include the identification of key initiating events leading to impaired development and tissue function. In the current study, we apply a well-validated animal model and robust microarray analysis to investigate hepatic transcriptional response to transient hyper- and hypothyroidism. The transient treatment time in combination with conditions of both hyper- and hypothyroidism provide a more efficient approach to identify direct hepatic targets of THs during liver development. In addition to identifying key genes that are directly controlled by THs, the work sheds light on the TRE sequences that direct TR binding. To identify TH targets, livers were collected from juvenile mice (postnatal day (PND) 15) rendered transiently hypo- or hyperthyroid. This developmental period corresponds to a dramatic increase in circulating THs. In PND 15 mice, circulating T4 is higher than at any other age [20] suggesting that chemically-induced disruption of TH levels should cause a marked response in transcription. Using DNA microarrays, we identify genes that are differentially expressed between control and TH-modulated animals, and provide supporting evidence for the presence of TREs in the regulatory regions of several of these genes. Collectively, the results provide an important addition to general knowledge of the genes that are directly regulated by THs and are useful for defining an improved TRE consensus sequence.

To identify direct transcriptional targets of TH and genes containing TREs, we produced global gene expression profiles in hypothyroid and hyperthyroid mouse models established by short-term exposure to MMI and perchlorate and/or TH injection. PND 15 livers were studied as animals of this age have circulating T4 levels higher than at any other age [20] and the liver is widely studied as a model of TH action. In the hypo pups TH levels were altered for three days, whereas in the hyperthyroid pups TH levels were altered for four hours immediately before tissue collection. Concentrations of anti-thyroid agents used in the hypo exposure and concentrations of injected THs were not unusually high when compared to other similar studies [17, 33]. Several hundred genes were responsive to this short-term treatment, with fold-changes as high as 13.8. However, the vast majority of these genes had relatively small fold-changes (i.e., < 3-fold). When filtered on a 2.0 fold-change cutoff, only 31 genes in females and 27 genes in males were perturbed, with an overlap of 14 genes.

PCA confirmed that the microarray data clustered by sex and by treatment type suggesting that there are significant differences at the level of gene expression between i) male and female pups and ii) control, hyper, hypo and hypo+ pups. GO analysis identified multiple pathways involved in oxidative stress and xenobiotic metabolism/signalling, as well as confirming the enrichment of genes in TR/RXR activation (Additional File 3). The focus of the present research was to identify direct TH-target genes, rather than conduct a detailed mechanistic analysis. Thus, the data were mined specifically to identify direct TH-regulation and search promoters for putative TREs.

As previously mentioned, studies characterizing TH hepatic gene regulation have been carried out by others [33–36] and by our group [18]. Although these studies were carried out using different experimental designs, animal models, and developmental stages, some important similarities in molecular functions and processes emerge including genes involved in TR and RXR activation (see Additional File 3). In addition, several specific genes previously shown to be involved in TH hepatic gene regulation, such as Me1, Thrsp and Dio1 were also found in the current study. However, our study design is unique relative to the other work described above because it is a short-term transient perturbation in TH that queries a specific developmental period to identify immediately responsive genes. The overarching goal of the work is to use a robust microarray analysis of gene expression profiling coupled with a novel TRE search algorithm to identify putative TREs. Thus, a detailed functional analysis of the perturbed genes in this study compared to the other studies has not been carried out.

To identify genes under direct TH-regulation we examined the expression of a subset of responsive genes in the hypo+ and hypo++ groups in detail. These groups were rendered hypothyroid for the majority of the 3 day exposure, and then received an i.p. TH injection four hours before decapitation and tissue collection. In other words, only a four hour period was allowed for transcriptional response to the TH surge. Where the direction of the transcriptional fold-change in the hypo+ group was the same as in the hyper group, but the opposite of the hypo group, we predicted that the gene in question was under direct TH-regulation. Transcript levels of some of these genes were analysed by RT-qPCR in the hypo+ and hypo++ groups (Table 3). RT-qPCR analysis confirmed a significant change in the hypo+ group for four out of the nine genes relative to vehicle controls. Analysis of the hypo++ group, which registered a much higher level of circulating serum T4 compared to hypo+, showed that seven out of the nine genes had significant fold-changes relative to control. Validation of the transcriptional response in the hypo+ and hypo++ groups help to further support our assumption that the identified genes (Table 1) are under direct TH-regulation. None of the genes identified as being directly regulated by THs have characterized TREs. Although specific genes that are commonly considered to be directly regulated by TH, such as Me1, Thrsp or Dio1, were confirmed to be differentially regulated in our experiment, they did not meet the stringent requirements set for this particular evaluation (i.e., our list of potentially directly regulated genes; Table 1). As previously described, genes were considered to be directly regulated by TH if they followed a specific expression pattern. Only genes exhibiting fold-changes with FDR-adjusted p ≤ 0.05 were considered for this analysis. However, the RT-qPCR analysis demonstrated that known TH responsive genes were responding as predicted in the animal model used (Table 2).

We retrieved the relevant DNA sequences and identified potential TREs within the promoter regions of the responsive genes (i.e., the genes in Table 1). We assumed that genes under direct TH-regulation should contain functional TREs in their promoter regions. For practical reasons, our search was limited to the immediate promoter region (from 8 kb upstream to 2 kb downstream of the TSS) even though not all response elements are located in this area. Some studies have shown that transcription factor binding sites can be exceptionally far from the chromosomal location of the genes they regulate. For example, some estrogen response elements (EREs) can be found 50 kb away from of the TSS [37], whereas others can be more than 100 kb away from the TSS [38]. In addition, interchromosomal regulatory interactions have also been documented [39]. Our gene expression profiling analysis identified 28 genes that appeared to be under direct TH-regulation. Within these 28 genes, bioinformatics mining for three different types of TREs (DR4, IR0 and ER6) revealed 196 candidate TREs. To help identify true TREs, we carried out a local alignment of promoters of these mouse genes with the promoter regions of Rattus norvegicus and Homo sapien. In the -8 kb to +2 kb promoter regions of these genes we were able to identified 33 TREs in 13 different genes. Thus, within the region examined approximately 46% of the genes possessed potentially conserved TREs that fall within the established half-site criteria.

Conservation of a stretch of DNA between all three species provides evidence of the selective importance of response element functionality. Two previously characterized TREs are presented at the bottom of Table 4. The TRE near Klf9 has been characterized in mice, rats, and humans [40]. Comparison of the TREs within the orthologous Klf9 gene of the three species shows that the TREs are within 2 kb of each other relative to the TSS. The TRE located in the Mbp gene promoter region has been characterized in mice and rats [41, 42]. The TREs for Mbp for both species are within 2 kb of each other relative to the TSS and show very good half-site sequence conservation. Conservation of response elements between species has been previously observed. Similar to the current study, Bourdeau et al. identified EREs in 660 different pairs of orthologs between mice and humans [43]. EREs identified in proximity of these orthologs were within 2 kb of their respective TSSs. Some EREs that were conserved between mice and humans were perfect matches, whereas others had multiple mismatches in their two half-site. EREs are very similar to TREs; they are both made up of AGGTCA half-sites. The classic ERE organization is an IR3, where the two half-sites are arranged in a palindromic organization separated by a 3 bp spacer. The classic TRE, on the other hand, is a DR4 with the two tandem repeated half-sites separated by a 4 bp spacer [44], although IR0, ER6 and other TREs have also been characterized. The similarities and differences that exist between EREs and TREs suggest that the underlying mechanisms that confer response element specificity are much more complicated than was once thought. A recent paper by Phan et al. points out some interesting findings about nuclear receptor DNA recognition specificity [45]. The authors looked at the half-site recognition for retinoic acid receptor (RAR) and TR, which both bind to the AGGTCA half-site sequences. It was once thought that the spacing between half-sites was the only aspect that conferred specificity [46–48] although the authors discovered that there were other properties that played an important role, which included naturally occurring non-consensus half-sites, flanking sequences, and auxiliary proteins produced by the cell. These additional properties make the identification of response elements even more complex, especially considering the high degree of degeneracy of TREs; nevertheless, they are important to take into consideration. It is also important to point out that our model did not consider flanking sequences and that the PWM was built with "known" TREs that biased towards the AGGTCA classic DR4 organization. We are hopeful that future response element search tools will be more flexible and take additional details such as those described above into consideration.

Previous studies have primarily used search algorithms that allow very little variation from the classic AGGTCA half-sites. Search algorithms often only allow substitution of a position towards one specific base pair, for example, position one of the half-site would only allow for an A or a G. This is a very restricted approach since TREs can be very degenerated and, based on previously characterized TREs, position one (as well as other positions) could potentially have any of the four nucleotides. Moreover, previous search algorithms often did not allow substitutions to occur at all half-site positions and often limited the number of "mismatches" when compared to the classic AGGTCA TRE half-site. Our TRE search algorithm was much more flexible when compared to these other approaches, as the PWM scored all possible TREs and applied a score based on the probability for a given base to be located at any position of the half-site based on known functional TREs. A cut-off value was applied to this criterion to minimize false-positives. Since our approach allowed for a great deal of half-site sequence variation, a filtering or refining step was needed to reduce false-positives. We found that a cross-species comparison was an effective method to reduce the amount of false-positives, and also increases the potential utility of the new TREs for future research as they may be relevant in multiple species. Our method allowed for a larger amount of variability and as a result we identified new candidate TREs that contain some features of the classic TRE, but are substantially different.

ChIP-PCR was used to validate some of the TREs identified by the bioinformatics search. Six candidate TREs were considered for this validation work. Four showed enrichment by TR-immunoprecipitation, whereas the two others did not show enrichment. This suggests that the FDR is approximately 33%, but our sample size for this estimate is quite low. The lack of identification of TREs in certain genes could be attributed to various factors. TREs could be present in the promoter regions of genes thought to be directly regulated by TH but are undetected by the bioinformatics tools currently available (i.e., do not exhibit the classic characteristics). Functional TREs have been described with a high level of degeneration, multiple spacer sizes and various half-site organizations [18, 49–52] with a structure divergent from the models on which our identification algorithm was based. Consequently our model may have been too conservative to identify other functional TREs in these genes. Alternatively, TREs may truly be absent from the promoter regions of these genes. Expression may be tied to TH action through intermediate regulatory mechanisms. Non-genomic actions of TH have also been characterized [53], by which TH activation of plasma membrane receptors induces signal transduction pathways leading to various genomic or cellular responses, although some claim that the non-genomic effects of THs do not play a significant role during vertebrate development [54]. MicroRNAs (miRNAs) could also potentially regulate TH-mediated mRNA expression. We recently identified significant alterations in 40 different miRNAs in the livers of PND 15 hypothyroid mice [32]. Lastly, within our list of significantly altered genes (FDR adjusted p ≤ 0.05), several transcription factors are present, including Rarβ, Esr1 and Nr3c1 (see Additional File 1 and 2). Thus, genes present within our list may be under the control of other transcription factors that are directly regulated by THs. For example, Esr1 expression is decreased in hypothyroid female mice. This could lead to many downstream transcriptional effects since Esr1 is a transcription factor that has been shown to interact with 18 different nuclear receptors (Nuclear Receptor Signaling Atlas, http://www.nursa.org).

ChIP-PCR revealed significant enrichment for some of the stretches of DNA thought to be involved in TR interactions. TRs are thought to be bound in both presence and absence of THs, thus ChIP-PCR was only performed with euthyroid animals. Importantly, the PCR primers targeted a DNA fragment with an average size of 168 bp which includes the putative TRE. In total, four TRE sites were validated by ChIP-PCR. The TREs identified in the promoter regions of Tor1a and Slc25a45 correspond to a DR4 organization. The TRE identified in the promoter region of Hectd3 corresponds to an ER6 organization, whereas the TRE identified in the promoter region of 2310003H01Rik corresponds to an IR0 organization. The genomic locations of these newly identified TREs all fall in upstream promoter regions, -2.3, -4.6, -7.4 and -2.3 kb from their respective TSS. Most TREs identified to date in the mouse genome are located in upstream promoter regions, and more then 70% are within the 0 to -2.5 kb upstream region. The TRE that is the furthest from a TSS that has been characterized to date in mice is in the promoter region Klf9, and is -3.8 kb from its TSS. This information is based on the 14 mouse TREs that we were able to find in the literature (list available on request). Since 50% of the TREs we have characterized here are outside these bounds, our findings suggest that TREs are not necessarily found more often within the 2 or 3 kb upstream promoter region. The TRE associated with 2310003H01Rik was not validated by EMSA. This could be due to the fact that recombinant purified protein samples were used, as opposed to tissue or cellular extracts. Extracts would include other transcription factors, such as TR's known heterodimer partner RXR, which may be required for the DNA-protein interaction to occur at a detectable level. In contrast, the TREs associated with Slc25a45, Hectd3 and Tor1a were validated by EMSA, each showing a shift in the presence of TR proteins and a decrease in detection when unlabelled probe or DR4 unlabelled probe was added. The addition of unlabelled probe demonstrates the specificity of the shifted band, whereas the addition of the unlabelled DR4 probe demonstrates that the shifted band was caused by an interaction with the TR protein, since a supershift using an antibody against TR was shown for the DR4 positive control probe. A correlation between the ChIP-PCR and EMSA results was observed. Scl25a45 has the strongest enrichment measured by ChIP-PCR, followed by Hectd3 and Tor1a. Slc25a45 again showed the strongest decrease in signal intensity when unlabelled specific probe was added in the EMSA, followed by Hectd3 and Tor1a.

Slc25a45 is involved in the transportation of molecules across the mitochondrial membrane [55]. THs have been shown to induce mitochondrial biogenesis and enhance ATP production [56]. The identification of a TRE in the vicinity of the TSS of this transporter could be an indication of its involvement in TH-dependant mitochondrial biogenesis. Hectd3 is a ubiquitin ligase. In humans, it has been shown to directly bind TARA, a guanine nucleotide exchange factor involved in regulating actin cytoskeletal reorganization, cell mobility, and cell growth [57]. THs have also been shown to play a role in cytoskeletal protein regulation [58]. Thus, a TRE in the promoter region of Hectd3 could link this gene to cytoskeletal regulation via THs. In mice 2310003H01Rik is the predicted ortholog of FAAP100, Fanconi anemia-associated protein, 100 kDa [59]. FAAP100 is involved in the Fanconi anemia (FA) core complex, which plays a role in the DNA damage response network [59]. FA is a genetic disease that results from defects in proteins involved in DNA repair. People affected by FA have physical anomalies, short stature, and are predisposed to cancer arising from chromosomal instability [60]. Children with FA have a high risk of endocrine abnormalities including hypothyroidism [61]. A recent study found that TH therapy could improve the linear growth of children suffering from FA [60]. The identification of a TRE in the promoter region of 2310003H01Rik could shed light on the underlying mechanisms leading to hypothyroidism in FA patients. Tor1a is an adenosine triphosphatase. Mutation of this gene has been linked to early onset dystonia and Parkinsonism [62]. Dystonia is a disorder characterized by sustained muscle contractions often causing twitching and repetitive movements or abnormal posture. THs are important for proper muscle development and maintenance [63]. Thus, the presence of a TRE in the promoter region of Tor1a could suggest a potential role for TH-Tor1a interaction in muscle development.

Using transcriptional profiling analysis we were been able to identify genes suspected of being directly regulated by THs in male and female mice. Bioinformatics mining of the promoter regions of these genes revealed 33 candidate mouse TREs potentially conserved between rats and humans. ChIP-PCR and EMSAs were employed to validate the mouse TREs. We provide evidence of four new TREs in the promoter regions of Tor1a, Slc25a45, Hectd3 and 2310003H01Rik. The results provide data to develop a stronger model for the TRE sequence motifs that direct TR binding, identify key genes that may be important for TH mediated effects, and help to further characterize the mechanisms by which TH directly regulate gene expression.