Tuning of liver circadian transcriptome rhythms by thyroid hormone state in male mice

To model hypothyroid conditions in mice, we supplemented their drinking water with methimazole and potassium perchlorate (MMI group), which suppresses TH biosynthesis through inhibition of thyroperoxidase in the thyroid gland and inhibits iodine uptake¹². MMI treatment strongly reduced T₃ and T₄ levels and led to strong Tshb mRNA induction in the pituitary (Supplmentary Fig. S1A). T₃ levels showed circadian rhythms in the control (CON) and MMI groups peaking in the late and mid-dark phase, respectively (Fig. 1A). T₄, on the other hand, was arrhythmic in both conditions (Fig. 1B; Supplementary file 1).

Reduced food and water intake were observed in MMI compared to CON mice (Supplementary Fig. 1B,C). As had been reported before¹³—and unlike hypothyroid humans—MMI mice had reduced body weight compared to CON mice (Fig. 1C). Compared to CON mice, locomotor activity rhythms were similar in MMI mice, but activity was reduced specifically during the dark phase (Fig. 1D; Supplement Fig. 1D). Body temperature was reduced in MMI during the dark phase (Fig. 1E; Supplementary Fig. 1E). Energy turnover, assessed by oxygen consumption, was slightly lower throughout the day (Fig. 1F; Supplementary Fig. 1F) while respiratory quotient profiles were largely unaltered in MMI mice (Fig. 1; Supplementary Fig. 1G,H).

By averaging profile data over the whole day and comparing data to high-T₃ conditions¹¹, clear systemic effects of TH state on metabolic homeostasis became apparent. TH state-dependent effects were observed for locomotor activity (r = 0.62, p = 0.0191), food intake (r = 0.85, p = 0.0002), body temperature (r = 0.75, p = 0.0042), and oxygen consumption (r = 0.79, p = 0.032, Fig. 1G–J). In sum, our findings confirm that a low TH state decreases systemic energy turnover, but this effect is more pronounced during the dark (active) phase.

In face of the marked effects of TH state on systemic energy metabolism, we focused our attention on the liver due to its major role as a metabolic organ. We collected tissues every 4 h over the whole day to allow the evaluation of time-of-day dependent effects of TH state on liver transcription. To assess T₃ state effects, we determined differentially expressed genes (DEGs) in T₃ and MMI conditions irrespective of sampling time. Surprisingly, T₃ mice showed a fivefold higher number of DEGs compared to MMI mice (Fig. 2A; Supplementary file 2). DEG analysis was performed separately for each sampling time, and the number of DEGs per timepoint is represented in UpSet plots. DEGs identified for at least one time point yielded 95 DEGs in MMI livers compared to 2200 DEGs (a factor of ca. 23) in high-T₃ mice¹¹ (Fig. 2B,C; Supplementary file 2). We previously identified 37 robust DEGs, i.e., genes consistently differentially expressed across all time points in high TH liver compared to CON¹¹. At a low TH state (MMI), however, no robust DEG was identified (Fig. 2D). Analysis of established liver TH/THR target genes and modulators confirmed a much stronger effect of high-T₃ than low TH conditions on liver TH state (Fig. 2E; Supplementary Fig. 2).

Together, these data show that despite marked systemic effects, MMI treatment had much lower effects on liver transcription. Vice versa, this suggests that, under physiological conditions (CON), the liver is already in a low-TH action.

Evaluation of core clock gene expression profiles showed little effect on rhythmic parameters between MMI and CON mice (Fig. 3A) in line with what was previously observed under high-T₃ conditions¹¹. To fully assess the rhythmic liver transcriptome, the JTK cycle algorithm¹⁴ was used. A total of 3,329 and 3,383 genes (3,354 and 3,397 probe sets) were classified as rhythmic (p < 0.05) in CON and MMI, respectively (Supplementary Fig. 3; Supplementary file 3). 1412 genes (1417 probes) were significantly rhythmic in both groups (Supplementary Fig. 3; Supplementary file 3). An average phase delay of 0.24 h was observed for these genes in MMI livers (Fig. 3B) compared to a 1 h phase advance in high-T₃ conditions¹¹. Expression peak phases of rhythmic genes were widely distributed, and amplitudes were overall similar between conditions (Supplementary Fig. 3).

To investigate the difference in rhythmic characteristics of the transcriptome regulation under low TH conditions, we performed differential circadian rhythm analysis using CircaCompare¹⁵. Of the 5297 genes (5334 probes) showing significant 24-h rhythmicity in at least one group, 1882 genes showed changes in mesor, 403 in amplitude, and 391 in phase between low-TH and control conditions (Fig. 4A,B; Supplementary file 4). Processes associated with cellular proliferation, mRNA processing, and macromolecule complex assembly were enriched in genes with mesor DOWN while processes associated with response to oxidative stress and xenobiotic metabolism processes were found in the mesor UP genes. Genes associated with response to hypoxia, exocytosis, and lipid transport showed increased amplitudes in MMI livers. Conversely, the expression of genes involved in oxidative stress and detoxification processes (cellular oxidant detoxification, hydrogen peroxide catabolism, glutathione, and ethanol metabolism) were phase-delayed in MMI mice (Fig. 4C; Supplementary file 4).

Lipid metabolism-associated processes were regulated for all three rhythm parameters (Fig. 4C). A dual effect (up- and down-regulation) on mesor was observed while most lipid metabolism genes showed reduced amplitudes and phase delays in MMI mice (Fig. 4D). The identified genes were manually inspected and only those genes participating directly in lipid metabolism pathways were included. In this refined gene set for mesor UP, approximately 50% of the genes were associated with acyl-CoA degradation (Acot1, 2, 3, 4, 7, 9, 13) and biosynthesis (Acsm 1, 3, 5) while others were associated with lipolysis (Lpl, Lipa) and TAG uptake (Vldr). In the mesor DOWN group, ca. 40% of the genes were involved in cholesterol biosynthesis (Aacs, Dhcr7, Hmgcr), uptake (Ldlr, Lrp10, Npc1, Lrp5, Pcsk9), and bile acid secretion (Abcb11, Cyp7a1, Abcg5, Slc10a1, Slc10a5, Slc10a2). Moreover, 35% of the genes were associated with FA biosynthesis (Acacb, Fasn), elongation (Elovl1, 2, 3, 5, Hacd1), transport (Cpt1a, Fabp5), and uptake (Slc27a2, Slc27a5). Around 30% of amplitude DOWN genes were associated with cholesterol biosynthesis (Aacs, Hsd17b7) and secretion (Akr1d1, Nr1h5) while genes with a gain in amplitude were linked with cholesterol uptake (Abcb1b) and internal trafficking (Stard3, Npc1) as well as TAG uptake (Vldlr). Genes with expression profile phase effects were more diverse in their functions. Phase-delayed genes were associated with cholesterol metabolization (Cyp8b1) and bile acid secretion (Cyp7a1, Abcb11), FA biosynthesis (Fads2), transport (Cpt1a) and elongation (Elovl2) (Supplementary file 4).

Focusing on cholesterol metabolism, we divided this pathway into three categories: uptake, biosynthesis, and bile acid secretion. Averaged diurnal gene expression from both groups was rhythmic with a strong mesor DOWN effect in MMI mice (Fig. 4E). Total liver cholesterol was reduced and arrhythmic in MMI mice compared to the CON group. On the other hand, serum total cholesterol was elevated in MMI compared to CON mice, indicating impairment in liver cholesterol uptake—as previously suggested. Liver TAG was less responsive to a low-TH state and showed a mesor DOWN effect (Fig. 4F).

Our findings show that a low-TH state results in distinct alterations in liver transcriptome rhythms. Rhythms in lipid and cholesterol metabolism-associated gene programs are affected by low TH levels and translate into differences mostly in cholesterol levels.

When assessing liver TH state across studies—particularly those based in clinical settings—it is not always possible to consider temporal dynamics in liver physiology. Therefore, it would be helpful to identify molecular markers that indicate liver TH state largely independent of sampling time. However, unlike what we had previously shown for high-TH conditions, no robust DEGs were identified in livers of MMI mice (Fig. 2B). While this would preclude reliable data analysis independent of sampling time, we circumvented this limitation by filtering our dataset for transcripts with marked TH state-dependent mesor effects in expression. This would allow for reliable detection of TH state at any time point if sampling time were consistent across experimental conditions. We identified 516 genes that showed consistent mesor effects in response to changes in TH state for both MMI and T₃ conditions (Fig. 5A). Gene Set Enrichment Analysis (GSEA) for genes in the DOWN-MMI/UP-T₃ group showed enrichment for metabolic processes such as steroid, cholesterol, FA, and bile acid metabolism. Conversely, genes in the UP-MMI/DOWN-T₃ group were associated with processes involved in FA, lipoprotein, carbohydrate, and xenobiotic metabolism, amongst others. Interestingly, glucose Glucose/glycogen metabolism was exclusively enriched in the UP in MMI/down in T₃ group (Fig. 5B,C; Supplementary file 5), suggesting a TH-driven overall inhibition of these pathways in line with our previous observations¹¹. Averaged expression of all genes pertaining to similar biological processes across the day showed that lipid metabolism—mostly comprised of FA biosynthesis, TAG uptake and biosynthesis—genes were highly responsive to TH levels, being up- and downregulated in T₃ and MMI mice, respectively. Genes associated with FA catabolism were up- and downregulated in MMI and T₃ mice, respectively. As already indicated in the GSEA, cholesterol and bile acid metabolism were markedly affected as MMI and T₃ mice showed down- and upregulation, respectively (Fig. 5C; Supplementary file 5). Validation of selected metabolic genes identified in the microarray was confirmed by qPCR (Supplementary Fig. 4).

Absolute mesor change was used to rank each gene for each comparison (CON vs. T3 and CON vs. MMI). The top-4 genes (Tlcd2, Stim2, Sdr9c7, Akr1c19) are shown (Fig. 5D) recapitulating low or high TH states at any fixed timepoint across the day. Independent qPCR validation confirmed the array findings (Supplementary Fig. 4). The fact that these genes are robustly rhythmic across conditions further emphasizes that sampling time should be kept consistent even in one-timepoint sampling studies (Fig. 5D; Supplementary file 5). To test the predictive efficiency of our approach, we used the published T₃ atlas that combines transcriptome and TH receptor ChipSeq data from different tissues¹⁶. We identified an overlap of ca. 70% (161 of 224 genes) between the two datasets, further supporting the validity of our approach. We further validated the predictive power of our tuning approach by comparing it against two studies. The first one¹⁷ consisted of rending mice hyperthyroid by supplementing T₄ in the drinking water for two weeks, followed by two weeks of recovery. From the TH-responsive DEGs identified, 42 overlapped with our tuning DEGs (20% overlap). The second study¹⁸ rendered THRβ^WT and THRβ^KO mice hypothyroid for four weeks, followed by one week of T₃ supplementation. This experiment design allowed the identification of THRβ-responsive genes. From our DEGs tunning set, a total of 128 genes were identified, representing 60% of overlap (Supplementary file 5).

Collectively, the predictive power of the tuning DEG set was validated in public datasets and represents an interesting gene set to estimate TH liver state quantitatively.

Two- to three-months-old male C57BL/6 J mice (Janvier Labs, Germany) were housed in groups of three under a 12-h light, 12-h dark (LD, ~ 300 lx) cycle at 22 ± 2 °C and a relative humidity of 60 ± 5% with ad-libitum access to food and water. Mice were treated with methimazole (MMI, 0.1%, Sigma-Aldrich, USA) potassium perchlorate (0.2%) and sucralose (1 tablet per 50 ml, Tafelsüss, Borchers) for three weeks.

During the treatment period, mice were monitored for body weight individually and food and water intake per cage. All in vivo experiments were ethically approved by the Animal Health and Care Committee of the Government of Schleswig–Holstein and were reported according to the ARRIVE international guidelines. All methods were carried out in accordance with relevant guidelines and regulations. The sample size was calculated using G-power software (version 3.1) and is shown as biological replicates in all graphs.

Euthanasia was carried out using cervical dislocation and tissues were collected every 4 h. Night experiments were carried out under dim red light. Tissues were immediately placed on dry ice and stored at  − 80 °C until further processing. Blood samples were collected from the trunk, and clotting was allowed for 20 min at room temperature. Serum was obtained after centrifugation at 2500 rpm, 30 min, 4 °C, and samples stored at  − 20 °C.

Serum quantification of T₃ and T₄ was performed using commercially available kits (NovaTec, Leinfelden-Echterdingen, DNOV053, Germany for T₃ and DRG Diagnostics, Marburg, EIA-1781, Germany for T₄) following the manufactures’ instructions.

TAG and total cholesterol evaluation were processed according to the manufacturer’s instructions (Sigma-Aldrich, MAK266 for TAG and Cell Biolabs, San Diego, USA, STA 384 for cholesterol).

Core body temperature and locomotor activity were monitored in a subset of single-housed animals using wireless transponders (E-mitters, Starr Life Sciences, Oakmont, USA). Probes were transplanted into the abdominal cavity of mice 7 days before starting the drinking water treatment. During the treatment period, mice were recorded once per week for at least two consecutive days. Recordings were registered in 1-min intervals using the Vital View software (Starr Life Sciences). Temperature and activity data were averaged over two consecutive days (treatment days: 19/20) and plotted in 60-min bins.

An open-circuit indirect calorimetry system (TSE PhenoMaster, TSE Systems, TSE Systems, USA) was used to determine respiratory quotient (RQ = carbon dioxide produced/oxygen consumed) and energy expenditure in a subset of single-housed mice during drinking water treatment. Mice were acclimatized to the system for one week prior to starting the measurement. Monitoring of oxygen consumption, water intake as well as activity took place simultaneously in 20-min bins. VO₂ and RQ profiles were averaged over two consecutive days (treatment days: 19/20) and plotted in 60-min bins. Energy expenditure was estimated by determining the caloric equivalent according to Heldmaier³⁶: heat production (mW) = (4.44 + 1.43*RQ)*VO₂ (ml O₂/h).

Total RNA was extracted using TRIzol (Thermofisher, Waltham, USA) and the Direct-zol RNA Miniprep kit (Zymo Research, Irvine, USA) according to the manufacturer’s instructions. Genome-wide expression analyses was performed using Clariom S arrays (Thermo Fisher Scientific) using 100 ng RNA of each sample according to the manufacturer’s recommendations (WT Plus Kit, Thermo Fisher Scientific). Data were analyzed using Transcriptome Analyses Console (Thermo Fisher Scientific, version 4.0) and expressed in log₂ values. Sample MMI_ZT06_a and MMI_ZT18_d were removed from the data due to low quality.

cDNA was generated by reverse transcription of total RNA using the high-capacity cDNA reverse transcription kit (Applied Biosystems, USA). Gene expression was determined by qPCR using the Go Taq qPCR master mix kit (Promega, USA) and a CFX-96 thermocycler (Bio-Rad, USA). Relative mRNA levels were obtained by analyzing gene expression with the Pfaffl Eq. ³⁷. Eef1α was used as a normalizer due to its robust expression across time and groups. Primer sequences are provided in supplementary file 5.

To identify global DEGs, all temporal data from each group was considered and analyzed by Student’s t test and corrected for false discovery rates (FDR < 0.1). Up- or downregulated DEGs were considered when a threshold of 1.5-fold (0.58 in log₂ values) regulation was met. As multiple probes can target a single gene, we curated the data to remove ambiguous genes. To identify DEGs at specific time points (ZTs—Zeitgeber time; ZT0 = “lights on”), the procedure described above for each ZT was performed separately. Time-independent DEGs were identified by finding consistent gene expression pattern across all ZTs.

To identify probes that showed diurnal (i.e., 24-h) oscillations, we employed the non-parametric JTK_CYCLE algorithm¹⁴ in the Metacycle package³⁸ with a set period of 24 h and an adjusted p-value (ADJ.P) cut-off of 0.05. Phase and amplitude parameter estimates from CircaSingle were used for rose plot visualizations¹⁵. To directly compare rhythm parameters (mesor and amplitude) in gene expression profiles between T₃ and CON, CircaCompare fits were used irrespective of rhythmicity thresholds. Phase comparisons were only performed when a gene was considered as rhythmic in both conditions (p < 0.05) as previously described¹¹. Temporal profiles were made using geom_smooth (ggplot2 package), method “lm”, and formula = y ~ 1 + sin(2*pi*x/24) + cos(2*pi*x/24). Small differences in rhythmic parameters can be present due to the different sine curve fitting between CircaCompare and ggplot curve fit.

Functional enrichment analysis of DEGs was performed using the Gene Ontology (GO) annotations for Biological Processes on the Database for Annotation, Visualization, and Integrated Discovery software (DAVID 6.8³⁹). Processes were considered significant for a biological process containing at least 5 genes (gene count) and a p-value < 0.05. To remove the redundancy of GSEA, we applied the REVIGO algorithm⁴⁰ using default conditions and a reduction of 0.5. For enrichment analyses from gene sets containing less than 100 genes, biological processes containing at least 2 genes were included. Overall gene expression evaluation of a given biological process was performed by normalizing each timepoint by CON mesor.

Samples were only excluded upon technical failure. Data from ZT0-12 were considered as light phase and from ZT 12–24 as dark phase. Day vs. night analyzes were performed by averaging the data and comparing using unpaired Student’s t test with Welch correction. Time-course data were analyzed by Two-way ANOVA for main treatment effects followed by Bonferroni post-test. When applicable, Two-Way for repeated measures was applied. Single timepoint data were evaluated by unpaired Student’s t test with Welch correction or Mann–Whitney test for parametric or non-parametric samples, respectively. ANOVA One-Way followed by Tukey was used for single timepoint data when CON, MMI, and T₃ groups were evaluated. Spearman’s correlation was used for all correlational analyzes. Analyzes were done in Prism 9.4 (GraphPad) and a p-value < 0.05 was used to reject the null hypothesis.

Statistical analyses were conducted using R 4.0.3 (R Foundation for Statistical Computing, Austria) or in Prism 9.4 (GraphPad). Rhythmicity was calculated using the JTK_CYLCE algorithm in meta2d, a function of the MetaCycle R package v.1.2.0. Rhythmic features were calculated and compared pairwise among the groups using the CircaCompare R package v.0.1.1. Data visualization was performed using the ggplot2 R package v.3.3.5, eulerr R package v.6.1.1, UpSetR v. 1.4.0, pheatmap v. 1.0.12, and Prism 9.4 (GraphPad).

Data from CON and T₃ were used and properly disclosed in each section. All experimental data from¹¹ was extracted from the Figshare depository, unless specified. Microarray data was extracted from Gene Expression Omnibus (GEO) database (GSE199998↗). All experimental data are already deposited in public depositories. No additional data was generated.