Henry C. Hollis
Julian N. Francis
Ron C. Anafi- 1School of Biomedical Engineering and Health Systems, Drexel University, Philadelphia, PA, United States
- 2Department of Mathematics, Howard University, Washington, DC, United States
- 3Division of Sleep Medicine and Chronobiology and Sleep Institute, University of Pennsylvania, Philadelphia, PA, United States
Objective: Both critical illness and current care have been hypothesized to upset daily rhythms and impair molecular circadian function. However, the influence of critical illness on clock function in different tissues and on circadian output genes are unknown. Here we evaluate the effect of critical care and illness on transcription, focusing on the functional organization of the core circadian oscillator.
Methods: We downloaded RNAseq count data from the Genotype-Tissue Expression (GTEx) project. Treating mechanical ventilation as a marker for intensive care, we stratified samples into acute death (AD) and intensive care (IC) groups based on the documented Hardy Death Scale. We restricted our analysis to the 25 tissues with >50 samples in each group. Using the edgeR package and controlling for collection center, gender, and age, we identified transcripts differentially expressed between the AD and IC groups. Overrepresentation and enrichment methods were used to identify gene sets modulated by intensive care across tissues. For each tissue, we then calculated the delta clock correlation distance (ΔCCD), a comparative measure of the functional organization of the core circadian oscillator, in the both the AD and IC groups. The statistical significance of the ΔCCD was assessed by permutation, modifying a pre-existing R package to control for confounding variables.
Results: Intensive care, as marked by ventilation, significantly modulated the expression of thousands of genes. Transcripts that were modulated in ≥75% of tissues were enriched for genes involved in mitochondrial energetics, cellular stress, metabolism, and notably circadian regulation. Transcripts that were more markedly affected, in ≥10 tissues, were enriched for inflammation, complement and immune pathways. Oscillator organization, as assessed by ΔCCD, was significantly reduced in the intensive care group in 11/25 tissues.
Conclusion: Our findings support the hypothesis that patients in intensive care have impaired molecular circadian rhythms. Tissues involved in metabolism and energetics demonstrated the most marked changes in oscillator organization. In adipose tissue, there was a significant overlap between transcripts previously established to be modulated by sleep deprivation and fasting with those modulated by critical care. This work suggests that intensive care protocols that restore sleep/wake and nutritional rhythms may be of benefit.
Introduction
Daily or circadian cycles in biology modulate much of physiology and behavior (Takahashi et al., 2008). Driven by endogenous biological clocks, circadian rhythms coordinate internal cellular processes, both with the external environment and with each other. In both mice (Zhang et al., 2014) and humans (Ruben et al., 2018), nearly 50% of protein encoding transcripts show daily rhythms in one or more tissues. Moreover, tissue specific circadian output genes are highly enriched for disease associated genes and drug targets. Circadian modulation of metabolism and immune function appear particularly robust (Zhang et al., 2016; Haspel et al., 2020).
In almost every human cell, molecular circadian rhythms are generated by a transcriptional translational feedback loop (TTFL). (Takahashi et al., 2008). Transcriptional activators CLOCK and ARNTL bind the promoters of EBOX containing targets including the Period (PER1, PER2, and PER3) and Cryptochrome (CRY1 and CRY2) genes. After protein translation and a series of post-translational modifications, the PER and CRY proteins are actively translocated, return to the nucleus, and ultimately inhibit their own translation by the CLOCK:ARNTL complex (Kume et al., 1999). Over time CRY and PER protein degradation removes their repressive influence, allowing the cycle to begin anew the next day. This simple delayed feedback loop has multiple levels of control (Takahashi, 2016): The stability of the CRY proteins is modulated by the opposing action of two ubiquitin ligases: F-Box and leucine rich repeat proteins 3 and 21 (FBXL3 and FBXL21) (Hirano et al., 2013; Yoo et al., 2013). Phosphorylation by kinases such as casein kinase 1 epsilon (CK1ε), casein kinase 1 delta (CK1δ) and cyclin-dependent kinase 5 (CDK5), orchestrate PER protein turnover (Lowrey et al., 2000) along with interactions between PER and CRY and their translocation back to the nucleus (Brenna et al., 2019). Additional regulation is provided by interacting feedback loops (Preitner et al., 2002; Bugge et al., 2012), allosteric protein-protein interactions (Cal-Kayitmazbatir et al., 2020), and epigenetic modifications. (Etchegaray et al., 2003; Nakahata et al., 2008; Katada and Sassone-Corsi, 2010; Koike et al., 2012; Greco et al., 2020; Li et al., 2020). Notably these additional levels of control have been hypothesized to link circadian rhythms with cellular metabolic state (Aguilar-Arnal and Sassone-Corsi, 2013; Milev and Reddy, 2015; Patke et al., 2019).
The utility of these clocks likely depends on the ability to entrain them to the external environment and with each other. The “master” clock in the supra chiasmatic nuclei of the hypothalamus receives photic information and helps coordinate peripheral clocks through a combination of autonomic and hormonal signals (Bass, 2012). Other tissues are more easily entrained by different time giving signals or zeitgebers. While the SCN, for example, is closely entrained to the light-dark cycle, clocks in the liver and even skin are more tightly entrained to feeding rhythms (Vollmers et al., 2009; Wang et al., 2017). The sleep/wake cycle, perhaps the most outwardly obvious of all circadian regulated behaviors, itself contributes substantially to daily rhythms in transcription (Möller-Levet et al., 2013; Archer et al., 2014) and may synchronize peripheral clocks through autonomic and temperature signals (Anafi et al., 2013; Hoekstra et al., 2019).
Patients with need of intensive care are, almost by definition, subject to a host of severe inflammatory and metabolic stressors. The decision to admit patients to intensive care is based on a need for aggressive interventions and/or continuous monitoring. These patients are typically at risk for hypoxia or decreased systemic perfusion, often related to trauma, acute neurologic events, or profound illness with systemic inflammatory response (Smith and Nielsen, 1999). Indeed, critical care and illness have been found to influence tissue physiology to an extent that gene expression can predict subjects’ pre-mortem critical care status (Somekh et al., 2022). Disease related disruption of the hypothalamic-pituitary-stress axis, along with associated changes in body temperature, autonomic tone, and hormone secretion are well established, and all have been implicated in peripheral circadian entrainment (Nicolaides et al., 2014, 2017; Astiz et al., 2019). Moreover, vasopressors, steroids, sedatives, and other treatments often used in intensive care may directly modulate these entraining signals (Brainard et al., 2015; Astiz et al., 2019; Maas et al., 2020b; Haspel et al., 2021). Of equal importance, contemporary intensive care medicine takes place in an environment often devoid of normal time-entraining cues with poor lighting, absent or irregular enteric nutrition, constant sedation, and/or sleep disruption (Freedman et al., 2001; Gabor et al., 2003; Brainard et al., 2015). Thus, several groups have hypothesized that patients in intensive care are likely to have disrupted molecular clocks and this disruption may, in theory, compound the existing immune and metabolic insults (Telias and Wilcox, 2019; Jobanputra et al., 2020; Haspel et al., 2021). Initial studies have demonstrated disruption of sleep/wake and core-body temperature rhythms (Gazendam et al., 2013) among IC patients and absent or altered circadian patterns of serum melatonin and cortisol (Mundigler et al., 2002; Kiss et al., 2015; Maas et al., 2020b). Similar results have been obtained when comparing the temporal patterns in whole blood gene expression among IC patient cohorts and controls, both at the level of select core clock genes (Diaz et al., 2020) and the larger transcriptome (Maas et al., 2020a). Indeed, agglomerating data from multiple patients to create a composite temporal profile, Diaz et al. found that blood born rhythms in CLOCK, ARNTL, and PER2 detectable upon admission were lost after 1 week in the ICU.
The influence of critical illness and critical care on circadian and transcriptional regulation in other tissues, however, remains unclear. Of course, practical safety limitations prohibit repeated sampling of internal organs in ill (or healthy) patients. However, recently several tools have been developed to characterize rhythms and circadian function in agglomerated non-time stamped data.
Materials and methods
RNAseq count data was downloaded from the Genotype-Tissue Expression project (GTEx Analysis V8) (Aguet et al., 2020). Using Hardy Death Scale scores documented in the GTEx data, we divided subjects into two groups, separating subject who died in an acute manner from those who died on a ventilator. Of the subjects ventilated, the median ventilation time was 76.39 h (median absolute deviation ± 32.4 h). A Hardy Death score of one describes patients who died in a violent or sudden manner due to accident or trauma, while a score of two describes patients who were reasonably healthy until a sudden terminal medical event (i.e., cardiac arrest). Subjects with scores of one or two are estimated to have a terminal phase <1 h and these subjects were combined to constitute the acute death (AD) group. Subjects with a Hardy Death score of zero died on a ventilator and constitute the Intensive Care (IC) group. Subjects who did not fall into the above-mentioned Hardy Death Scale categories were discarded for this analysis (categories three and four). While the GTEx dataset includes data from multiple tissues from each subject, different tissues were sampled from each subject. As a result, some tissues had very few samples from the AD or IC groups. As our analysis is focused on assessing difference between the IC and AD groups, we restricted our analysis to the 25 tissues with greater than 50 samples in each (IC and AD) group.
The edgeR package in R (R version 4.1.0) was used to assess differential expression between the AD and IC groups. Read counts were normalized prior to eliminating data from subjects with Hardy scale three and four deaths. Lowly or non-expressed genes were filtered out as detailed in the edgeR documentation (Robinson et al., 2010; Chen et al., 2016). The design matrix for differential expression analysis included collection center, age, sex and manner of death (IC vs. AD). The age assigned to each subject was the mean of the 10-years age range provided by GTEx. Expression values for each transcript were fit to a negative binomial distribution with a generalized log-linear model as per the edgeR user guide (Chen et al., 2022). Differentially expressed genes were found by performing a quasi-likelihood F test (edgeR: glmQLFTest). Two distinct criteria for differential expression were invoked: The more strict criteria required both a statistically significant difference in expression between the AD and IC groups (Benjaminni Hochberg q value BH. q < 0.050) (Benjamini and Hochberg, 1995) and a marked absolute change in expression of >50%. The more lenient criteria required a statistically significant difference in expression between the AD and IC groups (BH.q < 0.05) but did not require a specific fold change cutoff. Two different gene sets were thus identified for downstream analysis: The first list included all transcripts that met the strict criteria for differential expression in ten or more tissues. The second list included all genes that met the more lenient criteria but in a much larger number of tissues (75%, or 19 or more tissues).
The EnrichR R package was used to perform overrepresentation analysis and identify annotated gene sets overrepresented among transcripts differentially expressed in the IC group (Chen et al., 2013; Kuleshov et al., 2016). Annotated gene sets from the Gene Ontology (GO) GO_Biological_Process_2021 and KEGG_2021_Human gene set libraries were analyzed.
Clock correlation distance (CCD) is a metric designed to measure the functional organization of the core circadian oscillator (Shilts et al., 2018). The CCD method compares a matrix of Spearman correlation coefficients between 12 core clock genes to a reference correlation matrix provided by Shilts et al. To compare oscillator organization between two groups, the CCD of the treatment group is subtracted from the CCD of the control group yielding a ΔCCD. The statistical significance of this metric is assessed by bootstrap, shuffling the condition labels and recording the number of instances a greater or equal ΔCCD is observed (Supplementary Figure S1). A measure of significance is obtained using the method of Phipson and Smyth (2010).
We calculated the ΔCCD between acute (AD) and intensive care (IC) patients for each tissue. However, the ΔCCD metric can also be influenced by confounding variables. To account for confounding due to different collection center distributions in the AD and IC groups we modified the bootstrap procedure used to assess statistical significance. The modified test shuffles the condition labels (AD vs. IC) among subjects from each collection center independently. It then joins the permutations created for each collection center to create a combined permutation. This process maintains the same collection center distribution in each of the permutations (Supplementary Figure S2) used for comparison. This process thus directly accounts for the role of collection center in assessing the significance of any observed difference in the ΔCCD. We extended this procedure to account for the influence of both subject sex and collection center. One thousand bootstrap permutations were used to generate a null distribution simultaneously including both covariates. An alternative method to control for confounding variables was also performed using Combat-Seq in the sva R package (Leek et al., 2021). In this case, covariate normalized expression measures (adjusting for age range, gender, and collection center) were directly processed with the original ΔCCD method. An overview of the full informatic experimental design in shown in Supplementary Figure S3.
Results
We identified 25 tissues with more than 50 samples in each of the AD and IC groups (Table 1). Among those 25 tissues, a mean of 421.2 samples were collected. These samples were collected from two primary GTEx collection centers. The median age of the subjects ranged from 50–62 depending on tissue. A summary of the data is given in Table 1.
TABLE 1. Summary statistics of included subjects. AD and IC groups refers to Acute Death and Intensive Care groups, respectively. Samples collected from centers B1 and C1. MAD refers to the median absolute deviation of a variable from the median.
Using edgeR and a model that accounts for confounding categorical factors including sample collection center, subject sex, and age range, we analyzed the expression data for each tissue. We first identified transcripts that showed marked and significant changes in expression when comparing the AD and IC groups. Requiring strict statistical significance after adjusting for multiple testing (BH.q < 0.05) and an expression fold change of at least 50% between the two groups, we identified hundreds to thousands of differentially expressed genes in each tissue (Supplementary File S1). We then identified a set of these genes that were markedly differentially expressed in 10 or more tissues (Figure 1A). Notably these transcripts showed a consistent direction of expression change across tissues (Figure 1B). We used EnrichR to identify GO Biological Processes or KEGG annotated gene sets overrepresented among genes that markedly differed between the AD and IC group across ≥10 tissues. Consistent with the expectation that patients requiring critical care were under profound physiologic stress, pathways related to inflammation and complement activation along with metabolic pathways were overrepresented.
FIGURE 1. Analysis of genes strongly affected by critical care and illness. We analyzed the set of genes that were both significantly (BH.q < 0.05) and strongly (fold change of >1.5 or log2 fold change 0.585) differentially expressed between the IC and AD groups across tissues. (A) Histogram of the number of tissues for which any given gene meets these criteria. (B) Log2 fold changes of genes significantly differentially expressed between AD and IC groups (in 22 or more tissues) shown across tissue type. Transcripts appear to be consistently up or down regulated across tissue type. (C) GO BP and (D) KEGG pathways overrepresented in IC differentially expressed genes that met the above criteria in 10 or more tissues.
Applying more lenient criteria and identifying transcripts with a significant difference in expression between the AD and the IC groups (BH.q < 0.05), but not requiring a specific fold change, we again identified thousands of transcripts in each tissue (Figures 2A,B). We then identified a set of transcripts that, although showing these more modest differences between the AD and IC groups, were affected in an even greater number (≥75%) of tissues. Overrepresented KEGG and GO BP pathways identified among these genes reflected widespread changes in cellular metabolism, protein, and nucleic acid processing, mitophagy, and autophagy. Notably the KEGG pathway (Figure 2D) for circadian rhythm regulation was among the most significantly overrepresented pathways (BH.Q < 0.05).
FIGURE 2. Analysis of genes significantly affected by critical care across a wide variety of tissues. We analyzed the set of genes that were significantly (BH.q < 0.05) differentially expressed in the IC and AD groups. (A) Histogram of the number of tissues for which any given gene meets these criteria. (B) Log2 fold changes of genes significantly differentially expressed between AD and IC groups (in 25 tissues) shown across tissue type. Transcripts appear to be consistently up or down regulated across tissue type. (C) GO BP pathways and (D) KEGG pathways enriched in IC differentially expressed genes that met the above criteria in ≥75% of tissues. Note appearance of circadian rhythm among the most significantly (p < 0.0003, BH. q < 0.01) enriched pathways.
We defined a set of “core clock” genes combining the exemplar core clock components in (Anafi et al., 2014) with the clock genes used in (Shilts et al., 2018). In most tissues, the expression of genes in the core clock mechanism was modestly but significantly different in the IC vs. AD groups (Figure 3A). To help discern if the changes observed in circadian rhythm regulation reflect a significant change in oscillator function, we computed the ΔCCD, a comparative measure of the functional organization of the core circadian clock. The ΔCCD metric is based on a well conserved pattern of clock gene co-variance (Figure 3B Left). Tissues with a higher ΔCCD value display a greater disruption in the functional organization of the core oscillator when comparing the IC and AD groups. After pre-processing the expression data with Combat-Seq to adjust for age range, sex, and collection center, the ΔCCD suggested a significant decrease in oscillator function in the IC group in 14/25 tissues (Bonferroni corrected p ≤ 0.05). We then applied the more stringent, custom ΔCCD bootstrap procedure directly to the edgeR normalized count data. In 11/25 tissues, IC patients demonstrated disrupted core clock organization (Bonferroni corrected p ≤ 0.05). We did not have enough samples to adjust our liver analysis for gender. Thus, the liver analysis accounted only for collection center. Tissues with a higher ΔCCD showed a greater disruption in oscillator organization between IC and AD groups. Example CCD matrices showing clock organization in visceral adiposes, liver, and lung tissue–in both the IC and AD groups - are shown in Figure 3B (right). While the organization of the core clock in visceral adipose was markedly decreased in the IC group, liver provides an example of a tissue demonstrating a more moderate effect. Interestingly, while the IC group was defined by ventilator use, in the lung, core oscillator organization was relatively unchanged between the IC and AD groups. Regardless of the statistical method invoked, tissues associated with digestion and metabolism (Table 2) showed the largest decrease in clock function among IC patients.
FIGURE 3. Effect of critical care on core clock organization and circadian output genes. A set of core clock genes was taken from the literature (Anafi et al., 2014; Shilts et al., 2018). (A) Log2 fold changes of core circadian clock genes significantly differentially expressed between AD and IC groups shown across tissues. Note that in most tissues, most core clock genes show a very small fold change. CCD matrices assess core clock organization and co-expression. (B) Left: the reference correlation matrix, used to calculate the CCD metric. Right: Example Spearman correlation matrices of core circadian clock genes in visceral adipose, liver and lung tissue. Clock function in visceral adipose is markedly decreased in the IC group, liver tissue shows a more moderate effect. ΔCCD analysis does not suggest a change in pulmonary oscillator organization following IC. (C) We ranked genes based on the log fold change between the AD and IC groups in subcutaneous adipose tissue. This was compared to lists of core clock and robust circadian output genes using fGSEA. Enrichment plot for core clock set (top), does not show significant enrichment. Enrichment plot for the robust circadian output genes shows significant enrichment (p < 0.01) as clock outputs were frequently ranked among the genes most strongly influenced by IC.
TABLE 2. Influence of IC on clock correlation structure across tissues. Twelve core circadian clock genes have a well conserved signature correlation matrix in healthy mammalian tissue.
While the absolute expression changes in core clock genes were more modest, the altered correlation pattern captured by the ΔCCD suggest a functional change. In addition to defining a set of core clock genes, we defined a set of robust circadian output genes. Robust circadian outputs are human transcripts orthologous to transcripts that demonstrated statistically significant (JTK Cycle bHQ<0.05) circadian rhythms in ≥75% of tissues in the Mouse Circadian Atlas (Zhang et al., 2014). We then used Fischer’s exact test to see if either core clock or robust circadian output genes were overrepresented among those transcripts that were markedly (BH.q < 0.05 and >50% fold change) different between the IC and AD groups (Table 3). While core clock genes were not overrepresented among markedly IC affected genes in any tissue, 10 tissues showed an overrepresentation of robust clock output genes in IC affected genes (BH.q < 0.1)—consistent with altered oscillator function. In a related analysis we ranked the transcripts in each tissue by the log fold change observed between the IC and AD groups and used fGSEA (Korotkevich et al., 2021) to determine if core clock genes or robust circadian outputs were enriched among the transcripts most effected by critical illness (Figure 3C). The enrichment analysis, like the overrepresentation analysis, showed a stronger influence of IC on the expression of circadian outputs as compared to the absolute expression of core clock transcripts.
TABLE 3. Overexpression analysis of core clock and circadian output pathways in IC affected genes.
Discussion
Multiple studies have shown than circadian rhythms in plasma cortisol, melatonin, and blood-based gene expression are disrupted among IC patients. Light, noise, steroids, vasopressors anesthetics, irregular nutrition, and benzodiazepines, have all been hypothesized to disrupt circadian rhythms, weaken entraining signals, or disrupt behaviors that may synchronize rhythms (Brainard et al., 2015; Maas et al., 2020b). However, prior work has generally been limited to blood and little is known about internal tissue physiology.
Here we show that thousands of genes are differentially expressed (BH.q < 0.05) across tissues in patients who died on a mechanical ventilation as compared to those who died suddenly. Consistent with the extreme physiologic stress of critical illness and critical care, the subset of genes modulated in ≥75% of tissues is enriched for pathways involving mitochondrial energetics, cellular stress, and immune function. Notably circadian regulation was also among the most enriched pathways. While the expression changes in core circadian clock genes observed between the IC and AD groups were relatively modest, we used the ΔCCD metric to assess the influence of these changes on core clock organization. The organization of core oscillator appeared to be significantly disrupted 11/25 tissues after using our modified ΔCCD procedure to correct for covariates. Among the most highly dysregulated tissues were those involved in digestion and metabolism, including various types of adipose, thyroid, esophageal, and colon tissue.
It is important, however, to note the limits of this analysis. The de-identified GTEx data available for public download does not provide time-course information. As a result, our study does not purport to describe the circadian transcriptome. Even if we were to obtain time-of-death information for each autopsy sample, the ability to agglomerate these data to re-construct a circadian transcriptome would require that the IC subjects are all well entrained to the shared external environment. Indeed, it is the assumption of shared entrainment that justifies using a 2:00 a.m. sample from one patient as follow up for a 1:00 a.m. sample from another. Yet this assumption is clearly in doubt in the IC cohort. Thus, unlike previous assessments of blood-based rhythms in IC patients, we did not attempt to assess the entrainment.
The core circadian clock is composed of set of activators (e.g. CLOCK and ARNTL) that induce the transcription of a set of repressors (e.g. PER and CRY family genes). This internal architecture imparts a well conserved temporal ordering among core clock genes with the activators sharing similar temporal expression patterns and being out of phase with repressors. As a result, the expression of activators tends to be positively correlated across samples and negatively correlated with repressors. For example, even combining data from night-shift and day-shift workers, we would still expect to find the expression of ARNTL anticorrelated with the expression of PER2. The CCD metric evaluates this correlation structure and the ΔCCD metric provides a prima-facie, comparative assessment of clock organization. Our ΔCCD analysis shows evidence of disrupted oscillator organization in IC patients across a wide variety of tissues. Moreover, we find that the homologues of clock outputs that showed strong circadian rhythms across a wide variety of mouse tissues were highly enriched among those transcripts that differed most strongly between the AD and IC groups. This provides additional, suggestive evidence that the change in core clock organization among IC patients has a functional consequence.
Great effort was made in controlling for age, sex, and collection center. Our modified ΔCCD bootstrap procedure is far more stringent as compared to the original ΔCCD procedure that used Combat corrected expression (Table 2). Combat adjusts the expression of individual genes into a “batch-less distribution.” However, it does not address the concern that the correlation between genes may itself be disturbed by a confounding variable. Moreover, it assumes that the temporal sampling of both batches is similar. Our more conservative approach attempts to address these concerns. Regardless of methodology, IC appears to have a clear influence on oscillator organization and a strong influence on clock output genes. Nevertheless, listed confounders may still contribute to the differences we observe. Other confounders, for example, the post-mortem interval between death and sample collection may also contribute (Ferreira et al., 2018). Moreover, as with any retrospective study it is possible that there were unknown, pre-existing differences between the two groups that contribute to the observed differences.
A further limitation of this study is that we cannot conclusively separate the influence of critical illness from the influence of critical care. Critical illness is itself associated with decreased tissue perfusion, hypoxia, and autonomic dysregulation–all of which may potentially impact clock function. While many of the subjects who died on a ventilator had traumatic events or acute medical emergencies (e.g., intracerebral hemorrhage) as a primary cause of death, these emergencies are often coupled with hemodynamic collapse and tissue hypoperfusion. Sepsis and secondary infection are also common. Nonetheless, we can use other data to help evaluate the hypothesis that critical care was itself a contributing factor. Other investigators (Cedernaes et al., 2018; Defour et al., 2020) have studied the influence of sleep loss and fasting–two unfortunate biproducts of current critical care protocols - on transcript expression in human adipose tissues. Notably, the genes identified as being influenced by fasting and sleep deprivation had minimal overlap (Jaccard coefficient = 0.017). Yet both gene sets were heavily enriched among the genes most significantly differentially expressed between the AD and IC groups. We further separated these gene sets into groups that were distinctly up (or down) regulated in the previous sleep deprivation and fasting experiments. These subgroups were again significantly enriched (Table 4). The strong overlap of fasting related genes with those modulated in the IC group is also consistent with our finding that clock function was generally most disrupted in metabolic and digestive tissues.
TABLE 4. Potential influence of sleep loss and fasting in IC group: Previous studies (Cedernaes et al., 2018; Defour et al., 2020) identified genes differentially expressed in human adipose tissue following sleep loss or fasting.
Finally, we must emphasize that our analysis cannot demonstrate any prognostic difference related to circadian dysregulation. All subjects died prior to sample collection. A dedicated study would be required to show the impact of tissue specific clock dysfunction on mortality or other outcomes. However, our study does demonstrate that there is evidence of core clock dysfunction in the tissues of IC patients. Moreover, the transcriptomic changes we observe in the IC group mirror those observed following fasting and sleep deprivation. Taken together, these findings suggest that correcting circadian dysregulation has potential to improve critical care and that improving sleep and regularizing nutrition may advance that goal.
Data availability statement
Publicly available datasets were analyzed in this study. This data can be found here: https://www.gtexportal.org/home/datasets. The scripts used for data analysis, figure generation and the R files required for covariate corrected DCCD analysis are available at https://github.com/ranafi/ICU_clock_disruption.
Author contributions
RA conceived of this research. HH and JF performed the research. RA and HH wrote the manuscript.
Funding
HH was funded by NIA R01 AG068577. JF was funded NHGRI Diversity Action Plan R25HG010323 (MPI Bucan and Kim). RA was supported by NCI R01 CA227485 and NIA R01 AG068577.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2022.942704/full#supplementary-material
References
Aguet F., Barbeira A. N., Bonazzola R., Brown A., Castel S. E., Jo B., et al. (2020). The GTEx Consortium atlas of genetic regulatory effects across human tissues. Science 369, 1318–1330. doi:10.1126/SCIENCE.AAZ1776/SUPPL_FILE/AAZ1776_TABLESS10-S16.XLSX
Aguilar-Arnal L., Sassone-Corsi P. (2013). The circadian epigenome: How metabolism talks to chromatin remodeling. Curr. Opin. Cell Biol. 25, 170–176. doi:10.1016/j.ceb.2013.01.003
Anafi R. C., Lee Y., Sato T. K., Venkataraman A., Ramanathan C., Kavakli I. H., et al. (2014). Machine learning helps identify CHRONO as a circadian clock component. PLoS Biol. 12, e1001840. doi:10.1371/journal.pbio.1001840
Anafi R. C., Pellegrino R., Shockley K. R., Romer M., Tufik S., Pack A. I., et al. (2013). Sleep is not just for the brain: Transcriptional responses to sleep in peripheral tissues. BMC Genomics 14, 362. doi:10.1186/1471-2164-14-362
Archer S. N., Laing E. E., Möller-Levet C. S., van der Veen D. R., Bucca G., Lazar A. S., et al. (2014). Mistimed sleep disrupts circadian regulation of the human transcriptome. Proc. Natl. Acad. Sci. U.S.A. 111, E682–E691. doi:10.1073/pnas.1316335111
Astiz M., Heyde I., Oster H. (2019). Mechanisms of communication in the Mammalian Circadian timing system. Int. J. Mol. Sci. 20, 343. doi:10.3390/ijms20020343
Benjamini Y., Hochberg Y. (1995). Controlling the false discovery rate: A practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B 57, 289. doi:10.1111/j.2517-6161.1995.tb02031.x
Brainard J., Gobel M., Bartels K., Scott B., Koeppen M., Eckle T., et al. (2015). Circadian rhythms in anesthesia and critical care medicine. Semin. Cardiothorac. Vasc. Anesth. 19, 49–60. doi:10.1177/1089253214553066
Brenna A., Olejniczak I., Chavan R., Ripperger J., Langmesser A., Cameroni S., et al. (2019). Cyclin-dependent kinase 5 (CDK5) regulates the circadian clock. eLife 8, e50925. doi:10.7554/eLife.50925
Bugge A., Feng D., Everett L. J., Briggs E. R., Mullican S. E., Wang F., et al. (2012). Rev-erbα and Rev-erbβ coordinately protect the circadian clock and normal metabolic function. Genes Dev. 26, 657–667. doi:10.1101/GAD.186858.112
Cal‐Kayitmazbatir S., Kulkoyluoglu‐Cotul E., Growe J., Selby C. P., Rhoades S. D., Malik D., et al. (2020). CRY1‐CBS binding regulates circadian clock function and metabolism. Febs J. 288, 614–639. doi:10.1111/febs.15360
Cedernaes J., Schönke M., Orzechowski Westholm J., Mi J., Chibalin A., Voisin S., et al. (2018). Acute sleep loss results in tissue-specific alterations in genome-wide DNA methylation state and metabolic fuel utilization in humans. Sci. Adv. 4, eaar8590. doi:10.1126/sciadv.aar8590
Chen E. Y., Tan C. M., Kou Y., Duan Q., Wang Z., Meirelles G. v., et al. (2013). Enrichr: Interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinforma. 14, 128. doi:10.1186/1471-2105-14-128
Chen Y., Lun A. T. L., Smyth G. K. (2016). From reads to genes to pathways: Differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline. F1000Res 5, 1438. doi:10.12688/f1000research.8987.1
Chen Y., McCarthy D., Ritchie M., Robinson M., Smyth G., Hall E. (2022). edgeR: Differential analysis of sequence read count data User’s Guide.
Defour M., Michielsen C. C. J. R., O’Donovan S. D., Afman L. A., Kersten S. (2020). Transcriptomic signature of fasting in human adipose tissue. Physiol. Genomics 52, 451–467. doi:10.1152/physiolgenomics.00083.2020
Diaz E., Diaz I., del Busto C., Escudero D., Pérez S. (2020). Clock genes disruption in the intensive care unit. J. Intensive Care Med. 35, 1497–1504. doi:10.1177/0885066619876572
Etchegaray J.-P., Lee C., Wade P. A., Reppert S. M. (2003). Rhythmic histone acetylation underlies transcription in the mammalian circadian clock. Nature 421, 177–182. doi:10.1038/nature01314
Ferreira P. G., Muñoz-Aguirre M., Reverter F., Sá Godinho C. P., Sousa A., Amadoz A., et al. (2018). The effects of death and post-mortem cold ischemia on human tissue transcriptomes. Nat. Commun. 9, 490. doi:10.1038/s41467-017-02772-x
Freedman N. S., Gazendam J., Levan L., Pack A. I., Schwab R. J. (2001). Abnormal sleep/wake cycles and the effect of environmental noise on sleep disruption in the intensive care unit. Am. J. Respir. Crit. Care Med. 163, 451. doi:10.1164/ajrccm.163.2.9912128
Gabor J. Y., Cooper A. B., Crombach S. A., Lee B., Kadikar N., Bettger H. E., et al. (2003). Contribution of the intensive care unit environment to sleep disruption in mechanically ventilated patients and healthy subjects. Am. J. Respir. Crit. Care Med. 167, 708–715. doi:10.1164/rccm.2201090
Gazendam J. A. C., Van Dongen H. P. A., Grant D. A., Freedman N. S., Zwaveling J. H., Schwab R. J., et al. (2013). Altered circadian rhythmicity in patients in the ICU. Chest 144, 483–489. doi:10.1378/chest.12-2405
Greco C. M., Cervantes M., Fustin J. M., Ito K., Ceglia N., Samad M., et al. (2020). S-adenosyl-l-homocysteine hydrolase links methionine metabolism to the circadian clock and chromatin remodeling. Sci. Adv. 6, 5629–5645. doi:10.1126/SCIADV.ABC5629/SUPPL_FILE/ABC5629_SM.PDF
Haspel J. A., Anafi R., Brown M. K., Cermakian N., Depner C., Desplats P., et al. (2020). Perfect timing: Circadian rhythms, sleep, and immunity - An NIH workshop summary. JCI Insight 5, 131487. doi:10.1172/JCI.INSIGHT.131487
Haspel J., Kim M., Zee P., Schwarzmeier T., Montagnese S., Panda S., et al. (2021). A timely call to arms: COVID-19, the circadian clock, and critical care. J. Biol. Rhythms 36, 55–70. doi:10.1177/0748730421992587
Hirano A., Yumimoto K., Tsunematsu R., Matsumoto M., Oyama M., Kozuka-Hata H., et al. (2013). FBXL21 regulates oscillation of the circadian clock through ubiquitination and stabilization of cryptochromes. Cell 152, 1106–1118. doi:10.1016/j.cell.2013.01.054
Hoekstra M. M. B., Emmenegger Y., Hubbard J., Franken P. (2019). Cold-inducible RNA-binding protein (CIRBP) adjusts clock-gene expression and REM-sleep recovery following sleep deprivation. Elife 8, e43400. doi:10.7554/ELIFE.43400
Jobanputra A. M., Scharf M. T., Androulakis I. P., Sunderram J. (2020). Circadian disruption in critical illness. Front. Neurol. 11, 820. doi:10.3389/fneur.2020.00820
Katada S., Sassone-Corsi P. (2010). The histone methyltransferase MLL1 permits the oscillation of circadian gene expression. Nat. Struct. Mol. Biol. 17, 1414–1421. doi:10.1038/nsmb.1961
Kiss K., Földesi I., Kemény L., Csernus V., Molnár Z., Singer J., et al. (2015). Disturbed circadian rhythm in ICU patients as indicated by melatonin levels: A prospective pilot study. Crit. Care 19, P549. doi:10.1186/cc14629
Koike N., Yoo S.-H., Huang H.-C., Kumar V., Lee C., Kim T.-K., et al. (2012). Transcriptional architecture and chromatin landscape of the core circadian clock in mammals. Science 338, 349–354. doi:10.1126/science.1226339
Korotkevich G., Sukhov V., Budin N., Shpak B., Artyomov M. N., Sergushichev A. (2021). Fast gene set enrichment analysis. bioRxiv. doi:10.1101/060012
Kuleshov M. v., Jones M. R., Rouillard A. D., Fernandez N. F., Duan Q., Wang Z., et al. (2016). Enrichr: A comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 44, W90–W97. doi:10.1093/NAR/GKW377
Kume K., Zylka M. J., Sriram S., Shearman L. P., Weaver D. R., Jin X., et al. (1999). mCRY1 and mCRY2 are essential components of the negative limb of the circadian clock feedback loop. Cell 98, 193–205. doi:10.1016/S0092-8674(00)81014-4
Leek J. T., Johnson W. E., Parker H. S., Fertig E. J., Jaffe A. E., Zhang Y., et al. (2021). sva: Surrogate Variable Analysis.
Li Y., Shan Y., Kilaru G., Berto S., Wang G. Z., Cox K. H., et al. (2020). Epigenetic inheritance of circadian period in clonal cells. Elife 9, e54186. doi:10.7554/ELIFE.54186
Lowrey P. L., Shimomura K., Antoch M. P., Yamazaki S., Zemenides P. D., Ralph M. R., et al. (2000). Positional syntenic cloning and functional characterization of the mammalian circadian mutation tau. Science 288, 483–491. doi:10.1126/science.288.5465.483
Maas M. B., Iwanaszko M., Lizza B. D., Reid K. J., Braun R. I., Zee P. C., et al. (2020a). Circadian gene expression rhythms during critical illness. Crit. Care Med. 48, e1294–e1299. doi:10.1097/CCM.0000000000004697
Maas M. B., Lizza B. D., Abbott S. M., Liotta E. M., Gendy M., Eed J., et al. (2020b). Factors disrupting melatonin secretion rhythms during critical illness. Crit. Care Med. 48, 854–861. doi:10.1097/CCM.0000000000004333
Milev N. B., Reddy A. B. (2015). Circadian redox oscillations and metabolism. Trends Endocrinol. Metabolism 26, 430–437. doi:10.1016/j.tem.2015.05.012
Möller-Levet C. S., Archer S. N., Bucca G., Laing E. E., Slak A., Kabiljo R., et al. (2013). Effects of insufficient sleep on circadian rhythmicity and expression amplitude of the human blood transcriptome. Proc. Natl. Acad. Sci. U.S.A. 110, E1132–E1141. doi:10.1073/pnas.1217154110
Mundigler G., Delle-Karth G., Koreny M., Zehetgruber M., Steindl-Munda P., Marktl W., et al. (2002). Impaired circadian rhythm of melatonin secretion in sedated critically ill patients with severe sepsis*. Crit. Care Med. 30, 536–540. doi:10.1097/00003246-200203000-00007
Nakahata Y., Kaluzova M., Grimaldi B., Sahar S., Hirayama J., Chen D., et al. (2008). The NAD+-dependent deacetylase SIRT1 modulates CLOCK-mediated chromatin remodeling and circadian control. Cell 134, 329–340. doi:10.1016/j.cell.2008.07.002
Nicolaides N. C., Charmandari E., Chrousos G. P., Kino T. (2014). Circadian endocrine rhythms: The hypothalamic-pituitary-adrenal axis and its actions. Ann. N.Y. Acad. Sci. 1318, 71–80. doi:10.1111/nyas.12464
Nicolaides N. C., Charmandari E., Kino T., Chrousos G. P. (2017). Stress-related and circadian secretion and target tissue actions of glucocorticoids: Impact on health. Front. Endocrinol. 8, 70. doi:10.3389/fendo.2017.00070
Patke A., Young M. W., Axelrod S. (2019). Molecular mechanisms and physiological importance of circadian rhythms. Nat. Rev. Mol. Cell Biol. 21, 67–84. doi:10.1038/s41580-019-0179-2
Phipson B., Smyth G. K. (2010). Permutation p-values should never be zero: calculating exact p-values when permutations are randomly drawn. Stat. Appl. Genet. Mol. Biol. 9, 39.
Preitner N., Damiola F., Luis-Lopez-Molina L., Zakany J., Duboule D., Albrecht U., et al. (2002). The orphan nuclear receptor REV-erbα controls circadian transcription within the positive limb of the mammalian circadian oscillator. Cell 110, 251–260. doi:10.1016/S0092-8674(02)00825-5
Robinson M. D., McCarthy D. J., Smyth G. K. (2010). edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi:10.1093/BIOINFORMATICS/BTP616
Ruben M. D., Wu G., Smith D. F., Schmidt R. E., Francey L. J., Lee Y. Y., et al. (2018). A database of tissue-specific rhythmically expressed human genes has potential applications in circadian medicine. Sci. Transl. Med. 10, eaat8806. doi:10.1126/scitranslmed.aat8806
Shilts J., Chen G., Hughey J. J. (2018). Evidence for widespread dysregulation of circadian clock progression in human cancer. PeerJ 6, e4327. doi:10.7717/peerj.4327
Smith G., Nielsen M. (1999). ABC of intensive care: Criteria for admission. BMJ 318, 1544–1547. doi:10.1136/BMJ.318.7197.1544
Somekh J., Lotan N., Sussman E., Yehuda G. (2022). Predicting mechanical ventilation effects on six human tissue transcriptomes. PLoS ONE 17, e0264919. doi:10.1371/journal.pone.0264919
Takahashi J. S., Hong H.-K., Ko C. H., McDearmon E. L. (2008). The genetics of mammalian circadian order and disorder: Implications for physiology and disease. Nat. Rev. Genet. 9, 764–775. doi:10.1038/nrg2430
Takahashi J. S. (2016). Transcriptional architecture of the mammalian circadian clock. Nat. Rev. Genet. 18, 164–179. doi:10.1038/nrg.2016.150
Telias I., Wilcox M. E. (2019). Sleep and circadian rhythm in critical illness. Crit. Care 23, 82. doi:10.1186/s13054-019-2366-0
Vollmers C., Gill S., DiTacchio L., Pulivarthy S. R., Le H. D., Panda S., et al. (2009). Time of feeding and the intrinsic circadian clock drive rhythms in hepatic gene expression. Proc. Natl. Acad. Sci. U.S.A. 106, 21453–21458. doi:10.1073/pnas.0909591106
Wang H., van Spyk E., Liu Q., Geyfman M., Salmans M. L., Kumar V., et al. (2017). Time-restricted feeding shifts the skin circadian clock and alters UVB-induced DNA damage. Cell Rep. 20, 1061–1072. doi:10.1016/J.CELREP.2017.07.022
Yoo S.-H., Mohawk J. A., Siepka S. M., Shan Y., Huh S. K., Hong H.-K., et al. (2013). Competing E3 ubiquitin ligases govern circadian periodicity by degradation of CRY in nucleus and cytoplasm. Cell 152, 1091–1105. doi:10.1016/j.cell.2013.01.055
Zhang R., Lahens N. F., Ballance H. I., Hughes M. E., Hogenesch J. B. (2014). A circadian gene expression atlas in mammals: Implications for biology and medicine. Proc. Natl. Acad. Sci. U.S.A. 111, 16219–16224. doi:10.1073/pnas.1408886111
Keywords: circadian rhythm, circadian disruption, transcription, gene correlation analysis, critical care, critical illness, ICU, intensive care unit
Citation: Hollis HC, Francis JN and Anafi RC (2022) Multi-tissue transcriptional changes and core circadian clock disruption following intensive care. Front. Physiol. 13:942704. doi: 10.3389/fphys.2022.942704
Received: 12 May 2022; Accepted: 28 June 2022;
Published: 15 August 2022.
Edited by:
Roee Gutman, Migal—Galilee Research Institute, IsraelReviewed by:
Bala S. C. Koritala, Cincinnati Children’s Hospital Medical Center, United StatesAndrea Brenna, Université de Fribourg, Switzerland
Copyright © 2022 Hollis, Francis and Anafi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Ron C. Anafi, ron.anafi@pennmedicine.upenn.edu