Abstract

Background & Aims

   Liver homeostasis is ensured in part by time-of-day-dependent
   processes, many of them being paced by the molecular circadian clock.
   Liver functions are compromised in metabolic dysfunction-associated
   steatotic liver disease (MASLD) and metabolic dysfunction-associated
   steatohepatitis (MASH), and clock disruption increases susceptibility
   to MASLD progression in rodent models. We therefore investigated
   whether the time-of-day-dependent transcriptome and metabolome are
   significantly altered in human steatotic and MASH livers.

Methods

   Liver biopsies, collected within an 8 h-window from a carefully
   phenotyped cohort of 290 patients and histologically diagnosed to be
   either normal, steatotic or MASH hepatic tissues, were analyzed by RNA
   sequencing and unbiased metabolomic approaches. Time-of-day-dependent
   gene expression patterns and metabolomes were identified and compared
   between histologically normal, steatotic and MASH livers.

Results

   Herein, we provide a first-of-its-kind report of a daytime-resolved
   human liver transcriptome-metabolome and associated alterations in
   MASLD. Transcriptomic analysis showed a robustness of core molecular
   clock components in steatotic and MASH livers. It also revealed
   stage-specific, time-of-day-dependent alterations of hundreds of
   transcripts involved in cell-to-cell communication, intracellular
   signaling and metabolism. Similarly, rhythmic amino acid and lipid
   metabolomes were affected in pathological livers. Both TNFα and PPARγ
   signaling were predicted as important contributors to altered
   rhythmicity.

Conclusion

   MASLD progression to MASH perturbs time-of-day-dependent processes in
   human livers, while the differential expression of core molecular clock
   components is maintained.

Impact and implications

   This work characterizes the rhythmic patterns of the transcriptome and
   metabolome in the human liver. Using a cohort of well-phenotyped
   patients (n = 290) for whom the time-of-day at biopsy collection was
   known, we show that time-of-day variations observed in histologically
   normal livers are gradually perturbed in liver steatosis and metabolic
   dysfunction-associated steatohepatitis. Importantly, these
   observations, albeit obtained across a restricted time window, provide
   further support for preclinical studies demonstrating alterations of
   rhythmic patterns in diseased livers. On a practical note, this study
   indicates the importance of considering time-of-day as a critical
   biological variable which may significantly affect data interpretation
   in animal and human studies of liver diseases.

   Keywords: Liver homeostasis, daytime rhythmicity, gene expression,
   metabolomic, transcriptomic, NAFLD, NASH

Graphical abstract

   [63]graphic file with name ga1.jpg
   [64]Open in a new tab

Highlights

     * •
       Temporal patterns in mRNA transcripts and metabolites are observed
       in the human liver.
     * •
       Time-of-day-dependent patterns are lost or gained in human
       steatotic and MASH livers.
     * •
       Core clock gene expression is resilient to disease progression.
     * •
       Time-of-day variations in lipid and amino acid metabolites are most
       significant.

Introduction

   Normal tissue homeostasis requires a precisely timed expression of
   genes and proteins around the clock and its alignment with cycles of
   light/dark exposure, feeding periods and physical activity. The central
   clock, located in the suprachiasmatic nucleus, is light-entrained and
   connects with peripheral tissues to synchronize clock oscillators in
   these tissues. However, peripheral tissue clocks can operate
   autonomously, i.e. independently of the central hypothalamic clock. For
   example, the major Zeitgeber (“time giver”) setting the liver clock is
   food intake/nutrient availability rather than daylight.[65]^1^,[66]^2
   Studies of molecular mechanisms controlling these circadian regulations
   in nocturnal rodents generated a global picture defining a universal
   molecular clock machinery.[67]^3 This cell-autonomous circadian core
   clock is made of two autoregulatory loops comprising 14 transcription
   factors, encoded by so-called core clock genes (CCGs). Heterodimeric
   BMAL1 (ARNTL) and CLOCK (or NPAS2) transcriptional activators and PER
   and CRY transcriptional repressors, along with the nuclear receptors
   RORs, REV-ERBα and β constitute interlocked
   transcriptional-translational feedback loops. These interlocked
   transcriptional-translational feedback loops define a cell-autonomous
   clock machinery which controls clock output genes,[68]^4 in turn
   regulating multiple cellular functions.[69]^3

   Because most primates (including humans) are diurnal, there are likely
   important differences from rodents in circadian regulation that have
   yet to be explored. A limited number of human time-of-day-resolved
   transcriptomes is available, especially for internal organs.
   Transcriptomes from whole blood,[70]^5^,[71]^6 peripheral blood
   mononuclear cells,[72]^7 skin,[73]^8 subcutaneous white adipose
   tissue,[74]^9^,[75]^10 heart,[76]^11 or skeletal muscle[77][12],
   [78][13], [79][14] were analyzed for a relatively low number of
   individuals (n <30). A gene expression study in subcutaneous white
   adipose tissue and skin from 625 healthy volunteers enabled the
   identification of time-of-day-regulated genes strongly enriched in
   CCGs.[80]^15

   Circadian rhythm dyssynchrony is observed in, and likely causative of,
   various diseases such as obesity and its complications like metabolic
   dysfunction-associated steatotic liver disease (MASLD), formerly known
   as non-alcoholic fatty liver disease (NAFLD). MASLD is a spectrum of
   liver conditions characterized by hepatic steatosis combined with
   varying degrees of necroinflammation and excluding excessive alcohol
   consumption.[81]^16 Its more severe, yet generally asymptomatic form,
   termed metabolic dysfunction-associated steatohepatitis (MASH; formerly
   known as NASH, non-alcoholic steatohepatitis), may evolve towards liver
   fibrosis, cirrhosis and hepatocellular carcinoma. Large-scale patient
   cohort studies have reported differences in gene expression between
   disease stages without time-of-day information.[82][17], [83][18],
   [84][19], [85][20] While hepatic metabolomes and transcriptomes are
   deregulated in rodent models of MASH and fibrosis,[86]^21^,[87]^22
   whether this occurs similarly in human MASLD remains unknown. A
   24 h-circadian transcriptome atlas of 64 tissues from healthy baboons
   identified only a small set of robustly cycling genes in the liver,
   surprisingly not including CCGs.[88]^23 Rhythmic gene expression
   patterns were inferred from the analysis of tissues collected post
   mortem from 600 human donors. Relatively few genes (n = 648), including
   only a few CCGs, exhibited predicted time-of-day-dependent
   expression.[89]^24 Thus, both ethical and technical hurdles hinder the
   thorough investigation of time-of-day-dependent processes in healthy
   human liver and the deregulation thereof in MASLD. Importantly, this
   conclusion extends to the hepatic metabolome, which is clock-controlled
   and disturbed in various liver dysfunction models.[90]^22^,[91][25],
   [92][26], [93][27], [94][28] Disturbances in the hepatic
   chronometabolome observed in rodent MASLD models have not yet been
   reported in humans.[95]^29

   Considering these knowledge gaps, we asked whether the hepatic
   time-of-day-dependent transcriptome and metabolome are affected during
   MASLD progression. We leveraged a large cohort of morbidly obese
   patients undergoing bariatric surgery from whom liver biopsies were
   taken peri-operatively (Hôpital Universitaire de Lille [HUL]
   cohort[96]^18). Hepatic transcriptomes and metabolomes were obtained
   from a sub-cohort of 290 patients whose biopsies were histologically
   identified as either normal, steatotic, or MASH livers and for whom the
   exact time-of-day at biopsy was known ([97]Fig. 1A). In an original
   approach integrating multiple statistical tests, we provide the
   first-ever robust analysis of time-of-day-dependent gene expression and
   tissue metabolite abundance, as well as changes associated with the
   different stages of MASLD, in human livers.

Fig. 1.

   [98]Fig. 1
   [99]Open in a new tab

   Timed liver biopsies from a large cohort of humans with obesity. (A)
   Overall experimental strategy. (B) Decision tree to stratify the HUL
   sub-cohort. “HN” (histologically normal), “steatosis” (benign steatosis
   only) or “MASH” (steatosis+inflammation). (C) Inter-sample variation in
   gene expression. ANOVA was used to reveal the main sources of overall
   inter-sample variation. (D) Biopsy daytime distribution. (E) Gene
   expression analysis. A volcano plot was generated by comparing gene
   expression using DEseq2 from samples collected in the morning (AM) or
   in the afternoon (PM) regardless of the pathological state. X axis:
   log[2] (fold change), Y axis: -log[10](p values). (F) KEGG pathway
   enrichment analysis. Biological term enrichment was carried out using
   the 1,660 genes whose expression was significantly different as
   determined in (E) (FC >1.2, FDR <0.05). (E, F) Font size was adjusted
   for clarity purpose. FC, fold change; FDR, false discovery rate; HN,
   histologically normal; KEGG, Kyoto Encyclopedia of Genes and Genomes;
   MASH, metabolic dysfunction-associated steatohepatitis.

Materials and methods

Liver biopsies from the HUL cohort

   The HUL cohort, also known as the Biological Atlas of Severe Obesity
   (ABOS) cohort, was established in 2006 by the University Hospital of
   Lille, France (ClinicialTrials.gov: [100]NCT01129297) and comprised
   severely and morbidly obese patients visiting the Obesity Surgery
   Department. The study protocol conforms to the ethical guidelines of
   the 1975 Declaration of Helsinki. All patients of the cohort fulfilled
   criteria for, and were willing to undergo, bariatric weight-loss
   surgery (for details, see[101]^18). Written informed consent was
   obtained from each patient included in the study. The protocol required
   that patients fast from midnight until the time of surgery. During the
   surgical procedure, wedge biopsies were taken from the liver to be
   immediately snap-frozen and the exact time of the biopsy was noted. A
   total of >1,500 patients are currently included in the HUL cohort,
   amongst whom 319 were selected to build a sub-cohort with complete
   clinical, biometric parameters and a robust histological MASLD
   classification of quality-controlled biopsies eliminating all
   intermediary MASLD stages (see[102]^18 and [103]Fig. 1B for more
   details). Both transcriptomes and metabolomes were obtained for these
   319 patients with biopsy mass >100 mg. Out of these, the time-of-day at
   biopsy (ranging from 8am to 4pm) was known for 290 patients, who were
   included in this study. The main clinical and histological
   characteristics are provided in [104]Table 1. Patient clinical data
   shown in [105]Table 1 were analyzed using the package “gtsummary”
   (v1.5.0).

Table 1.

   Biometric and biochemical parameters of the HUL sub-cohort.
   MASLD group
     __________________________________________________________________

   p values
     __________________________________________________________________

   Variable N HN, n = 89 St., n = 122 MASH, n = 79 HN vs. St. HN vs. MASH
   St. vs. MASH
   Sex 290 0.003 0.008 0.7
    Female 75/89 (84%) 80/122 (66%) 49/79 (62%)
    Male 14/89 (16%) 42/122 (34%) 30/79 (38%)
   Age, yr 290 34.6 ± 11.4 41.8 ± 10.9 46.4 ± 10.4 <0.001 <0.001 0.050
   BMI 290 45.4 ± 6.8 47.3 ± 7.9 46.2 ± 7.9 0.2 0.5 0.2
   Steatosis score 290 <0.001 <0.001 <0.001
    0 89/89 (100%) 0/122 (0%) 0/79 (0%)
    1 0/89 (0%) 0/122 (0%) 18/79 (23%)
    2 0/89 (0%) 71/122 (58%) 29/79 (37%)
    3 0/89 (0%) 51/122 (42%) 32/79 (41%)
   Inflammation score 290 <0.001 <0.001 <0.001
    0 89/89 (100%) 70/122 (57%) 0/79 (0%)
    1 0/89 (0%) 44/122 (36%) 55/79 (70%)
    2 0/89 (0%) 8/122 (6.6%) 23/79 (29%)
     0/89 (0%) 0/122 (0%) 1/79 (1.3%)
   Ballooning score 290 0.54 <0.001 <0.001
    0 89/89 (100%) 115/122 (94%) 0/79 (0%)
    1 0/89 (0%) 7/122 (5.7%) 57/79 (72%)
    2 0/89 (0%) 0/122 (0%) 22/79 (28%)
   Fibrosis score (Kleiner) 282 0.009 <0.001 <0.001
    0 79/89 (89%) 78/118 (66%) 12/75 (16%)
    1 9/89 (10%) 28/118 (24%) 22/75 (29%)
    2 0/89 (0%) 7/118 (5.9%) 13/75 (17%)
    3 1/89 (1.1%) 5/118 (4.2%) 25/75 (33%)
    4 0/89 (0%) 0/118 (0%) 3/75 (4.0%)
   NAS score 290 <0.001 <0.001 <0.001
    0 89/89 (100%) 0/122 (0%) 0/79 (0%)
    2 0/89 (0%) 40/122 (33%) 0/79 (0%)
    3 0/89 (0%) 51/122 (42%) 11/79 (14%)
    4 0/89 (0%) 26/122 (21%) 19/79 (24%)
    5 0/89 (0%) 5/122 (4.1%) 29/79 (37%)
    6 0/89 (0%) 0/122 (0%) 17/79 (22%)
    7 0/89 (0%) 0/122 (0%) 2/79 (2.5%)
    8 0/89 (0%) 0/122 (0%) 1/79 (1.3%)
   HOMA-IR 279 3.8 ± 4.8 13.5 ± 58.8 24.2 ± 55.2 <0.001 <0.001 <0.001
   Biopsy time 290 0.13 0.3 0.4
    AM 62/89 (70%) 70/122 (57%) 46/79 (58%)
    PM 27/89 (30%) 52/122 (43%) 33/79 (42%)
   [106]Open in a new tab

   The main biometric, biochemical and liver histological features of
   selected patients are indicated. Continuous values are expressed as
   mean ± SD. Inter-group comparisons were performed using the unpaired
   Wilcoxon test for continuous variables (age, BMI, HOMA-IR) and Fisher’s
   exact test for the remaining categorical variables. HN, histologically
   normal; HOMA-IR, homeostatic model assessment for insulin resistance;
   MASH, metabolic dysfunction-associated steatohepatitis; NAS,
   NAFLD/MASLD activity score; St., steatosis.

Total RNA sequencing and data processing

   A detailed procedure can be found in the supplementary materials and
   methods.

Liver metabolomics by liquid chromatography-mass spectrometry

   All tissue samples were flash frozen and maintained at –80 °C until
   processing. Sample preparation was carried out as described
   previously[107]^30 at Metabolon, Inc. (Morrisville, NC, USA). A
   detailed procedure can be found in the supplementary materials and
   methods.

Bioinformatic analysis

Bulk RNA sequencing

   All analyses were carried out using RStudio (v1.4.1106) with R
   (v4.1.0). Data processing for differential expression as a function of
   time-of-day can be found in the supplementary materials and methods.

Single-cell RNA sequencing

   All analysis has been made under R (v 4.2.0). A detailed procedure for
   data extraction and processing can be found in the supplementary
   materials and methods.

Enrichment analysis

   Time-dependent gene lists were analyzed for enrichment of Kyoto
   Encyclopedia of Genes and Genomes (KEGG) pathways and gene ontology
   biological process terms using Metascape 3.5[108]^31 with default
   settings ([109]https://metascape.org/). A detailed procedure for data
   processing can be found in the supplementary materials and methods.

Data visualization and illustrations

   Graphs were generated as ∗.svg files using R packages mentioned above.
   Data were imported in CorelDraw2020 to assemble figures. Drawings in
   [110]Fig. 1A are from Renée Gordon, Victovoi, and Mikael Häggström,
   M.D. and were made available to the public domain via Wikimedia Commons
   with no restriction of use. Bubbleplots were generated in R studio
   using the ggplot2, plotly, reshape2, rcpp, and tidyverse packages as
   described in.[111]^32

Results

Time-of-day is a major factor affecting gene expression in human liver

   Human liver biopsies were collected from patients with obesity and
   undergoing bariatric surgery, for which the exact daytime of the liver
   biopsy was recorded ([112]Fig. 1A). Key clinical parameters of the 290
   patients are summarized in [113]Table 1. Based on histological features
   of liver biopsies (steatosis, hepatocyte ballooning, lobular
   inflammation), patients were grouped according to MASLD stages and
   labelled as histologically normal (HN), steatotic or MASH liver
   following the decision tree shown in [114]Fig. 1B. The proportion of
   men in this cohort increased with MASLD severity, rising from 16% (HN)
   to 38% (MASH) with an average Clinical Research Network NAFLD activity
   score rising from 0 to 5, respectively. Since sex is an important
   biological variable in this context,[115]^18 this was considered during
   further analysis (see below). Patients in the steatosis and MASH groups
   were slightly older than those in the HN group, and expectedly also had
   higher insulin resistance on average.

   The source of variation in gene expression levels was estimated by a
   multivariate analysis of variance ([116]Fig. 1C). The F-ratio (ratio of
   the between-group variance to the within-group variance) not only
   confirmed sex and group (i.e. MASLD stage) as the main sources of
   variation as previously reported,[117]^18 but very interestingly
   identified biopsy time (AM vs. PM) as the third most significant source
   of variation ([118]Fig. 1C). Age made only a minor contribution to
   signal variation. The distribution of biopsy times ([119]Fig. 1D)
   revealed a daytime window of about 8 h. Biopsies were predominantly
   (>60%) taken in the morning with a first peak around 9:30 AM and a
   second, lower peak around 2:30 PM, due to the logistical schedule of
   surgical interventions. There was, however, no significant difference
   in biopsy time distribution between histological groups ([120]Table 1).
   Since exclusion of biopsies collected between 11am and 1pm did not
   modify the outcomes of preliminary statistical analysis, the “AM“
   subgroup was defined as biopsies taken before noon (12:00), the “PM”
   subgroup as biopsies taken after noon.

   Gene expression profiles were thus compared between morning (AM) vs.
   afternoon (PM) samples. Differentially expressed genes were identified
   using DEseq2 independent of the MASLD status but correcting for sex.
   Transcript counts for 1,660 genes were significantly different
   (Benjamini-Hochberg adjusted p value <0.05) between AM or PM biopsies
   ([121]Fig. 1E). Among the 100 top hits were most of the CCGs (PER3,
   ARNTL/BMAL1, NPAS2, NR1D1, NR1D2, PER2, CRY1, PER1), clock-related
   genes (CIART, DBP, NFIL3) ([122]Fig. 1E) as well as circadian-regulated
   genes involved in lipid metabolism (PPARD, LIPG, LPIN2 …). Because
   patients were fasting from midnight irrespective of surgery time, genes
   implicated in hepatic gluconeogenesis (G6PC, PCK1, SGK2 …) were, as
   expected, expressed at a higher level in PM samples ([123]Fig. 1E and
   [124]Table S1). Contrary to nocturnal rodents, genes from the negative
   limb of the clock displayed lower expression in the afternoon (PER3,
   NR1D1, NR1D2, CIART …), whereas genes from the positive clock limb
   displayed higher expression in the afternoon (ARNTL/BMAL1, NPAS2 …)
   ([125]Fig. 1E). Globally, genes displaying AM vs. PM differential
   expression were significantly enriched for the KEGG pathways “circadian
   rhythm”, “PPAR signaling pathway”, carbohydrate and lipid metabolic
   pathways, as well as cellular architecture and communication, among
   others ([126]Fig. 1F).

   Thus, an 8-hour time frame allowed for the detection of significant
   changes in time-of-day-dependent liver gene expression, with a large
   proportion of transcripts functionally related to circadian
   rhythmicity.

Time-dependent genes vary between MASLD stages

   We next examined whether time-of-day-dependent distributions of gene
   expression would differ between the histological states “HN”,
   “steatosis” and “MASH”. In order to achieve statistical power and
   obtain robust and exhaustive lists of time-dependent gene expression
   over the available daytime window, we used three complementary
   statistical methods analyzing different aspects of gene expression
   distribution (differential expression, partial Spearman correlation,
   Kolmogorov-Smirnov test), the results of which were agglomerated by a
   Fisher test to yield a combined p value for each gene ([127]Fig. 2A-C
   and [128]Figs S1 and S2). The three types of analyses are graphically
   exemplified for the CCG ARNTL/BMAL1 ([129]Fig. 2A-C), which served as a
   positive control to validate our approach, as it is among the most
   highly time-dependent genes regardless of the histological group. First
   DEseq2, which relies on a negative binomial distribution of gene
   expression, was used to identify differential gene expression between
   two conditions (AM vs. PM as in [130]Fig. 1) ([131]Fig. 2A). Second,
   partial Spearman correlation was computed between ARNTL gene expression
   and time-of-day ([132]Fig. 2B). Both approaches integrated sex as a
   confounding factor. Third, the Kolmogorov-Smirnov test was employed to
   determine whether AM and PM ARNTL expression distribution followed a
   similar law and thus were similar in shape ([133]Fig. 2C). Finally, the
   Fisher combined probability test or “Fisher’s method” was used as a
   meta-analysis method for p value combination: individual raw p values
   resulting from each statistical test were agglomerated into a single p
   value per group ([134]Fig. 2D). The detected expression profile of
   ARNTL, of other CCGs ([135]Fig. S1) and of all other transcripts
   ([136]Fig. S2), clearly confirmed that the available time window was
   sufficient for robust time-of-day analysis of gene expression. A total
   of 1,427 genes with an absolute fold change greater than 1.2 (AM vs.
   PM) (false discovery rate [FDR] <0.01) was identified ([137]Fig. 2E).
   The vast majority of these time-dependent genes (TDGs) were strikingly
   distinct when comparing the three patient groups. Less than 10% (132
   genes) were indeed common to all three groups (“common TDGs”)
   ([138]Fig. 2E) and notably included most CCGs (ARNTL, NR1D1/2, NPAS2,
   CRY1, PER1/2/3, DBP, CIART) ([139]Table S1). TDG repartition outside of
   this core set was strongly unequal between groups, with ≈50% (558) of
   non-shared TDGs found in HN, ≈35% (392) in steatotic and less than 15%
   (177) in MASH livers. Along the same line, we found that AM to PM fold
   changes of common TDGs were, on average, decreased in steatotic and
   even more in MASH livers when compared to HN livers ([140]Fig. 2F).
   These differences are illustrated for a selection of TDGs with AM-PM
   differences either decreasing (PCAT18, ARMC4, SIK1B) or increasing in
   MASH (CYP4Z1, RHOBTB1, PPIAP71) ([141]Fig. 3A,B).

Fig. 2.

   [142]Fig. 2
   [143]Open in a new tab

   A multi-test method identifies time-of-day-dependent genes. BMAL1/ARNTL
   expression as a function of time and of the liver histological grade
   (A-C). (A) A violin plot is shown to illustrate DEseq2 results to
   compare AM vs. PM expression and correcting for sex as a confounding
   factor. (B) A dot plot with linear tendency lines is shown to
   illustrate partial Spearman correlation analysis between gene
   expression and biopsy daytime. (C) A density plot is shown to
   illustrate the output of the Kolmogorov-Smirnov test for comparison of
   AM or PM gene expression distributions. (D) Fisher’s agglomeration
   method. Frequency histograms of uncorrected Fisher p values are shown
   for each group and the first (colored) bar of each group histogram
   indicates genes with p <0.05. (E) TDG distribution within histological
   groups. A Venn diagram shows the unequal distribution of 1,427 TDGs (FC
   >1.2, FDR <0.01) amongst groups. Numbers indicate the numbers of
   identified transcripts. (F) Gene expression variation of 132 “common
   TDGs”. CCGs and other transcripts are indicated. ∗Indicates p <0.05
   using ANOVA and Fisher’s LSD post hoc test. CCGs, core clock genes; FC,
   fold change; FDR, false discovery rate; TDGs, time-dependent genes.

Fig. 3.

   [144]Fig. 3
   [145]Open in a new tab

   Example distributions of TDGs altered in MASH. (A) Violin plots showing
   AM vs PM gene expression variations for representative TDGs losing
   their time-dependency in MASH. (B) Violin plots showing AM vs. PM gene
   expression variations for representative TDGs gaining time-dependency
   in MASH. Statistical tests were Kruskal-Wallis tests followed by
   unpaired Wilcoxon post hoc test for AM-PM comparisons in each group (∗p
   <0.05, ∗∗p <0.01, ∗∗∗p <0.005, ∗∗∗∗p <0.001). MASH, metabolic
   dysfunction-associated steatohepatitis; TDGs, time-dependent genes;
   VST, variance-stabilizing transformation.

   Common TDGs were collectively enriched for KEGG terms like “circadian
   rhythms” as expected from the content in transcripts coding for CCGs,
   and for metabolic regulatory pathways like the PPAR and FoxO pathways
   ([146]Fig. 4, [147]Fig. S3), illustrated by genes such as S1PR1, G6PC1,
   PCK1, SGK2, FASN, AQP7 and PPARD. TDGs unique to the HN group were also
   enriched, albeit to a lesser extent, for pathways linked to circadian
   rhythm (notably including CLOCK) as well as to fatty acid and amino
   acid metabolism. The most highly represented pathway was “gap
   junctions” ([148]Figs 4A and S3), characterized by genes such as PDGFB,
   MAP2K1 and transcripts encoding for tubulins TUBA1C/8, TUBB, TUBB1/2B
   ([149]Table S1), suggesting that homeostasis of epithelial barrier
   integrity/permeability, which is known to be disturbed in
   MASLD,[150]^33 requires an oscillating expression of these genes in
   healthy conditions. Genes unique to the steatosis group ([151]Table S1)
   were mostly linked to metabolism of lipids and fatty acids (DGKG,
   PLA2G4B/5, LPIN2/3, ETNK2, ETNPPL, PLPP4, FADS1/2, CYP2C8, GDPD1, SCD)
   and also to metabolism of peptides and amino acids (DNMT3B, GCLM, SDS,
   PSAT1, GNMT, ALDOC, GPT2, CSAD, UPB1) ([152]Figs 4A and S3). Lastly,
   TDGs specific to the MASH group ([153]Table S1) were highly enriched
   for signaling by calcium, cAMP or neurotransmitters (ADRB2, DRD1, GRM1,
   NTRK1, NTSR1, P2RX7, RYR2, CACNA1H, SSTR5, TBXA2R, FFAR2, SUCNR1), as
   well as for lipolysis (ADRB2, IRS1, PNPLA2) ([154]Figs 4A and S3).The
   temporal pattern of gene expression in homeostatic conditions is thus
   strongly affected by the disease state and indicative of compromised
   cellular communication and metabolic pathways.

Fig. 4.

   [155]Fig. 4
   [156]Open in a new tab

   Term enrichment analysis of TDGs and upstream regulatory pathway
   prediction. (A) Biological term enrichment analysis. TDGs identified in
   [157]Fig. 2 were enriched for GO terms related to BP using Metascape.
   Significantly enriched GO-BP clusters were manually collapsed for
   visualization purposes and top hits are indicated. (B-D) Pathway
   activity ranking. Speed2 TDG enrichment for pathway signature genes.
   Each pathway is represented as a bar showing the mean rank of the query
   list. The “bar code” plot shows the distribution of genes from the
   query list in the ranked reference signatures. BP, biological process;
   GO, gene ontology; TDGs, time-dependent genes.

   Inferring upstream regulatory cues or altered biological processes may
   be achieved by comparing differentially expressed gene lists to
   consensus gene expression patterns induced by a given perturbagen
   ([158]Fig. 4B-D). Speed2 (Signaling Pathway Enrichment using
   Experimental Datasets[159]^34) analysis enables gene lists to be probed
   against ranked gene signatures for 16 signaling pathways, with the aim
   of identifying upstream signaling mediators. Ranked signatures
   suggested that cues in homeostatic (HN) conditions could be TGFβ, TNFα,
   oxidative stress, TLR and estrogen ([160]Fig. 4B). In steatosis and
   MASH conditions, this pattern shifted towards a more limited signaling
   pathway panel with similar statistical significance, which included
   either TLR and VEGF (steatosis) or TNFα and TLR (MASH) ([161]Fig. 4C,D,
   respectively). Although causative links cannot be proven, this data
   could reflect a loss of physiological rhythmic function(s) in steatotic
   and MASH livers, which in turn gain rhythmic functions associated to
   pathogenic immune and proliferative stimuli and responses.

MASLD stages correlates with time-of-day changes in liver metabolites

   Our results suggested that metabolic pathways are altered in a
   time-of-day-dependent manner as MASLD progresses, with an enrichment in
   amino acid- and lipid metabolism-regulating genes ([162]Fig. S3).
   Therefore, an unbiased tissue metabolomic study by liquid
   chromatography-mass spectrometry was performed on the same 290 liver
   samples. Similar to the gene expression analysis, a global approach was
   initially employed to evaluate overall time-of-day dependence of tissue
   metabolite levels regardless of the MASLD status. This global analysis
   identified ≈220 metabolites whose amounts were significantly different
   in AM and PM biopsies (DEseq2 corrected for sex, FDR <0.1)
   ([163]Fig. 5A, [164]Table S2). Visual inspection of the volcano plot
   highlighted intermediates of lipid β-oxidation (carnitine derivatives),
   amino acids (kynurenate, oxo-arginine …) as differentially detected in
   AM vs. PM livers ([165]Fig. 5A). It also confirmed the more marked
   fasting status of “PM” patients exhibiting an increased hepatic content
   in 3-hydroxybutyrate. A biological term enrichment analysis confirmed
   that the majority of the identified metabolites belonged to amino acid,
   lipid and fatty acid metabolic pathways ([166]Fig. 5B).

Fig. 5.

   [167]Fig. 5
   [168]Open in a new tab

   Identification of TDMs measured by liquid chromatography-mass
   spectrometry. (A) Volcano plot of metabolite time-of-day-dependent
   differential abundance in liver. Fold change values and corresponding p
   values obtained using DEseq2 considering sex as a confounding factor.
   (B) Enrichment of KEGG metabolic pathways. A biological term enrichment
   against the KEGG database was run using the differentially detected
   metabolites (TDMs) identified in (A) (DEseq2 adjusted p value <0.1).
   (C) p value agglomeration by Fisher’s method. As for gene expression
   analysis (see [169]Fig. 2), p value agglomeration from three separate
   statistical tests (DEseq2, partial Spearman correlation,
   Kolmogorov-Smirnov) was carried out using Fisher’s method and
   identified TDMs within each group. The resulting uncorrected Fisher p
   values are shown as frequency histograms and the first (colored) bar of
   each group histogram indicates metabolites with p <0.1. TDMs,
   time-dependent metabolites.

   To highlight a possible time-of-day differential representation of
   metabolites between MASLD groups, we again combined the three
   statistical approaches as described for gene expression analysis
   (DEseq2, Spearman correlation, Kolmogorov-Smirnov test) followed by
   Fisher’s agglomeration for a robust identification of time-dependent
   metabolites (TDMs) ([170]Fig. 5C). A total of 251 TDMs were identified
   using this method (combined FDR <0.1), out of which only 14 (6%) were
   common to all three MASLD groups ([171]Fig. 6A,B). These common TDMs
   included amino acids such as proline and threonine, several fatty acids
   and the ketone body component β-hydroxybutyrate ([172]Figs 6B and
   S4A-F). KEGG metabolic pathway enrichment analysis revealed that these
   common TDMs were most significantly associated with the metabolism of
   amino acids (arginine, proline, threonine) ([173]Figs 6B and S4A-C). In
   agreement with the identification of the PPAR pathway based on gene
   expression patterns ([174]Fig. 4) and the detection of
   3-hydroxybutyrate ([175]Fig. S4D), synthesis of ketone bodies was also
   identified as a relevant term ([176]Fig. 6B).

Fig. 6.

   [177]Fig. 6
   [178]Open in a new tab

   Distribution and characterization of TDMs. After correction for
   multiple testing, 238 metabolites with a Fisher FDR<0.1 were considered
   as robustly time-dependent. (A) Distribution of TDMs among liver
   histological groups. (B-E) KEGG metabolic pathway enrichment of common
   (central intersection) or group-specific TDMs using the online
   metabolomics analysis platform MetaboAnalyst. FDR, false discovery
   rate; TDMs, time-dependent metabolites.

   Among the 197 stage-specific TDMs, 31% were specific to HN, 47% to
   steatosis and 20% to MASH ([179]Fig. 6A and [180]Table S2). TDMs
   specific to HN livers, and thus lost at the steatosis and MASH stages,
   were mainly associated with the metabolism of sphingolipids
   ([181]Fig. 6C and [182]Table S2) such as CDP-choline and sphinganine
   ([183]Fig. S4G, H). Glycerophosphoethanolamines, glycerophosphocholines
   as well as derivatives of cholesterol, amino acids and pyrimidine were
   also identified as time-dependent in normal livers ([184]Table S2 and
   [185]Fig. S4I-L).

   A similar enrichment analysis identified amino acid metabolic pathways
   (branched-chain, sulfur-containing, arginine, taurine) ([186]Figs 6D
   and S5A,B and [187]Table S2) as time-dependent in steatotic livers.
   Visual inspection of steatosis TDMs also identified a carnitine
   precursor (N6,N6,N6 trimethyl-lysine, [188]Fig. S5C) and derivatives
   ([189]Table S2 and [190]Fig. S5D,E) which could reflect altered fatty
   acid oxidation activity. Finally, MASH-specific TDMs were enriched
   mainly for vitamin, glycan and glycosylphosphatidylinositol metabolic
   intermediates ([191]Figs 6E and S5, and [192]Table S2).

   Taken together, these analyses highlight the disruption of
   time-of-day-dependent bioactive phospholipid metabolism and of amino
   acid biotransformation pathways during MASLD progression. Intriguingly,
   PPARγ ligands of the linoleic acid class (9,10-dihydroxy-9-octadecenoic
   acid [DiHOME][193]^35 and 9- and 13-hydroxyoctadecadienoic acid
   [HODE][194]^36) displayed a differential abundance in AM vs. PM
   steatotic and MASH livers, with estimated concentrations in the
   10-100 μM range which are sufficient to activate PPARγ
   ([195]Fig. S5F,J).

Integrative analysis of time-dependent genes and metabolites

   We next performed an integrative analysis of the transcriptomic and
   metabolomic data at the pathway level using the KEGG database. This
   analysis combined TDGs and TDMs specific to either HN or MASH stages
   and common TDGs and TDMs, irrespective of their relative time-of-day
   direction of change, to identify associated transcriptomic and
   metabolomic conditions operating in normal and MASH livers
   ([196]Fig. 7A). TDGs and TDMs characterizing the HN stage were enriched
   for metabolic pathways related to lipid and amino acid metabolism,
   while most of them were not detected at the MASH stage, or with a
   decreased significance (arginine [Arg] and proline [Pro] metabolism,
   glycerophospholipid metabolism). Linoleic metabolism was associated
   with the MASH stage ([197]Fig. 7B,C). HN- or MASH-enriched pathways
   (glycerophospholipid and linoleic pathways, respectively) were further
   probed for daytime variation of associated TDGs and TDMs. The
   glycerophospholipid pathway was characterized by an increased abundance
   in AM livers of three out of four detected diacylglycerol (DAG) species
   specifically at the HN stage. Glycerophosphocholine (GPC) intermediates
   (CDP-choline, GPC) displayed stage-specific time-of-day variations that
   were not correlated to the occurrence of glycerophospholipid species
   (X-GPC) ([198]Fig. 7D). Higher abundance of DAG species in the morning
   did not correlate with HN TDG expression changes in transcripts
   encoding enzymes involved in this metabolic pathway, with the exception
   of the PNPLA3 gene. PNPLA3/adiponutrin has acyltransferase activity,
   increasing the formation of phosphatidic acid from lysophosphatidic
   acid, which may lead to more DAG synthesis. It might also reflect the
   hydrolytic activity of PNPLA3 on triacylglycerol molecules, favoring
   DAG accumulation. Along the same line, we compiled TDG and TDM data
   related to linoleic acid metabolism ([199]Fig. 7E). The increased
   abundance of the PPARγ ligands 9- and 13-HODE in the afternoon at the
   MASH stage was mirrored by the gene expression of ALOX15, a dioxygenase
   catalyzing the synthesis of these two hydroxyoctadecadienoic acids
   which was higher in the morning. A similar lack of correlation was
   observed between 9,10-DiHOME hepatic content and the expression of the
   linoleic acid-converting CYP2C8 at the steatosis stage. These results
   suggest that time-dependent metabolite variation in these pathways is
   delayed with respect to gene regulation, and/or controlled by
   post-transcriptional processes.

Fig. 7.

   [200]Fig. 7
   [201]Open in a new tab

   Integrative analysis of TDGs and TDMs. (A) Analysis strategy outline.
   Common and HN- or MASH-specific TDGs and TDMs were analyzed
   associatively for enrichment of KEGG metabolic pathways using the
   online metabolomics analysis platform MetaboAnalyst, revealing
   potential links between metabolic pathways as defined by time-dependent
   gene expression patterns and metabolomic profiling (B, C).
   Reconstitution (partial) of metabolic pathways displaying time-of-day
   dependency in the HN (D) or MASH (E) group. 9,10-DiHOME,
   9,10-dihydroxy-9-octadecenoic acid; 12,13-DiHOME,
   12,13-dihydroxy-9-octadecenoic acid; 9,10-EpOME,
   9(10)-epoxy-12Z-octadecenoic acid; 12(13)-EpOME,
   12(13)-epoxy-9Z-octadecenoic acid; 9-HODE, 9-hydroxyoctadecadienoic
   acid; 13-HODE, 13-hydroxyoctadecadienoic acid; 9-HpODE, 9-hydroperoxy
   octadecadienoic acid; 13-HpODE, 9-hydroperoxy octadecadienoic acid;
   CDP-choline, cytidine-diphosphate-choline; DAG, diacylglycerol; GPC,
   glycerophosphatidylcholine; HN, histologically normal; LPA,
   lysophosphatidic acid; MASH, metabolic dysfunction-associated
   steatohepatitis; PA, phosphatidic acid; PC, phosphatidylcholine; TDGs,
   TDGs, time-dependent genes; TDMs, time-dependent metabolites; X-GPC,
   acyl-conjugated GPC.

   Of note, mapping of human transcript expression to liver cell types
   using a reference single-cell RNA sequencing dataset[202]^37 suggested
   that the identified enzymatic pathways may follow a cell-specific
   expression pattern. They appear mainly restricted, but not limited to,
   hepatocytes. As an example, ALOX15 is detected in dendritic cells,
   whereas ALOX5 is also detected in monocytes, neutrophils and basophils
   ([203]Fig. S6). Therefore time-dependent metabolite variation may occur
   in either identical or distinct cell types, reflecting a functional
   compartmentalization. Yet these observations confirmed the presence of
   time-of day variation in hepatic gene expression and the hepatic
   metabolome, and identified several new oscillating metabolites at the
   MASH stage ([204]Fig. 7 and [205]Table S2).

Discussion

   Chronobiological studies require multiple replicates at 2-hour
   intervals over a total period of 24 or even 48 h, within a controlled
   environment including timed exposure to light and food. These
   conditions are not achievable in human studies, precluding the analysis
   of cyclic processes and particularly in internal organs.[206]^38 While
   human circadian rhythms are appreciated by genome-wide association
   studies and the phenotypic manifestation of disturbed cyclic processes
   such as sleep, light exposure and eating patterns,[207]^39 their study
   in healthy or pathological conditions is indeed hindered by ethical and
   technical constraints. Despite controlled experimental setup with
   regard to sleep behavior and food intake in previous studies,[208][5],
   [209][6], [210][7], [211][8], [212][9], [213][10], [214][11],
   [215][12], [216][13], [217][14] a limited number of samples of mostly
   healthy individuals were collected, thereby limiting data
   interpretation due to high inter-sample variability.

   Herein, we report the first-ever robust analysis of the
   time-of-day-resolved human liver transcriptome and associated liver
   tissue metabolites, using a 290-patient cohort. Although the available
   time window of liver biopsies was only about 8 h, this temporal window
   was sufficient to robustly identify TDGs and TDMs. Time-dependency of
   genes and metabolites was distinct between histologically defined MASLD
   groups. However, a small proportion of genes was identified as
   time-dependent in all three patient groups and included CCGs,
   indicating that the molecular clock is rather robust in pathological
   conditions. In contrast, the alignment of rhythmic biological processes
   such as intercellular communication (gap junctions) and metabolic
   regulations is disrupted upon MASLD progression. Interestingly, none of
   the detected TDGs in patients with MASH belonged to the human and mouse
   core set of MASH/fibrosis-associated genes,[218]^20 underlining the
   need to consider time as an important biological variable.
   Interestingly, the number of stage-specific TDGs ([219]Fig. 2E)
   decreases from the HN to the MASH stage, and concomitantly enriched
   pathways lessen ([220]Fig. S3), hinting at a loss of functional
   adaptability/(metabolic) flexibility. In high-fat diet-fed mice,
   transcripts gaining rhythmicity when compared to chow diet-fed mice are
   strongly enriched for glycerophospholid metabolism,[221]^22 as in
   steatotic patients ([222]Fig. S3), indicating convergent mechanisms for
   liver adaptation to dietary imbalance as often occurring in MASLD. At
   the cellular level, Ca^2+ fluxes are submitted to ultradian variations
   and coupled to metabolic regulations,[223]^40 and mishandled
   intracellular Ca^2+ stores in MASH can significantly impact parenchymal
   and non-parenchymal liver cellular functions.[224]^41 A number of genes
   encoding for Ca^2+ channel components or involved in intracellular
   Ca^2+ signaling (PKD1L1, TRPC1 P2RX7, FFAR2, ADRB2, CACNA1H, GRM1,
   NTRK1) exhibited differential AM vs. PM expression specifically in MASH
   livers. Thus, in addition to the lipid-induced dysfunction of
   endoplasmic reticulum Ca^2+ transport,[225]^42 de novo oscillation of
   the calcium handling process accompanies progression to MASH. CLOCK
   gene deletion affected metabolite oscillations in mouse livers[226]^25
   and we detected time-of-day-dependent differential expression of CLOCK
   exclusively in HN livers, which display a more diverse metabolic
   activity than MASH livers ([227]Fig. 7B). While hinting at a possible
   role of (the loss of) CLOCK, the examination of individual metabolites
   nevertheless showed little overlap between mouse CLOCK-dependent
   metabolites and HN-specific metabolites. This lack of clear concordance
   can be explained by species-specific mechanisms, distinct effects of
   gene deletion vs. loss of time-dependent expression and/or technical
   bias.

   While examining the coherence between MASLD state-specific liver
   transcriptomes and metabolomes, we observed little correlation between
   enzyme-encoding genes and metabolite abundance. This disconnection is
   not unprecedented and was also observed in a highly standardized mouse
   study which minimizes the variability typically observed in human
   samples.[228]^22 The narrow time window of our study may explain in
   part this lack of correlation as transcripts are likely to precede
   metabolite production, hence affecting our statistical approach for the
   PM sub-cohort. It may also indicate significant time-of-day-dependent
   translational control[229]^43 as well as post-translational
   modifications regulating enzyme activity which are not captured by our
   analysis.

   Finally, another point of convergence between our study and mouse
   studies is the differential abundance of linoleic acid derivatives and
   PPARγ ligands 9,10-DiHOME, 9-HODE and 13-HODE in steatotic and MASH
   livers, respectively. The PPARγ-encoding gene NR1C3/PPARG itself did
   not display significant oscillatory expression in the human liver,
   contrasting with high-fat diet-fed mouse livers.[230]^22 Whether the
   estimated concentrations of these compounds are indeed sufficient to
   differentially activate human liver PPARγ, which is mostly expressed in
   endothelial cells, hepatocytes and macrophages ([231]Fig. S6), and play
   a causative role in hepatic transcriptional reprograming requires
   further in-depth investigation.

   On the one hand, it is remarkable that many of the biological processes
   and pathways shown to be affected by MASLD in other studies, identified
   on the basis of changes in gene expression or metabolite abundance
   levels regardless of time, also display altered time-dependent
   expression profiles in our study. On the other hand, we identified many
   novel potential links between genes with deregulated timed expression
   and MASLD pathogenesis, which were not previously considered by
   standard analyses. It is probably the combination of both types of
   deregulations that underlies the deeply disturbed liver functions once
   MASH is declared. Conversely, some of the changes detected when
   time-of-day information is absent or ignored may turn out to be
   artefacts, as time-of-day as a biological variable might not be
   equilibrated between groups.

   Considering that the time-of-day dependence of
   transcript/protein/metabolite measurements was previously neglected or
   ignored in nearly all human MASLD studies, our findings reveal a
   significant impact of time-of-day on many relevant pathogenic
   processes. As such, differences found between sample groups in cohort
   studies might reflect a previously under-appreciated bias in sampling
   time between groups. Further investigation is needed in this regard,
   particularly when studying human material. In any case, a new
   generality should be that sampling daytimes (or Zeitgeber times) must
   be carefully recorded and included in post hoc analyses whenever
   possible.

   This study has both strengths and limitations. This first-of-its-kind
   study revealed time-of-day transcriptomic and metabolomic alterations
   in human livers as a function of the histologically proven MASLD stage.
   It used a large cohort allowing both the selection of biopsies to
   adequately encompass healthy control, steatosis and MASH cases and
   robust statistical analysis. There are however a number of limitations
   inherent to the observational nature of the study. The narrow time
   window for biopsy collection (8 h) precludes the assessment of a 24 h
   diurnal rhythmicity and hence of the integrity of the molecular clock.
   The unavoidable difference in fasting duration could also be a
   confounding factor.

   Our bulk RNA sequencing approach provided full coverage of the
   transcriptome, but precluded the identification of liver cell types
   expressing TDGs. We thus calculated a cell specificity index of TDGs in
   human resident parenchymal (PC, hepatocytes) and non-parenchymal (NPC)
   CD45^- cells. The tau (τ) index, calculated by a robust and simple
   method to assess cell type specificity, was used as a metric.[232]^44
   Using a reference single-cell RNA sequencing dataset for human livers
   collected at unknown times,[233]^37 τ could be calculated for 584
   transcripts out of the 1,481 identified TDGs ([234]Table S1). Using
   this metric, only 51 genes displayed a cell type-restricted expression
   pattern in CD45^- cells (τ >0.85), which was not limited to hepatocytes
   ([235]Fig. S7). The remaining 533 transcripts including CCGs could be
   mapped to two or more cell types. Of note, 15 out of the 51
   cell-restricted transcripts were mostly expressed in dendritic or
   Kupffer cells ([236]Fig. S6).

   Taken together, our bulk RNA sequencing approach detected a higher
   number of TDGs than previous efforts, including many pathways never
   previously reported as being time-dependent. This resource will serve
   as an important basis for further investigations that also consider the
   cell type specificity in time-of-day gene expression variation in human
   MASLD.

Financial support

   This work was supported by ANR (RHU PreciNASH 16-RHUS-0006, EGID
   ANR-10-LABX-0046 and DeCodeNASH ANR-20-CE14-0034). MJ was supported by
   Wallonie-Bruxelles International (WBI, Belgium, ref. SUB/2020/479801)
   and from European Association for the Study of the Liver (EASL, Sheila
   Sherlock fellowship). JTH and AB hold a European Research Council (ERC)
   grant (StG, Metabo3DC, contract number 101042759; CoG, OpiO, contract
   number 101043671 respectively). PL’s and BS’s teams are supported by
   Fondation pour la Recherche Médicale (Equipes labellisées FRM
   EQU202203014645 and FRM EQU202203014650 respectively).

Authors’ contributions

   Conceptualization: MJ, JTH, JV, FP, BS, PL; Software: MJ, JTH, JV, MD;
   Validation: MJ, JTH, JV, BD, VR, JE, BS, PL; Formal analysis: MJ, JTH,
   JV, JE, PL; Investigation: MJ, JTH, JV, BD, MD, AB; Resources: AB, PF,
   FP, BS, PL; Data curation: MJ, JTH, JV, BD, MD, VR, PL; Writing: MJ,
   JTH, JV, BS, PL; Visualization: MJ, JTH, JV, PL; Supervision: BS, PL;
   Project administration: BS, PL; Funding acquisition: FP, BS, PL, BD,
   VR.

Data availability statement

   Further information and requests for resources and data sets should be
   directed to and will be fulfilled by the corresponding author. Data
   exploration of the HUL cohort is still ongoing and restrictions apply
   to clinical data availability.

Disclaimer

   Funding sources were neither involved in the conduct of the research
   nor in the preparation of the article. They had no role in study
   design, in data collection, analysis and interpretation; in the writing
   of the report; and in the decision to submit the article for
   publication.

Conflicts of interest

   The authors of this study declare that they do not have any conflict of
   interest.

   Please refer to the accompanying ICMJE disclosure forms for further
   details.

Footnotes

   Author names in bold designate shared co-first authorship

   Supplementary data to this article can be found online at
   [237]https://doi.org/10.1016/j.jhepr.2023.100948.

Supplementary data

   The following are the supplementary data to this article:
   Multimedia component 1
   [238]mmc1.pdf^ (5.9MB, pdf)
   Multimedia component 2
   [239]mmc2.xlsx^ (201.6KB, xlsx)
   Multimedia component 3
   [240]mmc3.xlsx^ (5.8MB, xlsx)
   Multimedia component 4
   [241]mmc4.docx^ (35.5KB, docx)
   Multimedia component 5
   [242]mmc5.pdf^ (877.7KB, pdf)
   Multimedia component 6
   [243]mmc6.pdf^ (7.6MB, pdf)

References