BACKGROUND
the thyroid participates in the regulation of basic physiological processes by producing thyroid hormones , which include thyroxine  and triiodothyronine .
ths exert their effects on growth, development and metabolism of practically every cell and organ  <cit> .
their primary effect is the transcriptional regulation of target genes.
this occurs when ths interact with th receptors .
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 .
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  <cit> .
trs can bind to tres as monomers or homodimers, although it is thought that trs interact as heterodimers with the retinoid × receptor  the majority of the time  <cit> .
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  <cit> .
less severe and/or transient hypothyroidism has also been shown to cause adverse structural and functional effects  <cit> .
in humans, subclinical hypothyroidism during fetal development has been shown to result in reduced cognitive function  <cit> .
exposures to a broad array of substances, including both natural and synthetic chemicals, have been shown to alter th physiology .
environmental contaminants such as perchlorates  <cit> , polychlorinated biphenyls  <cit> , bisphenol a  <cit> , polybrominated diphenyl ethers  <cit> , triclosan  <cit>  and nitrate  <cit>  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  <cit> .
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  <cit> .
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 ).
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  15) rendered transiently hypo- or hyperthyroid.
this developmental period corresponds to a dramatic increase in circulating ths.
in pnd  <dig> mice, circulating t <dig> is higher than at any other age  <cit>  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.
methods
animals and exposures
all animal care and handling was in accordance with canadian council for animal care guidelines and was reviewed by the health canada animal care committee prior to commencement of the study.
c57bl/ <dig> mice were purchased from charles river  and were housed in hanging polycarbonate cages under a 12: <dig> hr light-dark cycle at 23°c.
animals were provided with shelters, nesting material, food ad libitum  and drinking water containing 1%  sucrose .
mice were acclimated for  <dig> days.
breeding was accomplished by transferring two sexually mature female mice  into the home cage of a sexually mature  male.
after four nights of co-housing, each female was transferred to a separate cage.
dams were allowed to litter naturally and pup numbers were not adjusted.
to produce hypothyroid  pups, dams were provided with drinking water that contained methimizole  and perchlorate  for three days .
to produce hyperthyroid  pups, intraperitoneal injections  of ths  were administered to pups on pnd  <dig>  four hours before decapitation and tissue collection.
for the hypothyroid/replacement group ; dams were provided with drinking water that contained mmi  and perchlorate  for three days .
pups then received i.p.
injections of ths  on pnd  <dig>  four hours before decapitation and tissue collection.
the control and hypo groups received saline i.p.
injections four hours before decapitation and tissue collection.
tissue collection, rna extraction and purification
on pnd  <dig>  pups were sacrificed by decapitation.
serum prepared from trunk blood using serum separator tubes  was used for t <dig> analyses by radioimmunoassay  kits  as per the manufacturer's instructions.
the hypo + group was subdivided into two sub groups - hereafter called hypo+ and hypo++ - based on measured serum t <dig> levels.
animals allotted to the hypo++ group had serum t <dig> levels that approximated those seen in the hyperthyroid group while those that remained as hypo+ had t <dig> levels near or slightly in excess of control levels.
pup livers were rapidly dissected and flash frozen in liquid nitrogen.
total rna was extracted from liver samples with trizol reagent  followed by rneasy mini kit  clean-up according to the manufacturer's instructions.
rna quality assessment was determined by nanodrop  and agilent  <dig> bioanalyzer and rna  <dig> nanolab chip kit .
all samples had 260/ <dig> ratios over  <dig> , and rna integrity numbers over  <dig> .
microarray hybridization
hepatic gene expression was assessed by hybridizing samples to agilent  <dig> ×  <dig> k whole genome microarrays  using a reference design .
briefly,  <dig> ng of total rna from the liver was used to synthesize double-stranded cdna and cyanine labelled crna.
experimental samples were labelled with cyanine 5-ctp and reference rna  with cyanine 3-ctp  according to the manufacturer's instructions .
cyanine labelled crna targets were in vitro transcribed using t <dig> rna polymerase and purified using rneasy mini kit .
experimental and reference samples  were hybridized to arrays at 60°c for  <dig> hours.
slides were washed and scanned on an agilent microarray scanner .
the data were acquired using agilent feature extraction software version  <dig> .
<dig> .
prior to statistical analysis, the scans were inspected using the agilent quality control report as well as an internal quality control metric including spike-in rna controls.
experimental design and statistical analysis of microarray data
the agilent  <dig> ×  <dig> k whole genome microarray features four sub-arrays per slide.
a reference design  <cit>  with sub-arrays as blocks of size two  was used to analyze the median signal intensities of the two-color microarray data.
the experiment included the main effects of treatment groups  plus the sub-array as a block term.
five biological replicates per condition were used for three of the four groups, while four replicates were used for the remaining hypo+ condition , yielding a total of  <dig> microarrays.
an identical experimental structure was followed for each gender.
separate statistical analyses were carried out for the male and female samples, respectively, to eliminate any correlation effect due to the dam as samples within genders were independent.
all pre-processing of the data was conducted using r. the data were normalized using loess normalization  <cit>  in the r library "maanova".
background fluorescence was measured using the 3xslv <dig> negative control probes; probes with median signal intensities less than the trimmed mean  plus three trimmed standard deviations of the 3xslv <dig> probes were flagged as absent .
ratio intensity plots  were constructed for the raw and normalized data for each array.
boxplots and hierarchical clustering using average linkage were generated to identify outlier arrays  <cit> .
two outlier microarrays were removed from each data set  based on separate hierarchical clustering of each dataset into sample = cy <dig> and reference = cy <dig> groups.
as a result, the final analysis was based on  <dig> microarrays for females, and  <dig> for males.
genes that were up- and down-regulated in any of the three treatment groups  relative to control were identified using the r library "maanova"  <cit> .
the required anova model was fitted to include the main effects of dose plus the sub-array as a block term.
the fs statistic  <cit> , a shrinkage estimator, was used for the gene-specific variance components, and the associated p-values for all the statistical tests were estimated using the permutation method .
these p-values were then adjusted for multiple comparisons using the false discovery rate  approach  <cit> .
all data are available through the gene expression omnibus  website, accession number: gse <dig>
the least squares mean  <cit> , a function of the model parameters, was used to estimate the fold-change for each pairwise comparison of interest
bioinformatics
promoter regions ) of genes whose expression profile suggested direct regulation by ths were downloaded from the ucsc genome browser .
using a list of validated mouse tres gathered from the literature, a position weight matrix  was developed to score the information content  for position one to six of tre half-sites.
in this way, each tre half-site was assigned a 'score' against the pwm.
the score for each tre half-site was obtained iteratively by leaving the subject out of the pwm construction.
previous bootstrapping analysis of the validated tres in mice revealed a difference  between the distribution of scores for the two half-site scores.
cross-validation of the tres from the literature was carried out with each half-site analyzed separately.
this analysis revealed that a low threshold score  for one half-site and a high threshold score  for the second half-site  allowed minimization of type i and type ii errors .
we therefore used the cut-off scores identified in this cross-validation exercise and scanned the promoter regions for putative tres.
the scan searched for three different types of tres: direct repeats with a four nucleotide spacer ; inverted repeats with no spacer ; and everted repeats with a spacer of six nucleotides.
.
orthologs in humans and in rats were matched using biomart http://www.biomart.org for the genes that were identified as having potential tres in mice.
promoter regions  for these orthologous genes were obtained.
a local alignment using a sliding window principle was used to align the tres identified in mice to the promoter regions of the corresponding orthologous genes in rats and humans.
a modified substitution matrix favouring conservation of the guanine in position two and three of the tre hexamer half-sites was used, since these two sites appear to be highly retained.
three additional features were used to filter identified tres: i) a maximum mismatch of two when comparing mouse to rat, and mouse to human, ii) a maximum distance of  <dig> kb between tres when comparing mouse to rat and mouse to human, and iii) a minimum score determined by the pwm scan.
gene ontology analysis and principal component analysis
genes with fdr corrected p-values smaller than  <dig>  in at least one experimental condition were used for principal component analysis  and pathway/ontology analysis.
pca was carried out using mean centering and scaling on data from genes that were significant in at least one treatment group using genespring gx  <dig> .
<dig> .
mean fold-change, averaged across all replicates per sex per treatment, was subjected to functional enrichment and gene ontology  analysis using ingenuity pathway analysis .
real-time quantitative pcr analysis
total rna was reverse transcribed into cdna using superscriptiii  as per the manufacturer's instructions.
iq sybr green supermix  was used with a reaction volume of  <dig> μl.
all primers were designed using beacon designer  <dig> .
hprt was used as the reference gene.
microarray analysis and rt-qpcr analysis confirmed that hprt expression was stable across all groups.
reactions were carried out in 96-well plates.
pcr were done in duplicate and each plate contained all samples for the gene of interest and reference gene for one sex.
analysis was carried out using a cfx <dig> real-time pcr detection systems .
melting curves were performed for each reaction to ensure primer specificity.
raw data were up-loaded into relative expression software tool  <dig> .
where fold-changes and p-values were calculated.
chip-pcr
livers from the male euthyroid control group were used for the chip using the ez chip kit , according to the manufacturer's instructions.
briefly, a small piece of liver was homogenized with a hand-held homogenizer in  <dig> μl pbs containing broad-spectrum protease inhibitors, and was then cross-linked with 1% formaldehyde.
cross-linking was stopped with glycine and nuclei were collected by adding lysis buffer.
to ensure that dna fragments ranged from  <dig> to  <dig> bp, the nuclear solution was sonicated in an ice bath with  <dig> second bursts at 28% amplitude.
fifteen bursts were completed, each separated by a  <dig> second period.
dna fragment size was verified by agarose gel electrophoresis.
six percent  of the sonicated solution was stored at -20°c as total input , while the remainder was incubated with anti-trβ polyclonal antibody  overnight with agitation at 4°c.
antibody-bound chromatin was precipitated with protein g conjugated agarose beads, washed with gradient stringent buffers, and eluted with elution buffer as per the manufacturer's instructions.
both the eluted solution and the stored ti were incubated at 65°c overnight to reverse cross-links.
immunoprecipitated  dna and ti dna were purified by treatment with rnase, proteinase k and multiple phenol:chloroform:isoamyl alcohol  extractions.
equivalent amounts of ip dna and ti dna were amplified in parallel, using a random primer method with genomeplex complete whole genome amplification kit .
amplified dna was then purified using genelute pcr clean-up kit .
primers targeting the tre identified by bioinformatics analysis were designed using beacon designer  <dig>  pcrs were performed using expand high fidelity pcr system .
β-actin  was used as a negative control and mlxipl  was used as a positive control  <cit> .
gel analyses and band quantifications were carried out using genetools .
binding ratios were calculated by dividing the ratio of ip/ti from the gene of interest by the ratio of ip/ti from the reference gene β-actin.
mean enrichment, standard deviations, and p-values were calculated using log ratios and were then back-transformed.
electrophoretic mobility shift assays
probes were designed against the specific genomic regions of the newly identified tres.
probe size ranged between  <dig> and  <dig> bp.
the probe contained the putative tre sequence plus six to seven nucleotides on each side.
oligonucleotides were labelled using biotin 3' end dna labeling kit  as per the manufacturer's instructions and were then annealed.
binding reaction and detection of complexed biotin-labeled dna-protein was accomplished using the lightshift chemiluminescent emsa kit  as per the manufacturer's instructions.
binding reactions were carried out in  <dig> μl volumes and had final concentrations of  <dig> % glycerol,  <dig> mn mgcl <dig>   <dig> ng/μl poly  and  <dig> % np- <dig>  one μg of chicken trα  was used per reaction.
an amino acid alignment of the trα dna-binding domains from mus musculus and gallus gallus showed an approximate 98% sequence homology.
amino acids involved in direct base contacts  <cit>  were specifically examined and were not found to be substituted between species.
we have previously found that trα and trβ interact with tres in a very similar manner  <cit> .
amino acid alignment of dna-binding domains showed an approximate 86% sequence homology between the tr isoforms.
amino acids involved in direct base contacts  <cit>  were specifically examined and showed 100% conservation between the two isoforms.
biotin end-labelled target dna was diluted 1: <dig> and  <dig> μl were used per reaction.
unlabeled probes were used at a 200-fold molar excess compared to labelled probe.
for the supershift,  <dig> μg of antibody trα/β  or normal mouse igg  were used per reaction.
samples were run on a 5% polyacrylamide gel, transferred onto biodyne precut a nylon membranes  by electrophoretic transfer and uv-cross-linked.
membranes were scanned using chemidoc xrs+  and images analysed using image lab version  <dig> .
<dig>
RESULTS
validation of the mmi/sodium perchlorate mouse model
serum t <dig> levels in all pnd  <dig> pups were measured .
a significant increase was observed in the hyper pup group compared to control levels, whereas a significant decrease was observed in the hypo pup group compared to control levels.
the hypo+ and hypo++ pup groups also showed a significant increase in serum t <dig> compared to control.
the hypo++ pup group received the same doses of ths in the i.p.
injection as the hypo+, but registered higher levels of circulating t <dig> .
levels of t <dig> in hyperthyroid and hypo++ groups were not significantly different.
part of this data has already been reported  <cit> .
genes significantly altered by th perturbations
maanova analysis identified approximately  <dig> significantly altered genes in male and/or female pups, with a fdr-adjusted p <  <dig>  in at least one of the treatment conditions .
fold-changes ranged from  <dig>  to  <dig> , and averaged approximately  <dig> .
in hypo mice,  <dig> genes were significantly altered relative to euthyroid controls.
of these, transcription levels of  <dig> genes were reduced and  <dig> were increased in hypo livers relative to controls.
in the hyperthyroid group,  <dig> genes were significantly altered.
of these, transcript levels of  <dig> genes were reduced and  <dig> were increased relative to controls.
in the hypo+ group,  <dig> genes were significantly altered with  <dig> reduced and  <dig> increased.
females exhibited more significantly altered genes than males.
females had  <dig> significantly altered genes whereas males had  <dig>  a comparison of these two lists revealed that there were  <dig> genes that were differentially expressed by at least one treatment condition in both females and males .
using the full set of significant genes  for all groups , a principal component analysis  was conducted .
this analysis revealed that expression patterns were highly correlated within sexes and treatment groups.
the pca revealed distinct male and female clusters  and four clusters were also found to separate control, hyper, hypo and hypo+ groups.
thus, the expression profiles reveal both an overall treatment as well as sex effect, with both sexes showing similarities in their response within treatment groups.
the list of all significant genes was subjected to go and pathway analysis.
the analysis revealed several pathways that were significantly altered by the treatment.
the five most affected pathways are presented in additional file  <dig>  along with the names of the genes that were affected by th treatment within the pathway.
many genes were found to be involved in oxidative stress response and xenobiotic metabolism/signalling.
a group of genes involved in tr/rxr activation was also identified.
by comparing the response profile of genes that were significantly altered in the hyper, hypo and hypo+ groups, we identified a subset of genes that are most likely to be directly regulated by th rather than by downstream effects.
for example, a directly regulated gene would be expected to show similar responses in the hyper and hypo+ groups, and should be in the opposite direction to the hypo group.
the hypo+ group is particularly informative, as these animals were sacrificed four hours after a single injection of th, but following a more prolonged state of hypothyroidism.
in addition, genes showing similar responses in hypo and hypo+ groups would be less likely to be directly regulated by th.
genes showing a pattern suggesting direct regulation by th are listed in table  <dig>
fold-changes for males and females are presented in comparison to control animals.
all fold-changes have an adjusted p ≤  <dig> .
asterisks  denote genes that were analysed by rt-qpcr.
see additional file  <dig> and  <dig> for p-value data.
microarray data validation and hypothyroid replacement analysis
microarray data were validated using real time rt-qpcr for subsets of significantly altered genes .
first, we investigated the responses of three genes that have previously been shown to be directly regulated by th .
some of these genes were only significant with rt-qpcr .
thus, the data produced by the microarray analysis are likely to be conservative, and include some false-negatives.
fold-changes are presented for treated animals relative to control animals; n =  <dig> for all groups except for hyper, where n =  <dig>  fold-changes  and p-values are in bold when p ≤  <dig> .
second, we validated several genes suspected of being directly regulated by th using rt-qpcr.
for this analysis, we compared gene expression changes in the four treatment groups, as well as in a group of animals subjected to the hypo+ treatment but whose circulating t <dig> levels were substantially higher than those in the hypo+ animals used for microarray analysis .
here we present the results from the hypo+ and hypo++ groups when compared to control data .
the direction of fold-change was consistent with the microarray data for each of the nine genes selected for rt-qpcr validation.
a significant fold-change in the hypo++ group was observed in all genes evaluated except for tor1a and trim <dig>
the hypo++ group was not included in the microarray analysis.
the fold-changes represent treated relative to control animals; for rt-qpcr data n =  <dig> for each group.
fold-changes  and p-values are in bold when p ≤  <dig> .
asterisks  denote that the data were taken from female pups due to the absence of significant data from male pups.
h01rik corresponds to 2310003h01rik.
absence of data for rt-qpcr values were because of undetectable signal due to low expression level.
bioinformatics analysis of direct th-regulated genes
a search for tres in the  <dig> regulatory regions of genes suspected to be under direct th-regulation  identified  <dig> tres in  <dig> different genes.
when locally aligned with the rat promoter region of the corresponding orthologs, the list was reduced to  <dig> tres in  <dig> different genes.
once this was accomplished, the remaining tres were locally aligned with the human promoter region of the corresponding orthologs to produce the final list of  <dig> tres found in  <dig> different genes .
the candidate tre sequences in rats or humans were within  <dig> kb of the relative distance of the putative tre identified in mice, and had two or less mismatches when compared to the mouse tre.
thus, this final list of  <dig> tres provides the most likely candidate tre-containing promoters that may be conserved between mice, rats and humans.
tres were found to be conserved between mice, rats and humans  when compared to the mouse tre).
mms were only considered in the bolded half-sites.
positions  of the tres are presented in reference to the transcription start site of the associated gene.
two previously characterized tres are presented at the bottom of the table.
* tres validated by chip-pcr.
‡ corresponds to 2310003h01rik.
a characterized by  <cit> .
b characterized by  <cit> ; bolding of half-sites and positions of tre were corrected.
chip-pcr analysis
chip-pcr, using euthyroid livers, was performed to validate tr-tre binding of six tres selected from the list of genes thought to contain conserved tres .
trβ- <dig> antibody was used to precipitate the protein-dna complexes followed by pcr analysis to compare ti  to ip samples.
antibody specificity has previously been demonstrated  <cit> .
mlxipl was used as a positive control as it binds tr in its promoter region, whereas β-actin was used as a negative control.
mean enrichment of mlxipl compared to β-actin is shown as well as the individual enrichment of the three biological samples .
tres identified in the promoter regions of tor1a, 2310003h01rik, hect3d and slc25a <dig> were enriched by 4-fold or more in the ip compared to ti .
candidate tres analysed in the promoter regions of ihh and arl4d showed no apparent enrichment in the ip samples .
electrophoretic mobility shift assays
emsas were used to demonstrate the ability of the suspected tres to bind to the tr.
a dr <dig> containing two perfect half-sites  was used as a positive control and showed a shift in electrophoretic mobility and a decrease in detection  when unlabelled probe was added .
the dr <dig> positive control tested with an antibody against trα/β exhibited a decrease in signal  of the shifted bands, and a slightly detectable supershifted band.
in contrast, the non-specific mouse igg antibody showed no reduction in detection of the shifted bands.
probes with sequences of suspected tres associated with the genes slc25a <dig>  hectd <dig> and tor1a all showed a shift in electrophoretic mobility, as well as a decrease in detection when unlabelled probe  or unlabelled dr <dig> positive control probe  was added .
thus, the emsa results help to support the idea that the tr can bind to putative tres in the promoter regions of slc25a <dig>  hectd <dig> and tor1a.
discussion
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  <dig> livers were studied as animals of this age have circulating t <dig> levels higher than at any other age  <cit>  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  <cit> .
several hundred genes were responsive to this short-term treatment, with fold-changes as high as  <dig> .
however, the vast majority of these genes had relatively small fold-changes .
when filtered on a  <dig>  fold-change cutoff, only  <dig> genes in females and  <dig> genes in males were perturbed, with an overlap of  <dig> 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 .
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  <cit>  and by our group  <cit> .
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 .
in addition, several specific genes previously shown to be involved in th hepatic gene regulation, such as me <dig>  thrsp and dio <dig> 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  <dig> 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 .
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 t <dig> 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  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 me <dig>  thrsp or dio <dig>  were confirmed to be differentially regulated in our experiment, they did not meet the stringent requirements set for this particular evaluation .
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 ≤  <dig>  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 .
we retrieved the relevant dna sequences and identified potential tres within the promoter regions of the responsive genes .
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  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  can be found  <dig> kb away from of the tss  <cit> , whereas others can be more than  <dig> kb away from the tss  <cit> .
in addition, interchromosomal regulatory interactions have also been documented  <cit> .
our gene expression profiling analysis identified  <dig> genes that appeared to be under direct th-regulation.
within these  <dig> genes, bioinformatics mining for three different types of tres  revealed  <dig> 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 - <dig> kb to + <dig> kb promoter regions of these genes we were able to identified  <dig> tres in  <dig> 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  <dig>  the tre near klf <dig> has been characterized in mice, rats, and humans  <cit> .
comparison of the tres within the orthologous klf <dig> gene of the three species shows that the tres are within  <dig> kb of each other relative to the tss.
the tre located in the mbp gene promoter region has been characterized in mice and rats  <cit> .
the tres for mbp for both species are within  <dig> 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  <dig> different pairs of orthologs between mice and humans  <cit> .
eres identified in proximity of these orthologs were within  <dig> 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 ir <dig>  where the two half-sites are arranged in a palindromic organization separated by a  <dig> bp spacer.
the classic tre, on the other hand, is a dr <dig> with the two tandem repeated half-sites separated by a  <dig> bp spacer  <cit> , although ir <dig>  er <dig> 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  <cit> .
the authors looked at the half-site recognition for retinoic acid receptor  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  <cit>  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 dr <dig> 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  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 .
functional tres have been described with a high level of degeneration, multiple spacer sizes and various half-site organizations  <cit>  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  <cit> , 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  <cit> .
micrornas  could also potentially regulate th-mediated mrna expression.
we recently identified significant alterations in  <dig> different mirnas in the livers of pnd  <dig> hypothyroid mice  <cit> .
lastly, within our list of significantly altered genes , several transcription factors are present, including rarβ, esr <dig> and nr3c <dig> .
thus, genes present within our list may be under the control of other transcription factors that are directly regulated by ths.
for example, esr <dig> expression is decreased in hypothyroid female mice.
this could lead to many downstream transcriptional effects since esr <dig> is a transcription factor that has been shown to interact with  <dig> different nuclear receptors .
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  <dig> 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 slc25a <dig> correspond to a dr <dig> organization.
the tre identified in the promoter region of hectd <dig> corresponds to an er <dig> organization, whereas the tre identified in the promoter region of 2310003h01rik corresponds to an ir <dig> organization.
the genomic locations of these newly identified tres all fall in upstream promoter regions, - <dig> , - <dig> , - <dig>  and - <dig>  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  <dig> to - <dig>  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 klf <dig>  and is - <dig>  kb from its tss.
this information is based on the  <dig> mouse tres that we were able to find in the literature .
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  <dig> or  <dig> 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 slc25a <dig>  hectd <dig> 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 dr <dig> unlabelled probe was added.
the addition of unlabelled probe demonstrates the specificity of the shifted band, whereas the addition of the unlabelled dr <dig> 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 dr <dig> positive control probe.
a correlation between the chip-pcr and emsa results was observed.
scl25a <dig> has the strongest enrichment measured by chip-pcr, followed by hectd <dig> and tor1a.
slc25a <dig> again showed the strongest decrease in signal intensity when unlabelled specific probe was added in the emsa, followed by hectd <dig> and tor1a.
slc25a <dig> is involved in the transportation of molecules across the mitochondrial membrane  <cit> .
ths have been shown to induce mitochondrial biogenesis and enhance atp production  <cit> .
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.
hectd <dig> 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  <cit> .
ths have also been shown to play a role in cytoskeletal protein regulation  <cit> .
thus, a tre in the promoter region of hectd <dig> could link this gene to cytoskeletal regulation via ths.
in mice 2310003h01rik is the predicted ortholog of faap <dig>  fanconi anemia-associated protein,  <dig> kda  <cit> .
faap <dig> is involved in the fanconi anemia  core complex, which plays a role in the dna damage response network  <cit> .
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  <cit> .
children with fa have a high risk of endocrine abnormalities including hypothyroidism  <cit> .
a recent study found that th therapy could improve the linear growth of children suffering from fa  <cit> .
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  <cit> .
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  <cit> .
thus, the presence of a tre in the promoter region of tor1a could suggest a potential role for th-tor1a interaction in muscle development.
CONCLUSIONS
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  <dig> 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, slc25a <dig>  hectd <dig> 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.
authors' contributions
conceived and designed the experiments: mp, hd, mw and cy.
performed the experiments: mp, hd and rg.
analyzed the data: mp, rg, aw and mm.
mp drafted the manuscript.
all authors revised and approved the final manuscript.
supplementary material
additional file 1
altered genes in male pups, with a fdr-adjusted p <  <dig>  in at least one of the treatment conditions.
click here for file
additional file 2
altered genes in female pups, with a fdr-adjusted p <  <dig>  in at least one of the treatment conditions.
click here for file
additional file 3
the five most affected pathways, using the list of all significant genes , along with the names of the genes that were affected by th treatment within the pathway.
click here for file
acknowledgements
this work was funded by the health canada genomics r&d initiative and the canadian regulatory system for biotechnology.
we would like to gratefully acknowledge the expert advice provided by dr. r.t. zoeller in the development of this animal model, and for very helpful discussions.
we also thank guillaume pelletier and lynn berndt-weis for helpful comments on the manuscript.

 the thyroid participates in the regulation of basic physiological processes by producing thyroid hormones, which include thyroxine and triiodothyronine. ths exert their effects on growth, development and metabolism of practically every cell and organ. their primary effect is the transcriptional regulation of target genes. this occurs when ths interact with th receptors. 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. 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. trs can bind to tres as monomers or homodimers, although it is thought that trs interact as heterodimers with the retinoid × receptor the majority of the time. 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. less severe and/or transient hypothyroidism has also been shown to cause adverse structural and functional effects. in humans, subclinical hypothyroidism during fetal development has been shown to result in reduced cognitive function. exposures to a broad array of substances, including both natural and synthetic chemicals, have been shown to alter th physiology. environmental contaminants such as perchlorates, polychlorinated biphenyls, bisphenol a, polybrominated diphenyl ethers, triclosan and nitrate 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. 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. 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 ). 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 15) rendered transiently hypo- or hyperthyroid. this developmental period corresponds to a dramatic increase in circulating ths. in pnd mice, circulating t is higher than at any other age 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. methods animals and exposures all animal care and handling was in accordance with canadian council for animal care guidelines and was reviewed by the health canada animal care committee prior to commencement of the study. c57bl/ mice were purchased from charles river and were housed in hanging polycarbonate cages under a 12: hr light-dark cycle at 23°c. animals were provided with shelters, nesting material, food ad libitum and drinking water containing 1% sucrose. mice were acclimated for days. breeding was accomplished by transferring two sexually mature female mice into the home cage of a sexually mature male. after four nights of co-housing, each female was transferred to a separate cage. dams were allowed to litter naturally and pup numbers were not adjusted. to produce hypothyroid pups, dams were provided with drinking water that contained methimizole and perchlorate for three days. to produce hyperthyroid pups, intraperitoneal injections of ths were administered to pups on pnd four hours before decapitation and tissue collection. for the hypothyroid/replacement group; dams were provided with drinking water that contained mmi and perchlorate for three days. pups then received i.p. injections of ths on pnd four hours before decapitation and tissue collection. the control and hypo groups received saline i.p. injections four hours before decapitation and tissue collection. tissue collection, rna extraction and purification on pnd pups were sacrificed by decapitation. serum prepared from trunk blood using serum separator tubes was used for t analyses by radioimmunoassay kits as per the manufacturer's instructions. the hypo + group was subdivided into two sub groups - hereafter called hypo+ and hypo++ - based on measured serum t levels. animals allotted to the hypo++ group had serum t levels that approximated those seen in the hyperthyroid group while those that remained as hypo+ had t levels near or slightly in excess of control levels. pup livers were rapidly dissected and flash frozen in liquid nitrogen. total rna was extracted from liver samples with trizol reagent followed by rneasy mini kit clean-up according to the manufacturer's instructions. rna quality assessment was determined by nanodrop and agilent bioanalyzer and rna nanolab chip kit. all samples had 260/ ratios over, and rna integrity numbers over. microarray hybridization hepatic gene expression was assessed by hybridizing samples to agilent × k whole genome microarrays using a reference design. briefly, ng of total rna from the liver was used to synthesize double-stranded cdna and cyanine labelled crna. experimental samples were labelled with cyanine 5-ctp and reference rna with cyanine 3-ctp according to the manufacturer's instructions. cyanine labelled crna targets were in vitro transcribed using t rna polymerase and purified using rneasy mini kit. experimental and reference samples were hybridized to arrays at 60°c for hours. slides were washed and scanned on an agilent microarray scanner. the data were acquired using agilent feature extraction software version.. prior to statistical analysis, the scans were inspected using the agilent quality control report as well as an internal quality control metric including spike-in rna controls. experimental design and statistical analysis of microarray data the agilent × k whole genome microarray features four sub-arrays per slide. a reference design with sub-arrays as blocks of size two was used to analyze the median signal intensities of the two-color microarray data. the experiment included the main effects of treatment groups plus the sub-array as a block term. five biological replicates per condition were used for three of the four groups, while four replicates were used for the remaining hypo+ condition, yielding a total of microarrays. an identical experimental structure was followed for each gender. separate statistical analyses were carried out for the male and female samples, respectively, to eliminate any correlation effect due to the dam as samples within genders were independent. all pre-processing of the data was conducted using r. the data were normalized using loess normalization in the r library "maanova". background fluorescence was measured using the 3xslv negative control probes; probes with median signal intensities less than the trimmed mean plus three trimmed standard deviations of the 3xslv probes were flagged as absent. ratio intensity plots were constructed for the raw and normalized data for each array. boxplots and hierarchical clustering using average linkage were generated to identify outlier arrays. two outlier microarrays were removed from each data set based on separate hierarchical clustering of each dataset into sample = cy and reference = cy groups. as a result, the final analysis was based on microarrays for females, and for males. genes that were up- and down-regulated in any of the three treatment groups relative to control were identified using the r library "maanova". the required anova model was fitted to include the main effects of dose plus the sub-array as a block term. the fs statistic, a shrinkage estimator, was used for the gene-specific variance components, and the associated p-values for all the statistical tests were estimated using the permutation method. these p-values were then adjusted for multiple comparisons using the false discovery rate approach. all data are available through the gene expression omnibus website, accession number: gse the least squares mean, a function of the model parameters, was used to estimate the fold-change for each pairwise comparison of interest bioinformatics promoter regions ) of genes whose expression profile suggested direct regulation by ths were downloaded from the ucsc genome browser. using a list of validated mouse tres gathered from the literature, a position weight matrix was developed to score the information content for position one to six of tre half-sites. in this way, each tre half-site was assigned a'score' against the pwm. the score for each tre half-site was obtained iteratively by leaving the subject out of the pwm construction. previous bootstrapping analysis of the validated tres in mice revealed a difference between the distribution of scores for the two half-site scores. cross-validation of the tres from the literature was carried out with each half-site analyzed separately. this analysis revealed that a low threshold score for one half-site and a high threshold score for the second half-site allowed minimization of type i and type ii errors. we therefore used the cut-off scores identified in this cross-validation exercise and scanned the promoter regions for putative tres. the scan searched for three different types of tres: direct repeats with a four nucleotide spacer; inverted repeats with no spacer; and everted repeats with a spacer of six nucleotides.. orthologs in humans and in rats were matched using biomart http://www.biomart.org for the genes that were identified as having potential tres in mice. promoter regions for these orthologous genes were obtained. a local alignment using a sliding window principle was used to align the tres identified in mice to the promoter regions of the corresponding orthologous genes in rats and humans. a modified substitution matrix favouring conservation of the guanine in position two and three of the tre hexamer half-sites was used, since these two sites appear to be highly retained. three additional features were used to filter identified tres: i) a maximum mismatch of two when comparing mouse to rat, and mouse to human, ii) a maximum distance of kb between tres when comparing mouse to rat and mouse to human, and iii) a minimum score determined by the pwm scan. gene ontology analysis and principal component analysis genes with fdr corrected p-values smaller than in at least one experimental condition were used for principal component analysis and pathway/ontology analysis. pca was carried out using mean centering and scaling on data from genes that were significant in at least one treatment group using genespring gx.. mean fold-change, averaged across all replicates per sex per treatment, was subjected to functional enrichment and gene ontology analysis using ingenuity pathway analysis. real-time quantitative pcr analysis total rna was reverse transcribed into cdna using superscriptiii as per the manufacturer's instructions. iq sybr green supermix was used with a reaction volume of μl. all primers were designed using beacon designer. hprt was used as the reference gene. microarray analysis and rt-qpcr analysis confirmed that hprt expression was stable across all groups. reactions were carried out in 96-well plates. pcr were done in duplicate and each plate contained all samples for the gene of interest and reference gene for one sex. analysis was carried out using a cfx real-time pcr detection systems. melting curves were performed for each reaction to ensure primer specificity. raw data were up-loaded into relative expression software tool. where fold-changes and p-values were calculated. chip-pcr livers from the male euthyroid control group were used for the chip using the ez chip kit, according to the manufacturer's instructions. briefly, a small piece of liver was homogenized with a hand-held homogenizer in μl pbs containing broad-spectrum protease inhibitors, and was then cross-linked with 1% formaldehyde. cross-linking was stopped with glycine and nuclei were collected by adding lysis buffer. to ensure that dna fragments ranged from to bp, the nuclear solution was sonicated in an ice bath with second bursts at 28% amplitude. fifteen bursts were completed, each separated by a second period. dna fragment size was verified by agarose gel electrophoresis. six percent of the sonicated solution was stored at -20°c as total input, while the remainder was incubated with anti-trβ polyclonal antibody overnight with agitation at 4°c. antibody-bound chromatin was precipitated with protein g conjugated agarose beads, washed with gradient stringent buffers, and eluted with elution buffer as per the manufacturer's instructions. both the eluted solution and the stored ti were incubated at 65°c overnight to reverse cross-links. immunoprecipitated dna and ti dna were purified by treatment with rnase, proteinase k and multiple phenol:chloroform:isoamyl alcohol extractions. equivalent amounts of ip dna and ti dna were amplified in parallel, using a random primer method with genomeplex complete whole genome amplification kit. amplified dna was then purified using genelute pcr clean-up kit. primers targeting the tre identified by bioinformatics analysis were designed using beacon designer pcrs were performed using expand high fidelity pcr system. β-actin was used as a negative control and mlxipl was used as a positive control. gel analyses and band quantifications were carried out using genetools. binding ratios were calculated by dividing the ratio of ip/ti from the gene of interest by the ratio of ip/ti from the reference gene β-actin. mean enrichment, standard deviations, and p-values were calculated using log ratios and were then back-transformed. electrophoretic mobility shift assays probes were designed against the specific genomic regions of the newly identified tres. probe size ranged between and bp. the probe contained the putative tre sequence plus six to seven nucleotides on each side. oligonucleotides were labelled using biotin 3' end dna labeling kit as per the manufacturer's instructions and were then annealed. binding reaction and detection of complexed biotin-labeled dna-protein was accomplished using the lightshift chemiluminescent emsa kit as per the manufacturer's instructions. binding reactions were carried out in μl volumes and had final concentrations of % glycerol, mn mgcl ng/μl poly and % np- one μg of chicken trα was used per reaction. an amino acid alignment of the trα dna-binding domains from mus musculus and gallus gallus showed an approximate 98% sequence homology. amino acids involved in direct base contacts were specifically examined and were not found to be substituted between species. we have previously found that trα and trβ interact with tres in a very similar manner. amino acid alignment of dna-binding domains showed an approximate 86% sequence homology between the tr isoforms. amino acids involved in direct base contacts were specifically examined and showed 100% conservation between the two isoforms. biotin end-labelled target dna was diluted 1: and μl were used per reaction. unlabeled probes were used at a 200-fold molar excess compared to labelled probe. for the supershift, μg of antibody trα/β or normal mouse igg were used per reaction. samples were run on a 5% polyacrylamide gel, transferred onto biodyne precut a nylon membranes by electrophoretic transfer and uv-cross-linked. membranes were scanned using chemidoc xrs+ and images analysed using image lab version.  validation of the mmi/sodium perchlorate mouse model serum t levels in all pnd pups were measured. a significant increase was observed in the hyper pup group compared to control levels, whereas a significant decrease was observed in the hypo pup group compared to control levels. the hypo+ and hypo++ pup groups also showed a significant increase in serum t compared to control. the hypo++ pup group received the same doses of ths in the i.p. injection as the hypo+, but registered higher levels of circulating t. levels of t in hyperthyroid and hypo++ groups were not significantly different. part of this data has already been reported. genes significantly altered by th perturbations maanova analysis identified approximately significantly altered genes in male and/or female pups, with a fdr-adjusted p < in at least one of the treatment conditions. fold-changes ranged from to, and averaged approximately. in hypo mice, genes were significantly altered relative to euthyroid controls. of these, transcription levels of genes were reduced and were increased in hypo livers relative to controls. in the hyperthyroid group, genes were significantly altered. of these, transcript levels of genes were reduced and were increased relative to controls. in the hypo+ group, genes were significantly altered with reduced and increased. females exhibited more significantly altered genes than males. females had significantly altered genes whereas males had a comparison of these two lists revealed that there were genes that were differentially expressed by at least one treatment condition in both females and males. using the full set of significant genes for all groups, a principal component analysis was conducted. this analysis revealed that expression patterns were highly correlated within sexes and treatment groups. the pca revealed distinct male and female clusters and four clusters were also found to separate control, hyper, hypo and hypo+ groups. thus, the expression profiles reveal both an overall treatment as well as sex effect, with both sexes showing similarities in their response within treatment groups. the list of all significant genes was subjected to go and pathway analysis. the analysis revealed several pathways that were significantly altered by the treatment. the five most affected pathways are presented in additional file along with the names of the genes that were affected by th treatment within the pathway. many genes were found to be involved in oxidative stress response and xenobiotic metabolism/signalling. a group of genes involved in tr/rxr activation was also identified. by comparing the response profile of genes that were significantly altered in the hyper, hypo and hypo+ groups, we identified a subset of genes that are most likely to be directly regulated by th rather than by downstream effects. for example, a directly regulated gene would be expected to show similar responses in the hyper and hypo+ groups, and should be in the opposite direction to the hypo group. the hypo+ group is particularly informative, as these animals were sacrificed four hours after a single injection of th, but following a more prolonged state of hypothyroidism. in addition, genes showing similar responses in hypo and hypo+ groups would be less likely to be directly regulated by th. genes showing a pattern suggesting direct regulation by th are listed in table fold-changes for males and females are presented in comparison to control animals. all fold-changes have an adjusted p ≤. asterisks denote genes that were analysed by rt-qpcr. see additional file and for p-value data. microarray data validation and hypothyroid replacement analysis microarray data were validated using real time rt-qpcr for subsets of significantly altered genes. first, we investigated the responses of three genes that have previously been shown to be directly regulated by th. some of these genes were only significant with rt-qpcr. thus, the data produced by the microarray analysis are likely to be conservative, and include some false-negatives. fold-changes are presented for treated animals relative to control animals; n = for all groups except for hyper, where n = fold-changes and p-values are in bold when p ≤. second, we validated several genes suspected of being directly regulated by th using rt-qpcr. for this analysis, we compared gene expression changes in the four treatment groups, as well as in a group of animals subjected to the hypo+ treatment but whose circulating t levels were substantially higher than those in the hypo+ animals used for microarray analysis. here we present the results from the hypo+ and hypo++ groups when compared to control data. the direction of fold-change was consistent with the microarray data for each of the nine genes selected for rt-qpcr validation. a significant fold-change in the hypo++ group was observed in all genes evaluated except for tor1a and trim the hypo++ group was not included in the microarray analysis. the fold-changes represent treated relative to control animals; for rt-qpcr data n = for each group. fold-changes and p-values are in bold when p ≤. asterisks denote that the data were taken from female pups due to the absence of significant data from male pups. h01rik corresponds to 2310003h01rik. absence of data for rt-qpcr values were because of undetectable signal due to low expression level. bioinformatics analysis of direct th-regulated genes a search for tres in the regulatory regions of genes suspected to be under direct th-regulation identified tres in different genes. when locally aligned with the rat promoter region of the corresponding orthologs, the list was reduced to tres in different genes. once this was accomplished, the remaining tres were locally aligned with the human promoter region of the corresponding orthologs to produce the final list of tres found in different genes. the candidate tre sequences in rats or humans were within kb of the relative distance of the putative tre identified in mice, and had two or less mismatches when compared to the mouse tre. thus, this final list of tres provides the most likely candidate tre-containing promoters that may be conserved between mice, rats and humans. tres were found to be conserved between mice, rats and humans when compared to the mouse tre). mms were only considered in the bolded half-sites. positions of the tres are presented in reference to the transcription start site of the associated gene. two previously characterized tres are presented at the bottom of the table. * tres validated by chip-pcr. ‡ corresponds to 2310003h01rik. a characterized by. b characterized by; bolding of half-sites and positions of tre were corrected. chip-pcr analysis chip-pcr, using euthyroid livers, was performed to validate tr-tre binding of six tres selected from the list of genes thought to contain conserved tres. trβ- antibody was used to precipitate the protein-dna complexes followed by pcr analysis to compare ti to ip samples. antibody specificity has previously been demonstrated. mlxipl was used as a positive control as it binds tr in its promoter region, whereas β-actin was used as a negative control. mean enrichment of mlxipl compared to β-actin is shown as well as the individual enrichment of the three biological samples. tres identified in the promoter regions of tor1a, 2310003h01rik, hect3d and slc25a were enriched by 4-fold or more in the ip compared to ti. candidate tres analysed in the promoter regions of ihh and arl4d showed no apparent enrichment in the ip samples. electrophoretic mobility shift assays emsas were used to demonstrate the ability of the suspected tres to bind to the tr. a dr containing two perfect half-sites was used as a positive control and showed a shift in electrophoretic mobility and a decrease in detection when unlabelled probe was added. the dr positive control tested with an antibody against trα/β exhibited a decrease in signal of the shifted bands, and a slightly detectable supershifted band. in contrast, the non-specific mouse igg antibody showed no reduction in detection of the shifted bands. probes with sequences of suspected tres associated with the genes slc25a hectd and tor1a all showed a shift in electrophoretic mobility, as well as a decrease in detection when unlabelled probe or unlabelled dr positive control probe was added. thus, the emsa results help to support the idea that the tr can bind to putative tres in the promoter regions of slc25a hectd and tor1a. discussion 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 livers were studied as animals of this age have circulating t levels higher than at any other age 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. several hundred genes were responsive to this short-term treatment, with fold-changes as high as. however, the vast majority of these genes had relatively small fold-changes. when filtered on a fold-change cutoff, only genes in females and genes in males were perturbed, with an overlap of 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. 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 and by our group. 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. in addition, several specific genes previously shown to be involved in th hepatic gene regulation, such as me thrsp and dio 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 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. 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 t 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 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 me thrsp or dio were confirmed to be differentially regulated in our experiment, they did not meet the stringent requirements set for this particular evaluation. 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 ≤ 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. we retrieved the relevant dna sequences and identified potential tres within the promoter regions of the responsive genes. 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 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 can be found kb away from of the tss, whereas others can be more than kb away from the tss. in addition, interchromosomal regulatory interactions have also been documented. our gene expression profiling analysis identified genes that appeared to be under direct th-regulation. within these genes, bioinformatics mining for three different types of tres revealed 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 - kb to + kb promoter regions of these genes we were able to identified tres in 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 the tre near klf has been characterized in mice, rats, and humans. comparison of the tres within the orthologous klf gene of the three species shows that the tres are within kb of each other relative to the tss. the tre located in the mbp gene promoter region has been characterized in mice and rats. the tres for mbp for both species are within 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 different pairs of orthologs between mice and humans. eres identified in proximity of these orthologs were within 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 ir where the two half-sites are arranged in a palindromic organization separated by a bp spacer. the classic tre, on the other hand, is a dr with the two tandem repeated half-sites separated by a bp spacer, although ir er 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. the authors looked at the half-site recognition for retinoic acid receptor 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 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 dr 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 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. functional tres have been described with a high level of degeneration, multiple spacer sizes and various half-site organizations 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, 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. micrornas could also potentially regulate th-mediated mrna expression. we recently identified significant alterations in different mirnas in the livers of pnd hypothyroid mice. lastly, within our list of significantly altered genes, several transcription factors are present, including rarβ, esr and nr3c. thus, genes present within our list may be under the control of other transcription factors that are directly regulated by ths. for example, esr expression is decreased in hypothyroid female mice. this could lead to many downstream transcriptional effects since esr is a transcription factor that has been shown to interact with different nuclear receptors. 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 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 slc25a correspond to a dr organization. the tre identified in the promoter region of hectd corresponds to an er organization, whereas the tre identified in the promoter region of 2310003h01rik corresponds to an ir organization. the genomic locations of these newly identified tres all fall in upstream promoter regions, -, -, - and - 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 to - 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 klf and is - kb from its tss. this information is based on the mouse tres that we were able to find in the literature. 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 or 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 slc25a hectd 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 dr unlabelled probe was added. the addition of unlabelled probe demonstrates the specificity of the shifted band, whereas the addition of the unlabelled dr 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 dr positive control probe. a correlation between the chip-pcr and emsa results was observed. scl25a has the strongest enrichment measured by chip-pcr, followed by hectd and tor1a. slc25a again showed the strongest decrease in signal intensity when unlabelled specific probe was added in the emsa, followed by hectd and tor1a. slc25a is involved in the transportation of molecules across the mitochondrial membrane. ths have been shown to induce mitochondrial biogenesis and enhance atp production. 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. hectd 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. ths have also been shown to play a role in cytoskeletal protein regulation. thus, a tre in the promoter region of hectd could link this gene to cytoskeletal regulation via ths. in mice 2310003h01rik is the predicted ortholog of faap fanconi anemia-associated protein, kda. faap is involved in the fanconi anemia core complex, which plays a role in the dna damage response network. 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. children with fa have a high risk of endocrine abnormalities including hypothyroidism. a recent study found that th therapy could improve the linear growth of children suffering from fa. 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. 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. 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 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, slc25a hectd 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. authors' contributions conceived and designed the experiments: mp, hd, mw and cy. performed the experiments: mp, hd and rg. analyzed the data: mp, rg, aw and mm. mp drafted the manuscript. all authors revised and approved the final manuscript. supplementary material additional file 1 altered genes in male pups, with a fdr-adjusted p < in at least one of the treatment conditions. click here for file additional file 2 altered genes in female pups, with a fdr-adjusted p < in at least one of the treatment conditions. click here for file additional file 3 the five most affected pathways, using the list of all significant genes, along with the names of the genes that were affected by th treatment within the pathway. click here for file acknowledgements this work was funded by the health canada genomics r&d initiative and the canadian regulatory system for biotechnology. we would like to gratefully acknowledge the expert advice provided by dr. r.t. zoeller in the development of this animal model, and for very helpful discussions. we also thank guillaume pelletier and lynn berndt-weis for helpful comments on the manuscript.
a bioinformatics search identified thyroid response elements in the promoter regions of different genes thought to be directly regulated by thyroid hormones. transcriptional profiling of thyroid hormone disrupted animals coupled with a novel bioinformatics search revealed new thyroid response elements associated with genes previously unknown to be responsive to thyroid hormone. thyroid hormones exert their effects through interactions with thyroid receptors that directly bind thyroid response elements and can alter transcriptional activity of target genes. a subgroup of thyroid hormone-disrupted pups received intraperitoneal injections of replacement thyroid hormones four hours prior to sacrifice. thyroid response elements found in the promoter regions of tor1a, 2310003h01rik, hect3d and slc25a were further validated by confirming that the thyroid receptor is associated with these sequences in vivo and that it can bind directly to these sequences in vitro. the effects of short-term thyroid hormone perturbation on hepatic mrna transcription in juvenile mice were evaluated, with the goal of identifying genes containing active thyroid response elements. the regulatory regions of the genome adjacent to these genes were examined for half-site sequences that resemble known thyroid response elements. three different arrangements of thyroid response elements were identified. some of these thyroid response elements were located far up-stream of the transcription start site of the regulated gene. an additional group received only thyroid hormones four hours prior to sacrifice.
 disruption of thyroid hormone signalling can alter growth, development and energy metabolism. thyroid hormones exert their effects through interactions with thyroid receptors that directly bind thyroid response elements and can alter transcriptional activity of target genes. the effects of short-term thyroid hormone perturbation on hepatic mrna transcription in juvenile mice were evaluated, with the goal of identifying genes containing active thyroid response elements. thyroid hormone disruption was induced from postnatal day to by adding goitrogens to dams' drinking water. a subgroup of thyroid hormone-disrupted pups received intraperitoneal injections of replacement thyroid hormones four hours prior to sacrifice. an additional group received only thyroid hormones four hours prior to sacrifice. hepatic mrna was extracted and hybridized to agilent mouse microarrays.  transcriptional profiling enabled the identification of genes that appeared to be under direct thyroid hormone-regulation. the regulatory regions of the genome adjacent to these genes were examined for half-site sequences that resemble known thyroid response elements. a bioinformatics search identified thyroid response elements in the promoter regions of different genes thought to be directly regulated by thyroid hormones. thyroid response elements found in the promoter regions of tor1a, 2310003h01rik, hect3d and slc25a were further validated by confirming that the thyroid receptor is associated with these sequences in vivo and that it can bind directly to these sequences in vitro. three different arrangements of thyroid response elements were identified. some of these thyroid response elements were located far up-stream of the transcription start site of the regulated gene.  transcriptional profiling of thyroid hormone disrupted animals coupled with a novel bioinformatics search revealed new thyroid response elements associated with genes previously unknown to be responsive to thyroid hormone. the work provides insight into thyroid response element sequence motif characteristics. 